How these scenario elements are represented and which values are
chosen depends on what one would like to achieve. Simulating the growth
of Lemna in a controlled lab environment will likely require
different inputs than Lemna growing in an outdoor water body,
for example.
Core functions
To make functions and sample datasets of the lemna package
available in your R workspace, load the library first:
The package function param_defaults() provides a list
with all suggested default parameters. Some parameter values will be
missing, i.e. set to NA, because they are substance
specific and default values would not be meaningful for these:
# get list of default parameters
params <- param_defaults()
params$k_photo_max
#> [1] 0.47
params$EC50_int # substance specific
#> [1] NA
# get default parameters and set a custom parameter value
myparam <- param_defaults(c(EC50_int = 42))
myparam$EC50_int
#> [1] 42
The growth of a Lemna population is simulated using the
lemna() function. The required scenario data are either
supplied individually on function call or are passed as a pre-defined
scenario object, such as the metsulfuron sample
scenario:
lemna(metsulfuron)
#> time BM M_int C_int FrondNo
#> 1 0 0.001200000 0.0000000000 0.000000000 12.00000
#> 2 1 0.001725634 0.0063939067 0.221871311 17.25634
#> 3 2 0.002003771 0.0120701443 0.360701550 20.03771
#> 4 3 0.002171628 0.0164518190 0.453640790 21.71628
#> 5 4 0.002320205 0.0198229397 0.511593653 23.20205
#> 6 5 0.002467715 0.0225374005 0.546880329 24.67715
#> 7 6 0.002619619 0.0248568637 0.568187670 26.19619
#> 8 7 0.002778289 0.0269567504 0.580996626 27.78289
#> 9 8 0.002999069 0.0167083478 0.333603522 29.99069
#> 10 9 0.003711271 0.0103269725 0.166622594 37.11271
#> 11 10 0.005507806 0.0063830367 0.069395642 55.07806
#> 12 11 0.008377291 0.0039454957 0.028202105 83.77291
#> 13 12 0.012750953 0.0024390409 0.011454074 127.50953
#> 14 13 0.019407273 0.0015074611 0.004651201 194.07273
#> 15 14 0.029537211 0.0009316238 0.001888664 295.37211
lemna() returns a table which describes the change of
state variables over time. In addition, some supporting derived
variables such as internal toxicant concentration (C_int)
and the number of fronds (FrondNo) will be returned by
default.
A visual description of the simulated scenario and its results can be
created by running the plot() function. The
plot() function requires a simulation result as its first
argument:
The effect of the toxicant on the Lemna population can be
calculated using the effect() function. It requires
scenario data the same way as lemna() does. For the sample
metsulfuron scenario, the effects of the toxicant are as
follows:
effect(metsulfuron)
#> BM r
#> 93.12935 45.53310
In this scenario, exposure to the toxicant resulted in an 93%
decrease of population size (BM) and a 46% decrease in
average growth rate (r) until the end of the simulation.
Effects are always calculated relative to an identical control scenario
which contains no toxicant exposure.
For more information on the metsulfuron sample scenario,
please refer to the help files:
Tutorial
Simulate the Lemna growth model
To simulate a Lemna population, one has to pass the four
mandatory scenario elements to the lemna() function:
# initial state of the model system: 1.0 g dw biomass, 0.0 ug/m2 internal toxicant
myinit <- c(BM=1, M_int=0)
# simulated period and output time points: each day for 7 days
mytimes <- 0:7
# default model parameters + substance specific values
myparam <- param_defaults(c(
EC50_int = 4.16,
b = 0.3,
P = 0.0054
))
# constant environmental conditions, including exposure
myenvir <- list(
tmp = 18, # 18 °C ambient temperature
irr = 15000, # 15,000 kJ m-2 d-1 irradiance
P = 0.3, # 0.3 mg L-1 Phosphorus concentration
N = 0.6, # 0.6 mg L-1 Nitrogen concentration
conc = 1 # 1 ug/L toxicant concentration
)
lemna(
init = myinit,
times = mytimes,
param = myparam,
envir = myenvir
)
#> time BM M_int C_int FrondNo
#> 1 0 1.000000 0.000000 0.0000000 10000.00
#> 2 1 1.198841 5.028960 0.2511888 11988.41
#> 3 2 1.410120 9.480796 0.4025985 14101.20
#> 4 3 1.648541 13.628139 0.4950173 16485.41
#> 5 4 1.921189 17.702777 0.5517657 19211.89
#> 6 5 2.234682 21.897913 0.5867735 22346.82
#> 7 6 2.596018 26.379982 0.6084857 25960.18
#> 8 7 3.012897 31.299206 0.6220605 30128.97
The init argument controls at which system state the
simulation starts. The times argument defines the length of
the simulated period and for which time points results are returned. The
temporal resolution of results can be increased by specifying additional
output times:
simresult <- lemna(
init = myinit,
times = seq(0, 7, 0.1), # a step length of 0.1 days = ~2 hours
param = myparam,
envir = myenvir
)
tail(simresult)
#> time BM M_int C_int FrondNo
#> 66 6.5 2.797005 28.77600 0.6160568 27970.05
#> 67 6.6 2.838957 29.26990 0.6173707 28389.57
#> 68 6.7 2.881514 29.76904 0.6186253 28815.14
#> 69 6.8 2.924684 30.27355 0.6198233 29246.84
#> 70 6.9 2.968477 30.78358 0.6209675 29684.77
#> 71 7.0 3.012902 31.29926 0.6220605 30129.02
The resulting table now contains ten times as much rows because we
decreased the step length by a factor of ten but simulated the same
period, i.e. seven days. It can be observed that the state-variables
differ slightly at the end of the simulation although the scenarios were
otherwise identical. The differences originate from small numerical
errors introduced by the solver of the model’s Ordinary Differential
Equations (ODE). The step-length in time can have influence on the
precision of simulation results. To decrease the solver’s step length
without increasing the number of result time points, make use of the
optional argument hmax. The smaller hmax, the
more precise the results:
# hmax=0.01 forces a maximum step length of 0.01 days = ~15 minutes
lemna(myinit, mytimes, myparam, myenvir, hmax = 0.01)
#> time BM M_int C_int FrondNo
#> 1 0 1.000000 0.000000 0.0000000 10000.00
#> 2 1 1.198845 5.028969 0.2511883 11988.45
#> 3 2 1.410128 9.480842 0.4025982 14101.28
#> 4 3 1.648549 13.628208 0.4950173 16485.49
#> 5 4 1.921199 17.702865 0.5517656 19211.99
#> 6 5 2.234694 21.898028 0.5867734 22346.94
#> 7 6 2.596031 26.380122 0.6084857 25960.31
#> 8 7 3.012913 31.299374 0.6220605 30129.13
By default, simulation results contain supporting variables such as
internal toxicant concentration and total frond number. These are
calculated from simulation results and model parameters for reasons of
convenience. If these variables are not required, they can be disabled
by setting the optional argument nout = 0:
lemna(myinit, mytimes, myparam, myenvir, nout = 0)
#> time BM M_int
#> 1 0 1.000000 0.000000
#> 2 1 1.198841 5.028960
#> 3 2 1.410120 9.480796
#> 4 3 1.648541 13.628139
#> 5 4 1.921189 17.702777
#> 6 5 2.234682 21.897913
#> 7 6 2.596018 26.379982
#> 8 7 3.012897 31.299206
Using environmental time-series
The previous examples mostly assumed that environmental variables
stay constant in time. To simulate a scenario with changing
environmental variables, such as a temperature curve or exposure
pattern, one has to define or load a data time-series. The model accepts
time-series for all environmental variables, i.e. exposure
concentration, temperature, irradiation, phosphorus concentration, and
nitrogen concentration.
Within the scope of this package, time-series are represented by a
data.frame containing exactly two numerical columns: the
first column for time, the second for the variable’s value. The column
names are irrelevant but sensible names may help documenting the data.
As an example, the metsulfuron sample scenario contains a
step-function as its exposure time-series: seven days of 1 ug/L
metsulfuron-methyl starting at time point zero
(0.0), followed by seven days of recovery (no
exposure).
metsulfuron$envir$conc
#> time conc
#> 1 0.00 1
#> 2 7.00 1
#> 3 7.01 0
#> 4 14.00 0
Time points of the time-series and time points processed by the ODE
solver may not always match. To derive environmental variable values
which are not explicitly part of the time-series, variable values are
interpolated with a linear function. If the time-series does not cover
the full simulation period, the closest value from the time-series is
used. In the case of the metsulfuron sample scenario, the
step function will effectively extend to infinity, i.e. any time point
before day 7.0 will have 1 ug/L of exposure and any time
point after 7.01 will have no exposure.
As an example, we will modify the metsulfuron sample
scenario to use an exposure time-series that declines linearly between
start and day seven:
# define start and end points for the exposure series, the values
# in between will be interpolated
myexpo <- data.frame(time=c(0, 7), conc=c(1, 0))
# modify the sample scenario's exposure series
myenvir <- metsulfuron$envir
myenvir$conc <- myexpo
# simulate the sample scenario with modified environmental variables
plot(lemna(metsulfuron, envir=myenvir))
Time-series and data.frame objects can be stored
conveniently as .csv files which can be created and edited
by common spreadsheet programs such as Microsoft Excel. Be
aware that the separator character used by R and your
spreadsheet program may differ depending on your computer’s locale
settings.
set.seed(23)
# define a random time-series, values will be uniformly distributed between
# the values 0.1 and 3.0, e.g to represent an exposure time-series
myexpo <- data.frame(time = 0:14,
conc = round(runif(15, 0.1, 3.0), 1))
# plot the time-series
plot(myexpo, main="Random exposure time-series")
lines(myexpo)
# write data to .csv file in working directory
write.csv(myexpo, file="random_series.csv", row.names=FALSE)
# write data using semicolons as separating character
write.csv2(myexpo, file="random_series2.csv", row.names=FALSE)
# read file from working directory
myimport <- read.csv(file="random_series.csv")
# check that written and read data are identical
myexpo$conc == myimport$conc
#> [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
Time-series can be imported manually as in the previous example or
they can be imported automatically by the lemna() function
for convenience. If an environmental variable is set to a string, it
will be interpreted as a file path and lemna() will try to
import the time-series using read.csv():
# automatically load the exposure time-series from a file
myenvir <- metsulfuron$envir
myenvir$conc <- "random_series.csv"
# simulate the sample scenario with the exposure series loaded from a .csv file
plot(lemna(metsulfuron, envir=myenvir), legend=FALSE)
For a more complex scenario that uses hourly and daily time-series of
exposure and temperature/irradiance, respectively, please have a look at
e.g. the focusd1 scenario:
myenvir <- focusd1$envir
myenvir$conc
myenvir$tmp
myenvir$irr
Using simulation results
Simulation results are returned as a table, i.e. a
data.frame object. The table will contain the state
variables biomass (BM) and internal toxicant mass
(M_int) for each requested output time point. The table may
also contain additional columns for other supporting variables. The data
can be processed like any other dataset in R to e.g. create
plots, derive other values, or to perform statistical tests:
myresult <- lemna(focusd1)
head(myresult)
#> time BM M_int C_int FrondNo
#> 1 0 80.00000 0.0000000 0.0000000000 200000.0
#> 2 1 79.48177 0.6614503 0.0004983256 198704.4
#> 3 2 79.00414 2.2061048 0.0016720907 197510.3
#> 4 3 78.57538 4.1750616 0.0031817050 196438.5
#> 5 4 78.23362 6.2899397 0.0048143381 195584.0
#> 6 5 78.01035 8.3810150 0.0064332127 195025.9
To get an initial impression of a scenario and its results, simply
pass the simulation result to the plot() function:
As an example, we will analyze if and how the internal toxicant
concentration (C_int) correlates with the internal toxicant
mass (M_int):
summary(lm(C_int ~ M_int, myresult))
#>
#> Call:
#> lm(formula = C_int ~ M_int, data = myresult)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -0.0128702 -0.0017575 -0.0000844 0.0023533 0.0112755
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) -3.071e-04 3.301e-04 -0.93 0.353
#> M_int 7.750e-04 1.878e-05 41.27 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.004307 on 364 degrees of freedom
#> Multiple R-squared: 0.8239, Adjusted R-squared: 0.8234
#> F-statistic: 1703 on 1 and 364 DF, p-value: < 2.2e-16
The linear model indicates a strong correlation of internal toxicant
mass and concentration which intuitively makes sense. The correlation is
not a 100% because biomass is a confounding factor in the model
equations.
Derive effect endpoints
To quantify the influence a toxicant exerts on a Lemna
population, use the effect() function. It works similar to
lemna() and accepts the same arguments in order to specify
a scenario:
# calculate effects on biomass in sample scenario
effect(metsulfuron)
#> BM r
#> 93.12935 45.53310
The return values describe the effect in percent (%) on the
respective effect endpoint. Effects are calculated relative to a control
scenario which exhibits no exposure. By default, the effect refers to
the reduction in biomass (BM) or average growth rate
(r) at the end of the simulation. In the example above,
biomass was reduced by 93% and the growth rate was reduced by 46% in the
Lemna population due to exposure to the toxicant.
If a scenario covers a long time period but effects are desired for
an earlier time point, the scenario can be cut short by using the
duration argument. If duration is set, the
scenario will be clipped to the time period from t0 to
t0 + duration:
# calculate effects on biomass after 7 days, instead of 14
effect(metsulfuron, duration=7)
#> BM r
#> 87.77416 71.45610
In this example, the effect on biomass is smaller after 7 days
compared to the effects after 14 days. However, the average growth rate
experienced a strong decrease from 46 to 71%.
Create scenario objects
A Lemna growth scenario consists of the following four
mandatory scenario elements: model parameters, environmental variables,
initial state, and output times. The elements can be passed to
lemna() and effect() separately or they can be
combined to a compact scenario object. All sample scenarios which were
used in this tutorial are scenario objects:
# list properties of the sample scenario object
metsulfuron
Scenario objects are basically just a base R
list object with some additional metadata. If correctly
defined, scenario objects fully describe a scenario and can be passed to
e.g. lemna() without additional arguments. It is, however,
possible to override a scenario object’s data by passing an alternative
dataset:
# custom output times and time period:
# four days with a 12 hour time step
mytimes <- seq(0, 4, 0.5)
# simulate sample scenario with custom output times & period
lemna(metsulfuron, times=mytimes)
#> time BM M_int C_int FrondNo
#> 1 0.0 0.001200000 0.000000000 0.0000000 12.00000
#> 2 0.5 0.001473874 0.003215291 0.1306301 14.73874
#> 3 1.0 0.001725627 0.006393890 0.2218716 17.25627
#> 4 1.5 0.001891853 0.009391186 0.2972463 18.91853
#> 5 2.0 0.002003761 0.012070096 0.3607018 20.03761
#> 6 2.5 0.002092620 0.014411810 0.4123934 20.92620
#> 7 3.0 0.002171612 0.016451729 0.4536416 21.71612
#> 8 3.5 0.002246676 0.018240740 0.4861670 22.46676
#> 9 4.0 0.002320186 0.019822776 0.5115938 23.20186
A custom scenario object can be created by passing the scenario
elements to new_lemna_scenario():
myscenario <- new_lemna_scenario(
init = c(BM=1, M_int=0),
times = 0:7,
param = param_defaults(c(EC50_int = 4.16, b = 0.3, P = 0.0054)),
envir = list(
tmp = 18, # 18 °C ambient temperature
irr = 15000, # 15,000 kJ m-2 d-1 irradiance
P = 0.3, # 0.3 mg L-1 Phosphorus concentration
N = 0.6, # 0.6 mg L-1 Nitrogen concentration
conc = 1 # 1 ug/L toxicant concentration
)
)
lemna(myscenario)
#> time BM M_int C_int FrondNo
#> 1 0 1.000000 0.000000 0.0000000 10000.00
#> 2 1 1.198841 5.028960 0.2511888 11988.41
#> 3 2 1.410120 9.480796 0.4025985 14101.20
#> 4 3 1.648541 13.628139 0.4950173 16485.41
#> 5 4 1.921189 17.702777 0.5517657 19211.89
#> 6 5 2.234682 21.897913 0.5867735 22346.82
#> 7 6 2.596018 26.379982 0.6084857 25960.18
#> 8 7 3.012897 31.299206 0.6220605 30128.97
Speed up simulations with compiled code
The Lemna growth model is simulated by default using model
equations implemented in pure R. In case many simulations have
to be conducted or the time required to get results becomes an issue,
the compiled code feature can be used. The lemna package
provides an alternative implementation of the Klein et al.
model equations using C code. The C code executes
significantly faster than the pure R alternative.
# use model implemented in pure R
tail(lemna(metsulfuron, ode_mode="r"), n = 1)
#> time BM M_int C_int FrondNo
#> 15 14 0.02953721 0.0009316238 0.001888664 295.3721
# use model implemented in C
tail(lemna(metsulfuron, ode_mode="c"), n = 1)
#> time BM M_int C_int FrondNo
#> 15 14 0.02953721 0.0009316238 0.001888664 295.3721
Simulation results of R and C code will be
identical as far as numerical precision allows. The speed increase of
using C will range from a factor of 3 to 5 for short scenarios
and up to 50+ for longer scenarios:
# Benchmark the shorter metsulfuron scenario
microbenchmark::microbenchmark(
lemna(metsulfuron, ode_mode="r"),
lemna(metsulfuron, ode_mode="c")
)
#> Unit: milliseconds
#> expr min lq mean median uq
#> lemna(metsulfuron, ode_mode = "r") 5.9443 6.03445 6.261556 6.08350 6.14925
#> lemna(metsulfuron, ode_mode = "c") 1.6148 1.66900 1.787848 1.71915 1.84595
#> max neval
#> 9.8552 100
#> 4.8319 100
# Benchmark the more complex and longer focusd1 scenario
microbenchmark::microbenchmark(
lemna(focusd1, ode_mode="r"),
lemna(focusd1, ode_mode="c"),
times = 10
)
#> Unit: milliseconds
#> expr min lq mean median
#> lemna(focusd1, ode_mode = "r") 1652.9946 1662.8406 1683.02801 1678.3198
#> lemna(focusd1, ode_mode = "c") 34.4333 34.8228 35.84511 35.1424
#> uq max neval
#> 1699.6687 1733.9992 10
#> 37.2548 38.7207 10
There is however a small disadvantage to using the C model:
if there are any issues stemming from, for example, invalid parameters,
the error messages raised by the C code might be less
descriptive than those from R. On the other hand, the
C code can output on demand almost all intermediary model
variables which can support debugging and model understanding:
# simulate and request all additional output variables
lemna(metsulfuron, ode_mode="c", nout=18)
#> time BM M_int C_int FrondNo f_loss f_photo
#> 1 0 0.001200000 0.0000000000 0.000000000 12.00000 1 1.0000000
#> 2 1 0.001725634 0.0063939067 0.221871311 17.25634 1 0.6194142
#> 3 2 0.002003771 0.0120701443 0.360701550 20.03771 1 0.3125173
#> 4 3 0.002171628 0.0164518190 0.453640790 21.71628 1 0.2562340
#> 5 4 0.002320205 0.0198229397 0.511593653 23.20205 1 0.2409028
#> 6 5 0.002467715 0.0225374005 0.546880329 24.67715 1 0.2350152
#> 7 6 0.002619619 0.0248568637 0.568187670 26.19619 1 0.2322779
#> 8 7 0.002778289 0.0269567504 0.580996626 27.78289 1 0.2308635
#> 9 8 0.002999069 0.0167083478 0.333603522 29.99069 1 0.3435404
#> 10 9 0.003711271 0.0103269725 0.166622594 37.11271 1 0.8253400
#> 11 10 0.005507806 0.0063830367 0.069395642 55.07806 1 0.9941664
#> 12 11 0.008377291 0.0039454957 0.028202105 83.77291 1 0.9998612
#> 13 12 0.012750953 0.0024390409 0.011454074 127.50953 1 0.9999967
#> 14 13 0.019407273 0.0015074611 0.004651201 194.07273 1 0.9999999
#> 15 14 0.029537211 0.0009316238 0.001888664 295.37211 1 1.0000000
#> fT_photo fI_photo fP_photo fN_photo fBM_photo fCint_photo C_int_unb
#> 1 0.2410198 0.775 0.9858692 0.9463722 0.9999932 1.0000000 0.000000000
#> 2 0.2410198 0.775 0.9858692 0.9463722 0.9999902 0.6194142 0.295828415
#> 3 0.2410198 0.775 0.9858692 0.9463722 0.9999886 0.3125173 0.480935400
#> 4 0.2410198 0.775 0.9858692 0.9463722 0.9999877 0.2562340 0.604854387
#> 5 0.2410198 0.775 0.9858692 0.9463722 0.9999868 0.2409028 0.682124871
#> 6 0.2410198 0.775 0.9858692 0.9463722 0.9999860 0.2350152 0.729173772
#> 7 0.2410198 0.775 0.9858692 0.9463722 0.9999851 0.2322779 0.757583560
#> 8 0.2410198 0.775 0.9858692 0.9463722 0.9999842 0.2308635 0.774662168
#> 9 0.2410198 0.775 0.9858692 0.9463722 0.9999830 0.3435404 0.444804696
#> 10 0.2410198 0.775 0.9858692 0.9463722 0.9999789 0.8253400 0.222163459
#> 11 0.2410198 0.775 0.9858692 0.9463722 0.9999687 0.9941664 0.092527522
#> 12 0.2410198 0.775 0.9858692 0.9463722 0.9999524 0.9998612 0.037602806
#> 13 0.2410198 0.775 0.9858692 0.9463722 0.9999276 0.9999967 0.015272099
#> 14 0.2410198 0.775 0.9858692 0.9463722 0.9998897 0.9999999 0.006201602
#> 15 0.2410198 0.775 0.9858692 0.9463722 0.9998322 1.0000000 0.002518218
#> C_ext Tmp Irr Phs Ntr dBM dM_int
#> 1 1 12 15000 0.3 0.6 0.0005040000 0.0064800000
#> 2 1 12 15000 0.3 0.6 0.0004160929 0.0062420735
#> 3 1 12 15000 0.3 0.6 0.0001941316 0.0050129590
#> 4 1 12 15000 0.3 0.6 0.0001529477 0.0038111983
#> 5 1 12 15000 0.3 0.6 0.0001466934 0.0029915453
#> 6 1 12 15000 0.3 0.6 0.0001491910 0.0024820689
#> 7 1 12 15000 0.3 0.6 0.0001550044 0.0021863653
#> 8 1 12 15000 0.3 0.6 0.0001625461 0.0020328523
#> 9 0 12 15000 0.3 0.6 0.0003342882 -0.0080390165
#> 10 0 12 15000 0.3 0.6 0.0012540748 -0.0049686961
#> 11 0 12 15000 0.3 0.6 0.0022981775 -0.0030711198
#> 12 0 12 15000 0.3 0.6 0.0035179156 -0.0018983268
#> 13 0 12 15000 0.3 0.6 0.0053553805 -0.0011735146
#> 14 0 12 15000 0.3 0.6 0.0081510539 -0.0007252964
#> 15 0 12 15000 0.3 0.6 0.0124056285 -0.0004482394