lemna {lemna}R Documentation

Simulate a Lemna scenario

Description

The function numerically integrates the Lemna model using the supplied parameters, environmental factor time series, output time points, and other settings. Numerical integration is handled by deSolve::lsoda() by default which is well suited to handle stiff and non-stiff numerical problems.

Usage

lemna(...)

## Default S3 method:
lemna(
  init = c(BM = 0, M_int = 0),
  times,
  param,
  envir,
  ode_mode = c("r", "c"),
  nout = 2,
  ...
)

## S3 method for class 'lemna_scenario'
lemna(x, init, times, param, envir, ...)

Arguments

...

additional parameters passed on to deSolve::ode()

init

initial state of the model variables

times

numeric vector, output times for which model results are returned

param

named list, Lemna model parameters

envir

named list, contains time series data for each of the five environmental variables

ode_mode

character, switch controls if ODE functions implemented in R or C are used, defaults to R

nout

numeric, controls the number of additional output variables of the model, such as C_int (internal concentration) or FrondNo (the number of fronds), see e.g. deSolve::lsoda() for details. Defaults to 2

x

a lemna_scenario object

Value

data.frame with simulation results

Methods (by class)

State variables

The model has two state variables:

Output times

The function will only return results for the requested time points. The solver may make additional time steps which are not returned, depending on the solver settings such as numerical tolerance. The number of output times and the distance between them can have influence on the numerical precision of results. Please refer to the manual of the deSolve package for more information on that topic.

Parameters

The list of available model parameters, their units, and suggested default values is documented in param_defaults().

Environmental variables

The model requires a time series for each of the five environmental variables For ease of use, a time series can be represented by either a constant numeric value, a data.frame containing a time series, a character string with a path to a CSV file containing a time series, or a function returning a value for a certain time point.

The five environmental variables are as follows:

A data.frame containing a time series must consist of exactly two columns: The first column contains the time in days, the second column the environmental factor in question. The data.frame must contain at least one row. Time points are not required to coincide with the solver's output times but the user should take care that the solver's time step length is sufficiently small so that it can consider the relevant dynamics of the time series. This problem can be avoided by calling the model ODEs implemented in C, see argument ode_mode.

If a string is passed as an environmental factor, the function will interpret the string as a file path and will try to load the contents of a CSV file. The same requirements as for a data.frame apply to a time series loaded from a file.

When using the pure R code ODE, it is also possible to pass a function as an environmental factor. The function must accept exactly one argument which represents time and must return a single scalar value. Using a function in combination with the compiled code solver will raise an error.

R vs. compiled code

The model can be simulated using pure R code (the default) or an implementation that uses the compiled code feature of deSolve. Compiled code is almost always significantly faster than pure R. The solver for compiled ODEs also handles environmental variables better than the pure R version. For optimal performance, the time series of environmental variables should always be as short as possible and as long as needed.

To use the compiled code feature, set the argument ode_mode = "c".

Simulation output

For reasons of convenience, the return value contains by default two additional variables derived from simulation results: the internal concentration C_int as well as the number of fronds FrondNo. These can be disabled by setting the argument nout = 0.

The compiled code model can return up to 18 additional variables which represent intermediary variables, environmental variables, and derivatives calculated within the model equations. Please refer to the model description of Klein et al. for more information on these variables and how they influence model behavior.

The available output levels are as follows:

References

Klein J., Cedergreen N., Heine S., Kehrein N., Reichenberger S., Rendal C., Schmitt W., Hommen U., 2022: Refined description of the Lemna TKTD growth model based on Schmitt et al. (2013) – equation system and default parameters, implementation in R. Report of the working group Lemna of the SETAC Europe Interest Group Effect Modeling. Version 1.1, uploaded on 09 May 2022. https://www.setac.org/group/SEIGEffectModeling

Examples

# Simulate the metsulfuron example scenario
lemna(metsulfuron)

# Create a simple plot of the number of fronds
lemna(metsulfuron) -> result
plot(result$time, result$FrondNo)

# Create a nicer plot using a dedicated plotting routine
plot(result)

# Simulate the example scenario for a period of 30 days
lemna(metsulfuron, times=0:30) -> result30
plot(result30)

##
## Create a custom Lemna scenario from scratch
##

# Initial state: 12 fronds, no toxicant mass
myinit <- c(BM=0.0012, M_int=0)

# Output times: simulate 7 days with a ~2 hour time step
mytime <- seq(0, 7, 0.1)

# Default model parameters + TKTD parameters of a hypothetical substance
myparam <- param_defaults(list(
  EC50_int = 0.1, # 50% effect level (ug L-1)
  b = 0.7,        # slope parameter (-)
  P = 0.01        # permeability (cm d-1)
))

# Custom environmental variables
myenvir <- list(
  # exposure step function:
  # 3 days no exposure, followed by 4 days of 10 ug L-1
  conc = data.frame(t=c(0,3,4,7), conc=c(0,0,10,10)),
  tmp = 18,    # constant temperature of 18 °C
  irr = 15000, # constant irradiance of 15,000 kJ m-2 d-1
  N = 0.6,     # constant Nitrogen concentration of 0.6 mg L-1
  P = 0.3      # constant Phosphorus concentration of 0.3 mg L-1
)

# Simulate the custom scenario and plot results
lemna(init=myinit, times=mytime, param=myparam, envir=myenvir) -> result_custom
plot(result_custom)

# Simulate again, forcing the solver to use smaller time steps of hmax=0.001.
# The resulting curves are almost identical for this example.
lemna(init=myinit, times=mytime, param=myparam, envir=myenvir, hmax=0.001) -> result_custom2
library(ggplot2)
ggplot(result_custom, aes(time, FrondNo)) +
  geom_line() +
  geom_line(data=result_custom2, color="red", style="dashed")

# Combine all settings into a scenario object and simulate it
scenario <- new_lemna_scenario(
 init = myinit,
 param = myparam,
 times = mytime,
 envir = myenvir
)
lemna(scenario)
[Package lemna version 1.0.1 Index]