From f656e88dc9973dbfa1527157840774100a8f9da4 Mon Sep 17 00:00:00 2001 From: Nils Kehrein Date: Tue, 10 May 2022 10:15:08 +0200 Subject: [PATCH] References updated to report v1.1, package version set to v1.0.0 --- DESCRIPTION | 15 +- LICENSE | 2 +- LICENSE.md | 2 +- NEWS.md | 11 + R/lemna.R | 9 +- R/param.R | 11 +- README.Rmd | 14 +- README.md | 19 +- dev/release.R | 2 +- docs/lemna-introduction.html | 463 ++++++++++++++++++++++---- docs/lemna-verification.html | 398 +++++++++++++++++++--- man/lemna.Rd | 9 +- man/param_defaults.Rd | 11 +- src/lemna.c | 2 +- tests/testthat/test-verify-variable.R | 8 +- vignettes/lemna-introduction.Rmd | 11 +- vignettes/lemna-verification.Rmd | 11 +- 17 files changed, 846 insertions(+), 152 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 520c5e3..c2b0065 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,16 +1,17 @@ Package: lemna Title: Lemna Ecotox Effect Model -Version: 0.9.2 +Version: 1.0.0 Authors@R: c( person("Nils", "Kehrein", , "nils.kehrein@gmail.com", c("aut", "cre")), person("SETAC Europe IG Effect Modeling", role="ccp") ) -Description: An implementation of model equations and default parameters for the - toxicokinetic-toxicodynamic (TKTD) model of the Lemna (duckweed) aquatic - plant. Lemna is a standard test macrophyte used in ecotox effect studies. - The model was described and published by the SETAC Europe Interest Group - Effect Modeling. It is a refined description of the Lemna TKTD model - published by Schmitt et al. (2013) . +Description: The reference implementation of model equations and default + parameters for the toxicokinetic-toxicodynamic (TKTD) model of the Lemna + (duckweed) aquatic plant. Lemna is a standard test macrophyte used in ecotox + effect studies. The model was described and published by the SETAC Europe + Interest Group Effect Modeling. It is a refined description of the Lemna + TKTD model published by Schmitt et al. (2013) + . URL: https://github.com/nkehrein/lemna BugReports: https://github.com/nkehrein/lemna/issues License: MIT + file LICENSE diff --git a/LICENSE b/LICENSE index bd1b4c9..bbbcd72 100644 --- a/LICENSE +++ b/LICENSE @@ -1,2 +1,2 @@ -YEAR: 2021 +YEAR: 2022 COPYRIGHT HOLDER: Nils Kehrein diff --git a/LICENSE.md b/LICENSE.md index f26f2f5..5713be0 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -1,6 +1,6 @@ # MIT License -Copyright (c) 2021 Nils Kehrein +Copyright (c) 2022 Nils Kehrein Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/NEWS.md b/NEWS.md index 67405ca..e0e521a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,14 @@ +# lemna 1.0.0 + +* Documentation adapted to reference report version 1.1 +* Added a warning message in case removed parameter `BM_threshold` is used + +# lemna 0.9.2 + +* Minor adaption of the biomass ODE according to the draft report version 1.1, + handling of low biomass densities was simplified to use only one parameter, + `BM_min` + # lemna 0.9.1 * Documentation improved and typos fixed diff --git a/R/lemna.R b/R/lemna.R index c00543f..41615b5 100644 --- a/R/lemna.R +++ b/R/lemna.R @@ -121,11 +121,12 @@ #' @export #' #' @references -#' Klein J., Cedergreen N., Heine S., Reichenberger S., Rendal C., -#' Schmitt W., Hommen U., 2021: Refined description of the *Lemna* TKTD growth model -#' based on *Schmitt et al.* (2013) – equation system and default parameters. +#' 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, uploaded on 22. Sept. 2021. +#' Modeling. Version 1.1, uploaded on 09 May 2022. #' #' #' @examples diff --git a/R/param.R b/R/param.R index c0979ce..4f23d10 100644 --- a/R/param.R +++ b/R/param.R @@ -1,6 +1,6 @@ #' Default parameters #' -#' Returns the default Lemna model parameters as reported by Klein et al. (2021). +#' Returns the default Lemna model parameters as reported by Klein *et al.* (2022). #' #' ## Model parameters #' @@ -50,11 +50,12 @@ #' @export #' #' @references -#' Klein J., Cedergreen N., Heine S., Reichenberger S., Rendal C., -#' Schmitt W., Hommen U., 2021: Refined description of the *Lemna* TKTD growth model -#' based on *Schmitt et al.* (2013) – equation system and default parameters. +#' 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, uploaded on 22. Sept. 2021. +#' Modeling. Version 1.1, uploaded on 09 May 2022. #' #' #' @examples diff --git a/README.Rmd b/README.Rmd index 4fa5829..0b0283e 100644 --- a/README.Rmd +++ b/README.Rmd @@ -26,8 +26,9 @@ It implements model equations and default parameters to simulate the toxicokinetic-toxicodynamic (TKTD) model of the *Lemna* aquatic plant. *Lemna* is a standard test macrophyte used in ecotox effect studies. The model was described and published by the *SETAC Europe Interest Group Effect Modeling* -(Klein *et al.* 2021). It is a refined description of the *Lemna* TKTD model -published by Schmitt *et al.* (2013). +(Klein *et al.* 2022). It is a refined description of the *Lemna* TKTD model +published by Schmitt *et al.* (2013). This package contains the model's reference +implementation which is provided by the *SETAC* interest group. ## Installation @@ -81,11 +82,12 @@ If you find any issues or bugs within the package, please create a ## References -- Klein J., Cedergreen N., Heine S., Reichenberger S., Rendal C., - Schmitt W., Hommen U., 2021: Refined description of the *Lemna* TKTD growth model - based on *Schmitt et al.* (2013) – equation system and default parameters. +- 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, uploaded on 22. Sept. 2021. + Modeling. Version 1.1, uploaded on 09 May 2022. https://www.setac.org/group/SEIGEffectModeling - Schmitt W., Bruns E., Dollinger M., Sowig P., 2013: Mechanistic TK/TD-model simulating the effect of growth inhibitors on *Lemna* populations. Ecol Model diff --git a/README.md b/README.md index 942429e..819cff4 100644 --- a/README.md +++ b/README.md @@ -14,8 +14,10 @@ default parameters to simulate the toxicokinetic-toxicodynamic (TKTD) model of the *Lemna* aquatic plant. *Lemna* is a standard test macrophyte used in ecotox effect studies. The model was described and published by the *SETAC Europe Interest Group Effect Modeling* (Klein -*et al.* 2021). It is a refined description of the *Lemna* TKTD model -published by Schmitt *et al.* (2013). +*et al.* 2022). It is a refined description of the *Lemna* TKTD model +published by Schmitt *et al.* (2013). This package contains the model’s +reference implementation which is provided by the *SETAC* interest +group. ## Installation @@ -80,12 +82,13 @@ issue](https://github.com/nkehrein/lemna/issues) on GitHub. ## References -- Klein J., Cedergreen N., Heine S., Reichenberger S., Rendal C., - Schmitt W., Hommen U., 2021: Refined description of the *Lemna* TKTD - growth model based on *Schmitt et al.* (2013) – equation system and - default parameters. Report of the working group *Lemna* of the SETAC - Europe Interest Group Effect Modeling. Version 1, uploaded on 22. - Sept. 2021. +- 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. + - Schmitt W., Bruns E., Dollinger M., Sowig P., 2013: Mechanistic TK/TD-model simulating the effect of growth inhibitors on *Lemna* populations. Ecol Model 255, pp. 1-10. DOI: diff --git a/dev/release.R b/dev/release.R index de905ef..864d4cb 100644 --- a/dev/release.R +++ b/dev/release.R @@ -14,7 +14,7 @@ devtools::document() devtools::install(upgrade="never", quick=TRUE) devtools::build_vignettes(quiet=FALSE, install=FALSE) if(!dir.exists("docs")) dir.create("docs") -file.copy(list.files("doc", pattern="*.html", full.names=TRUE), "docs") +file.copy(list.files("doc", pattern="\\.html$", full.names=TRUE), "docs", recursive=TRUE) unlink("doc", recursive = TRUE) # build package files diff --git a/docs/lemna-introduction.html b/docs/lemna-introduction.html index 3ac272c..e08b8bb 100644 --- a/docs/lemna-introduction.html +++ b/docs/lemna-introduction.html @@ -15,7 +15,19 @@ Introduction to the Lemna package - + @@ -140,7 +332,7 @@

Introduction to the Lemna package

Nils Kehrein

-

27 January, 2022

+

10 May, 2022

-

The lemna package provides model equations and some useful helpers to simulate the growth of Lemna (duckweed) aquatic plant populations. Lemna is a standard test macrophyte used in ecotox effect studies. The model was described and published by the SETAC Europe Interest Group Effect Modeling (Klein et al. 2021).

-

The model’s main state variable is biomass, or BM for short, of the simulated Lemna population. Growth of Lemna is influenced by environmental variables such as temperature, irradiation, nutrient concentrations, population density, and toxicant concentration in the surrounding medium. To consider the influence of toxicants on the plants, a one-compartment model was assumed by the authors for the mass-balance of internal toxicant mass. The total amount of internal toxicant mass is represented by state-variable M_int. The combination of state variables BM and M_int fully describe the state of the model system at any point in time.

-

To simulate a Lemna population, one has to define a scenario that consists of the following data:

+

The lemna package provides model equations and some useful +helpers to simulate the growth of Lemna (duckweed) aquatic +plant populations. Lemna is a standard test macrophyte used in +ecotox effect studies. The model was described and published by the +SETAC Europe Interest Group Effect Modeling (Klein et al. +2022).

+

The model’s main state variable is biomass, or BM for +short, of the simulated Lemna population. Growth of +Lemna is influenced by environmental variables such as +temperature, irradiation, nutrient concentrations, population density, +and toxicant concentration in the surrounding medium. To consider the +influence of toxicants on the plants, a one-compartment model was +assumed by the authors for the mass-balance of internal toxicant mass. +The total amount of internal toxicant mass is represented by +state-variable M_int. The combination of state variables +BM and M_int fully describe the state of the +model system at any point in time.

+

To simulate a Lemna population, one has to define a +scenario that consists of the following data:

  • Initial system state
  • A time period to simulate
  • Model parameters
  • Environmental variables
-

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.

+

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:

+

To make functions and sample datasets of the lemna package +available in your R workspace, load the library first:

library(lemna)
-

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:

+

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
@@ -185,7 +405,11 @@ 
Core functions

 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:

+

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
@@ -203,23 +427,40 @@ 
Core functions

 #> 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:

+

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:

plot(lemna(metsulfuron))

-

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:

+

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:

+

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:

?metsulfuron

Tutorial

Simulate the Lemna growth model

-

To simulate a Lemna population, one has to pass the four mandatory scenario elements to the lemna() function:

+

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
@@ -254,7 +495,11 @@ 
Simulate the Lemna growth model

 #> 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:

+

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
@@ -269,7 +514,17 @@ 
Simulate the Lemna growth model

 #> 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:

+

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
@@ -281,7 +536,11 @@ 
Simulate the Lemna growth model

 #> 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:

+

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
@@ -295,16 +554,40 @@ 
Simulate the Lemna growth model

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).

+

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:

+

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))
@@ -315,7 +598,12 @@ 
Using environmental time-series

 # 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.

+

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
@@ -335,7 +623,11 @@ 
Using environmental time-series

 # 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():

+

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"
@@ -343,7 +635,9 @@ 
Using environmental time-series

 # 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:

+

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
@@ -351,7 +645,13 @@ 
Using environmental time-series

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:

+

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
@@ -361,10 +661,13 @@ 
Using simulation results

 #> 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:

+

To get an initial impression of a scenario and its results, simply +pass the simulation result to the plot() function:

plot(myresult)

-

As an example, we will analyze if and how the internal toxicant concentration (C_int) correlates with the internal toxicant mass (M_int):

+

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:
@@ -384,29 +687,57 @@ 
Using simulation results

 #> 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.

+

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:

+

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:

+

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%.

+

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:

+

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:

+

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)
@@ -423,7 +754,8 @@ 
Create scenario objects

 #> 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():

+

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,
@@ -449,7 +781,13 @@ 
Create scenario objects

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.

+

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
@@ -459,19 +797,22 @@ 
Speed up simulations with compiled code

 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:

+

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
-#>  lemna(metsulfuron, ode_mode = "r") 12.3105 12.52985 14.294595 12.73895
-#>  lemna(metsulfuron, ode_mode = "c")  3.4660  3.62260  4.411199  3.73155
-#>        uq     max neval
-#>  15.51680 22.4626   100
-#>   4.92065 10.9884   100
+#>                                expr    min      lq     mean median      uq
+#>  lemna(metsulfuron, ode_mode = "r") 6.2229 6.51900 7.425877 6.7554 7.84990
+#>  lemna(metsulfuron, ode_mode = "c") 1.7881 1.89195 2.258368 1.9977 2.24125
+#>      max neval
+#>  11.6633   100
+#>   6.7320   100
 
 # Benchmark the more complex and longer focusd1 scenario
 microbenchmark::microbenchmark(
@@ -480,13 +821,18 @@ 
Speed up simulations with compiled code

   times = 10
 )
 #> Unit: milliseconds
-#>                            expr       min        lq      mean    median
-#>  lemna(focusd1, ode_mode = "r") 3608.1600 3684.6830 3811.9995 3800.6699
-#>  lemna(focusd1, ode_mode = "c")   71.4479   73.1658   78.4811   75.7557
+#>                            expr       min        lq       mean    median
+#>  lemna(focusd1, ode_mode = "r") 1668.6090 1686.6235 1819.10066 1820.5735
+#>  lemna(focusd1, ode_mode = "c")   35.8434   37.8879   39.79211   38.9555
 #>         uq       max neval
-#>  3879.3347 4137.7390    10
-#>    85.5938   92.8994    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:

+#> 1923.1008 1980.6964 10 +#> 39.9389 49.0085 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
@@ -542,7 +888,12 @@ 
Speed up simulations with compiled code

 

 
References

 
    
-
  • Klein J., Cedergreen N., Heine S., Reichenberger S., Rendal C., Schmitt W., Hommen U., 2021: Refined description of the Lemna TKTD growth model based on Schmitt et al. (2013) – equation system and default parameters. Report of the working group Lemna of the SETAC Europe Interest Group Effect Modeling. Version 1, uploaded on 22. Sept. 2021. https://www.setac.org/group/SEIGEffectModeling
    • 
+
  • 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
    • 
 

 

 
diff --git a/docs/lemna-verification.html b/docs/lemna-verification.html
index b445f29..13ec1ad 100644
--- a/docs/lemna-verification.html
+++ b/docs/lemna-verification.html
@@ -15,7 +15,19 @@
 
 Model verification
 
-
+
 
 
 
 
 
@@ -140,7 +332,7 @@
 
 
Model verification

 
Nils Kehrein

-
27 January, 2022

+
10 May, 2022

 
 
 

@@ -149,11 +341,15 @@ 
27 January, 2022

 
  • Schmitt et al. model
    • 
 
  • Unlimited growth scenarios
 
    • 
-
  • Environmental variability scenarios
+
  • Environmental
+variability scenarios
 
      
 
    • FOCUS D1 Ditch
      • 
 
    • FOCUS D2 Ditch
      • 
@@ -167,13 +363,29 @@ 
      27 January, 2022
      
 
 
      
 
      Introduction
      
-
      This document aims at verifying the implementation of the Lemna toxicokinetic-toxicodynamic (TKTD) model as described by Klein et al. (2021). The model is a refined description of the Lemna TKTD model based on Schmitt et al. (2013).
      
-
      First, this model’s output is compared to published results which were created with the reference implementation of the Schmitt et al. model. For unlimited growth scenarios, both models are expected to yield identical results for identical inputs. Second, model dynamics are compared for additional scenarios which simulate changing environmental conditions, such as under field study conditions. In this case, the models are not expected to yield identical but equivalent behavior due to differences between the models’ response functions.
      
+
      This document aims at verifying the implementation of the
+Lemna toxicokinetic-toxicodynamic (TKTD) model as described by
+Klein et al. (2022). The model is a refined description of the
+Lemna TKTD model based on Schmitt et al. (2013).
      
+
      First, this model’s output is compared to published results which
+were created with the reference implementation of the Schmitt et
+al. model. For unlimited growth scenarios, both models are expected
+to yield identical results for identical inputs. Second, model dynamics
+are compared for additional scenarios which simulate changing
+environmental conditions, such as under field study conditions. In this
+case, the models are not expected to yield identical but equivalent
+behavior due to differences between the models’ response functions.
      
 
      
 
      
 
      Schmitt et al. model
      
-
      The outputs of the Lemna TKTD model by Schmitt et al. are used to verify the implementation of the Klein et al. model. The reference implementation of Schmitt et al. is included in this package’s GitHub repository but is not an official part of the package.
      
-
      Sourcing files mmc2.r and mmc3.r gives access to the Schmitt et al. model equations as well as to the fitted parameter set for the substance metsulfuron-methyl.
      
+
      The outputs of the Lemna TKTD model by Schmitt et
+al. are used to verify the implementation of the Klein et
+al. model. The reference implementation of Schmitt et al.
+is included in this package’s GitHub repository but is not an official
+part of the package.
      
+
      Sourcing files mmc2.r and mmc3.r gives
+access to the Schmitt et al. model equations as well as to the
+fitted parameter set for the substance metsulfuron-methyl.
      
 
      library(lemna)   # lemna model
 # load packages to ease plotting and data wrangling
 library(dplyr)   # for data preparation
@@ -200,10 +412,17 @@ 
      Schmitt et al. model
      
 
      
 
      
 
      Unlimited growth scenarios
      
-
      This section simulates several Lemna growth scenarios under unlimited growth conditions using the Klein et al. model. Model dynamics and numerical results are compared to results of the Schmitt et al. model.
      
+
      This section simulates several Lemna growth scenarios under
+unlimited growth conditions using the Klein et al. model. Model
+dynamics and numerical results are compared to results of the
+Schmitt et al. model.
      
 
      
 
      Seven days exposure, seven days recovery
      
-
      Figure 2 in Schmitt et al. (2013) presents several scenarios which simulate the growth of Lemna exposed to several concentrations of metsulfuron-methyl. Exposure lasts for seven days, followed by a seven day recovery period. Symbols show observed data and lines as calculated with fitted model parameters.
      
+
      Figure 2 in Schmitt et al. (2013) presents several scenarios
+which simulate the growth of Lemna exposed to several
+concentrations of metsulfuron-methyl. Exposure lasts for seven
+days, followed by a seven day recovery period. Symbols show observed
+data and lines as calculated with fitted model parameters.
      
 
      # simulate the 7d exposure, 7d recovery experiment
 simulate_7d <- function(model=c("schmitt","klein"), ...) {
   model <- match.arg(model)
@@ -261,7 +480,9 @@ 
      Seven days exposure, seven days recovery
      
   labs(x="Time (days)", y="Number of Fronds", color="Conc. (ug/L)",
        title="Schmitt et al. model: 7d exposure, 7d recovery")
      
 
      
-
      The depicted model dynamics of the Schmitt et al. are identical to published results. Simulating the same scenarios using the Klein et al. model yields identical dynamics:
      
+
      The depicted model dynamics of the Schmitt et al. are
+identical to published results. Simulating the same scenarios using the
+Klein et al. model yields identical dynamics:
      
 
      # simulate using the Klein et al. model
 df_k <- simulate_7d(model="klein")
 
@@ -282,7 +503,9 @@ 
      Seven days exposure, seven days recovery
      
   labs(x="Time (days)", y="Number of Fronds", color="Conc. (ug/L)",
        title="Klein et al. model: 7d exposure, 7d recovery")
      
 
      
-
      A comparison of numerical values between the two models shows that simulated biomass deviates at most by 0.15%. A plot of relative errors describes small but consistent deviations in biomass results:
      
+
      A comparison of numerical values between the two models shows that
+simulated biomass deviates at most by 0.15%. A plot of relative errors
+describes small but consistent deviations in biomass results:
      
 
      df_s %>%
   mutate(error = 100 * (df_k$BM/BM - 1)) %>%
   ggplot(aes(time, error, color=level)) +
@@ -291,7 +514,14 @@ 
      Seven days exposure, seven days recovery
      
   labs(x="Time (days)", y="Relative error (%)", color="Conc. (ug/L)",
        title="Deviation in biomass in Klein model relative to Schmitt")
      
 
      
-
      The deviations are a result of imprecisions introduced by the numerical ODE solver. Numerical precision can be increased by, for example, decreasing the solver’s step length in time. This can be achieved by decreasing the solver parameter hmax from the model’s default value of hmax=0.1 to e.g. hmax=0.01. The value of 0.01 is equivalent to a step length of about 15 minutes. This decreases the observed deviations significantly to a maximum of about 0.007%:
      
+
      The deviations are a result of imprecisions introduced by the
+numerical ODE solver. Numerical precision can be increased by, for
+example, decreasing the solver’s step length in time. This can be
+achieved by decreasing the solver parameter hmax from the
+model’s default value of hmax=0.1 to
+e.g. hmax=0.01. The value of 0.01 is
+equivalent to a step length of about 15 minutes. This decreases the
+observed deviations significantly to a maximum of about 0.007%:
      
 
      df_s <-  simulate_7d(model="schmitt", hmax=0.01)
 df_k <-  simulate_7d(model="klein", hmax=0.01)
 
@@ -303,11 +533,17 @@ 
      Seven days exposure, seven days recovery
      
   labs(x="Time (days)", y="Relative error (%)", color="Conc. (ug/L)",
        title="Deviation in biomass in Klein model relative to Schmitt, hmax=0.01")
      
 
      
-
      It can be concluded that model dynamics and numerical results are identical as far as numerical precision allows.
      
+
      It can be concluded that model dynamics and numerical results are
+identical as far as numerical precision allows.
      
 
      
 
      
 
      Two days of exposure, twelve days of recovery
      
-
      The right-hand side of Figure 3 in Hommen et al. (2015) presents additional scenarios which simulate the growth of Lemna exposed to several concentrations of metsulfuron-methyl. Exposure lasts for two days, followed by a twelve day recovery period. Symbols show observed data and lines as calculated with fitted model parameters:
      
+
      The right-hand side of Figure 3 in Hommen et al. (2015)
+presents additional scenarios which simulate the growth of
+Lemna exposed to several concentrations of
+metsulfuron-methyl. Exposure lasts for two days, followed by a
+twelve day recovery period. Symbols show observed data and lines as
+calculated with fitted model parameters:
      
 
      # simulate the 2d exposure, 12d recovery experiment
 simulate_2d <- function(model=c("schmitt","klein"), ...) {
   model <- match.arg(model)
@@ -367,7 +603,10 @@ 
      Two days of exposure, twelve days of recovery
      
   labs(x="Time (days)", y="Number of Fronds", color="Conc. (ug/L)",
        title="Klein et al. model: 2d exposure, 12d recovery")
      
 
      
-
      Observed model dynamics of the Klein et al. model are identical to published results and fit equally well with the observed frond numbers reported by Hommen et al. A comparison of numerical values of simulated biomass yields:
      
+
      Observed model dynamics of the Klein et al. model are
+identical to published results and fit equally well with the observed
+frond numbers reported by Hommen et al. A comparison of
+numerical values of simulated biomass yields:
      
 
      # simulate using the Schmitt et al. model
 df_s <-  simulate_2d(model="schmitt", hmax=0.01)
 # simulate using the Klein et al. model
@@ -380,11 +619,17 @@ 
      Two days of exposure, twelve days of recovery
      
   abs() %>%
   max() 
 #> [1] 0.004476714
      
-
      Numerical biomass values deviate at most by about 0.004% between the two models. It can be concluded that model dynamics and numerical results of both models are identical for the presented scenarios.
      
+
      Numerical biomass values deviate at most by about 0.004% between the
+two models. It can be concluded that model dynamics and numerical
+results of both models are identical for the presented scenarios.
      
 
      
 
      
 
      Complex exposure patterns
      
-
      Figure 4 in Schmitt et al. (2013) depicts the growth of Lemna exposed to three exposure patterns of several different concentrations. The patterns include constant exposure, repeating intervals of four days out of seven of exposure, as well as repeating intervals of two days out of seven of exposure.
      
+
      Figure 4 in Schmitt et al. (2013) depicts the growth of
+Lemna exposed to three exposure patterns of several different
+concentrations. The patterns include constant exposure, repeating
+intervals of four days out of seven of exposure, as well as repeating
+intervals of two days out of seven of exposure.
      
 
      First, the constant exposure pattern (Figure 4a) is simulated:
      
 
      # simulate a specific exposure pattern
 simulate_pattern <- function(model=c("schmitt","klein"), levels, expo_fun, ...) {
@@ -446,7 +691,8 @@ 
      Complex exposure patterns
      
   labs(x="Time (days)", y="Frond area (cm²)", color="Conc. (ug/L)",
        title="Klein et al. model: constant exposure")
      
 
      
-
      Next, the four out of seven days exposure pattern (Figure 4b) is simulated:
      
+
      Next, the four out of seven days exposure pattern (Figure 4b) is
+simulated:
      
 
      # exposure pattern: repeats 4d exposure, 3d recovery for 42 days,
 # followed by recovery for 8 days
 expo_4of7 <- function(conc=1) {
@@ -476,7 +722,8 @@ 
      Complex exposure patterns
      
   labs(x="Time (days)", y="Frond area (cm²)", color="Conc. (ug/L)",
        title="Klein et al. model: 4 out of 7d of exposure")
      
 
      
-
      Finally, the two out of seven days exposure pattern (Figure 4c) is simulated:
      
+
      Finally, the two out of seven days exposure pattern (Figure 4c) is
+simulated:
      
 
      # exposure pattern: repeats 2d exposure, 5d recovery for 42 days,
 # followed by recovery for 8 days
 expo_2of7 <- function(conc=1) {
@@ -506,7 +753,9 @@ 
      Complex exposure patterns
      
   labs(x="Time (days)", y="Frond area (cm²)", color="Conc. (ug/L)",
        title="Klein et al. model: 2 out of 7d of exposure")
      
 
      
-
      Comparing numerical results of the Klein et al. model with the Schmitt et al. implementation yields the following relative errors in percent:
      
+
      Comparing numerical results of the Klein et al. model with
+the Schmitt et al. implementation yields the following relative
+errors in percent:
      
 
      # relative error (%) for the constant exposure pattern
 simulate_pattern(model="schmitt", levels=levels_const, expo_fun=expo_const) %>%
   mutate(error = 100 * (df_p1$BM/BM - 1)) %>%
@@ -530,16 +779,31 @@ 
      Complex exposure patterns
      
   abs() %>%
   max() 
 #> [1] 0.02197463
      
-
      It can be observed that the Klein et al. model reproduces the same model dynamics and biomass amounts as reported by Schmitt et al. under unlimited growth conditions.
      
+
      It can be observed that the Klein et al. model reproduces
+the same model dynamics and biomass amounts as reported by Schmitt
+et al. under unlimited growth conditions.
      
 
      
 
      
 
      
 
      Environmental variability scenarios
      
-
      In this section, simulated biomass dynamics of the Klein et al. model are compared to published results under conditions of changing environmental variables. In contrast to the Schmitt et al. model, it applies Liebig’s Law so that growth is controlled by the most limiting resource, i.e. the limiting factor. The environmental variables temperature, global radiation, phosphorus, and nitrate are considered for Liebig’s law (Klein et al. 2015).
      
-
      The following graphs will replicate Figure 5 from Hommen et al. (2015). The simulations make use of exposure time series and time series of environmental variables that were extracted from the three FOCUS Surface Water exposure scenarios D1 Ditch, D2 Ditch, and R3 Stream. Substance specific parameters follow the fitted parameters reported by Schmitt et al. (2013).
      
+
      In this section, simulated biomass dynamics of the Klein et
+al. model are compared to published results under conditions of
+changing environmental variables. In contrast to the Schmitt et
+al. model, it applies Liebig’s Law so that growth is controlled by
+the most limiting resource, i.e. the limiting factor. The environmental
+variables temperature, global radiation, phosphorus, and nitrate are
+considered for Liebig’s law (Klein et al. 2015).
      
+
      The following graphs will replicate Figure 5 from Hommen et
+al. (2015). The simulations make use of exposure time series and
+time series of environmental variables that were extracted from the
+three FOCUS Surface Water exposure scenarios D1 Ditch, D2
+Ditch, and R3 Stream. Substance specific parameters follow
+the fitted parameters reported by Schmitt et al. (2013).
      
 
      
 
      FOCUS D1 Ditch
      
-
      The following code simulates Lemna population growth for the example FOCUS exposure pattern D1 Ditch multiplied by factors from 1 to 100:
      
+
      The following code simulates Lemna population growth for the
+example FOCUS exposure pattern D1 Ditch multiplied by factors from 1 to
+100:
      
 
      library(tidyr)  # for additional data preparation
 library(furrr)  # for multi-threading capabilities
 
@@ -592,7 +856,10 @@ 
      FOCUS D1 Ditch
      
 df_k <- future_map_dfr(factors, simulate_focus, model="klein", scenario=focusd1)
 rainbow_plot(df_k, focusd1$envir$conc, 2000, "D1 Ditch: Klein et al. model, L. gibba")
      
 
      
-
      It can be observed that the Klein et al. model generates biomass dynamics that are very similar to the Schmitt et al. results. Comparing biomass values for each time point and multiplication factor yields the following graph:
      
+
      It can be observed that the Klein et al. model generates
+biomass dynamics that are very similar to the Schmitt et al.
+results. Comparing biomass values for each time point and multiplication
+factor yields the following graph:
      
 
      bm_dev_plot <- function(res_klein, res_schmitt, title) {
   fs <- unique(res_klein$factor)
   
@@ -609,8 +876,22 @@ 
      FOCUS D1 Ditch
      
 
 bm_dev_plot(df_k, df_s, "D1 Ditch")
      
 
      
-
      For D1 Ditch, the Klein et al. model shows a trend to a slightly larger biomass over the course of the simulated year. Achieved biomass levels are nevertheless very similar in both models. The perceived differences seem to originate from the fact that model results are slightly shifted in time. This observation can be explained by the photosynthesis response function of the Klein et al. model, where growth is controlled by the most limiting factor. Unlike the Schmitt et al. model, where all response functions are multiplied. This allows the population to react quicker to changing environmental conditions and achieve slightly larger biomass in the Klein et al. model. Generally, both models generate similar growth dynamics.
      
-
      The following table reproduces Table 2 from Hommen et al. (2015). It lists the number of days with deviations from controls exceeding different magnitudes (% dev) depending on the exposure multiplication factor in the D1 Ditch scenario:
      
+
      For D1 Ditch, the Klein et al. model shows a trend to a
+slightly larger biomass over the course of the simulated year. Achieved
+biomass levels are nevertheless very similar in both models. The
+perceived differences seem to originate from the fact that model results
+are slightly shifted in time. This observation can be explained by the
+photosynthesis response function of the Klein et al. model,
+where growth is controlled by the most limiting factor. Unlike the
+Schmitt et al. model, where all response functions are
+multiplied. This allows the population to react quicker to changing
+environmental conditions and achieve slightly larger biomass in the
+Klein et al. model. Generally, both models generate similar
+growth dynamics.
      
+
      The following table reproduces Table 2 from Hommen et al.
+(2015). It lists the number of days with deviations from controls
+exceeding different magnitudes (% dev) depending on the exposure
+multiplication factor in the D1 Ditch scenario:
      
 
      df_k %>%
   pivot_wider(id_cols=time, names_from=factor, values_from=BM) -> df.wide
 
@@ -787,11 +1068,18 @@ 
      FOCUS D1 Ditch
      
 
 
 
-
      Predictions of the Klein et al. model are again very similar to the results reported by Hommen et al. (2015). The number of days with deviations do appear to have a small trend to lower numbers compared to the results of the Schmitt et al. model. As an example: for a multiplication factor of 4, 47 days experience a deviation in biomass of at least 1% compared to 50 days as reported by Hommen et al.
      
+
      Predictions of the Klein et al. model are again very similar
+to the results reported by Hommen et al. (2015). The number of
+days with deviations do appear to have a small trend to lower numbers
+compared to the results of the Schmitt et al. model. As an
+example: for a multiplication factor of 4, 47 days experience a
+deviation in biomass of at least 1% compared to 50 days as reported by
+Hommen et al.
      
 
      
 
      
 
      FOCUS D2 Ditch
      
-
      Simulating Lemna population growth for the example FOCUS exposure pattern D2 Ditch multiplied by factors from 1 to 100:
      
+
      Simulating Lemna population growth for the example FOCUS
+exposure pattern D2 Ditch multiplied by factors from 1 to 100:
      
 
      # simulate using the Schmitt et al. model
 df_s <- future_map_dfr(factors, simulate_focus, model="schmitt", scenario=focusd2)
 rainbow_plot(df_s, focusd2$envir$conc, 100, "D2 Ditch: Schmitt et al. model, L. gibba")
      
@@ -802,11 +1090,18 @@ 
      FOCUS D2 Ditch
      
 
      
 
      bm_dev_plot(df_k, df_s, "D2 Ditch")
      
 
      
-
      Biomass dynamics and absolute numbers are again very similar to the behavior reported by Hommen et al. (2015). Generally, the Klein et al. model shows a trend towards temporary larger biomass values which can be explained by the model reacting quicker to changes in environmental variables. Biomass levels at the end of the simulated year are slightly lower for low and moderate exposure levels than predictions of the Schmitt et al. model.
      
+
      Biomass dynamics and absolute numbers are again very similar to the
+behavior reported by Hommen et al. (2015). Generally, the
+Klein et al. model shows a trend towards temporary larger
+biomass values which can be explained by the model reacting quicker to
+changes in environmental variables. Biomass levels at the end of the
+simulated year are slightly lower for low and moderate exposure levels
+than predictions of the Schmitt et al. model.
      
 
      
 
      
 
      FOCUS R3 Stream
      
-
      Simulating Lemna population growth for the example FOCUS exposure pattern R3 Stream multiplied by factors from 1 to 1000:
      
+
      Simulating Lemna population growth for the example FOCUS
+exposure pattern R3 Stream multiplied by factors from 1 to 1000:
      
 
      # multiplication factors for the exposure series
 factors <- c(0, round(10^seq(0,3,0.3), 1))
 
@@ -822,24 +1117,45 @@ 
      FOCUS R3 Stream
      
 
      
 
      # disable multithreading
 future::plan("sequential")
      
-
      Again, biomass dynamics and absolute numbers are very similar to the behavior reported by Hommen et al. (2015) with a small trend towards a larger biomass in the Klein et al. model.
      
+
      Again, biomass dynamics and absolute numbers are very similar to the
+behavior reported by Hommen et al. (2015) with a small trend
+towards a larger biomass in the Klein et al. model.
      
 
      
 
      
 
      
 
      Compiled code
      
-
      Model equations implemented as compiled code in C are not part of this verification document for reasons of brevity. However, results of the C code are identical as far as numerical precision allows. This is ensured by automated unit tests which are shipped as part of the lemna source package. The suite of verification tests can be run manually by:
      
+
      Model equations implemented as compiled code in C are not
+part of this verification document for reasons of brevity. However,
+results of the C code are identical as far as numerical
+precision allows. This is ensured by automated unit tests which are
+shipped as part of the lemna source package. The suite of
+verification tests can be run manually by:
      
 
      devtools::test(pkg="lemna", filter="verify")
      
 
      
 
      
 
      Acknowledgements
      
-
      The author would like to thank U. Hommen and S. Heine for providing the original data sets and scripts which were used for past publications, as well as J. Klein and J. Witt for their feedback and suggestions.
      
+
      The author would like to thank U. Hommen and S. Heine for providing
+the original data sets and scripts which were used for past
+publications, as well as J. Klein and J. Witt for their feedback and
+suggestions.
      
 
      
 
      
 
      References
      
 
        
-
      • Hommen U., Schmitt W., Heine S., Brock Theo C.M., Duquesne S., Manson P., Meregalli G., Ochoa-Acuña H., van Vliet P., Arts G., 2015: How TK-TD and Population Models for Aquatic Macrophytes Could Support the Risk Assessment for Plant Protection Products. Integr Environ Assess Manag 12(1), pp. 82-95. DOI: 10.1002/ieam.1715
        • 
-
      • Klein J., Cedergreen N., Heine S., Reichenberger S., Rendal C., Schmitt W., Hommen U., 2021: Refined description of the Lemna TKTD growth model based on Schmitt et al. (2013) – equation system and default parameters. Report of the working group Lemna of the SETAC Europe Interest Group Effect Modeling. Version 1, uploaded on 22. Sept. 2021. https://www.setac.org/group/SEIGEffectModeling
        • 
-
      • Schmitt W., Bruns E., Dollinger M., Sowig P., 2013: Mechanistic TK/TD-model simulating the effect of growth inhibitors on Lemna populations. Ecol Model 255, pp. 1-10. DOI: 10.1016/j.ecolmodel.2013.01.017
        • 
+
      • Hommen U., Schmitt W., Heine S., Brock Theo C.M., Duquesne S.,
+Manson P., Meregalli G., Ochoa-Acuña H., van Vliet P., Arts G., 2015:
+How TK-TD and Population Models for Aquatic Macrophytes Could Support
+the Risk Assessment for Plant Protection Products. Integr Environ Assess
+Manag 12(1), pp. 82-95. DOI: 10.1002/ieam.1715
        • 
+
      • 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
        • 
+
      • Schmitt W., Bruns E., Dollinger M., Sowig P., 2013: Mechanistic
+TK/TD-model simulating the effect of growth inhibitors on Lemna
+populations. Ecol Model 255, pp. 1-10. DOI: 10.1016/j.ecolmodel.2013.01.017
        • 
 
      
 
      
 
diff --git a/man/lemna.Rd b/man/lemna.Rd
index adb25c0..ae0590b 100644
--- a/man/lemna.Rd
+++ b/man/lemna.Rd
@@ -250,10 +250,11 @@ scenario <- new_lemna_scenario(
 lemna(scenario)
 }
 \references{
-Klein J., Cedergreen N., Heine S., Reichenberger S., Rendal C.,
-Schmitt W., Hommen U., 2021: Refined description of the \emph{Lemna} TKTD growth model
-based on \emph{Schmitt et al.} (2013) – equation system and default parameters.
+Klein J., Cedergreen N., Heine S., Kehrein N., Reichenberger S., Rendal C.,
+Schmitt W., Hommen U., 2022: Refined description of the \emph{Lemna} TKTD growth model
+based on \emph{Schmitt et al.} (2013) – equation system and default parameters,
+implementation in R.
 Report of the working group \emph{Lemna} of the SETAC Europe Interest Group Effect
-Modeling. Version 1, uploaded on 22. Sept. 2021.
+Modeling. Version 1.1, uploaded on 09 May 2022.
 \url{https://www.setac.org/group/SEIGEffectModeling}
 }
diff --git a/man/param_defaults.Rd b/man/param_defaults.Rd
index 60de7d7..49939e5 100644
--- a/man/param_defaults.Rd
+++ b/man/param_defaults.Rd
@@ -17,7 +17,7 @@ defaults}
 named \code{list}
 }
 \description{
-Returns the default Lemna model parameters as reported by Klein et al. (2021).
+Returns the default Lemna model parameters as reported by Klein \emph{et al.} (2022).
 }
 \details{
 \subsection{Model parameters}{
@@ -110,10 +110,11 @@ param_defaults(list(
 param_new()
 }
 \references{
-Klein J., Cedergreen N., Heine S., Reichenberger S., Rendal C.,
-Schmitt W., Hommen U., 2021: Refined description of the \emph{Lemna} TKTD growth model
-based on \emph{Schmitt et al.} (2013) – equation system and default parameters.
+Klein J., Cedergreen N., Heine S., Kehrein N., Reichenberger S., Rendal C.,
+Schmitt W., Hommen U., 2022: Refined description of the \emph{Lemna} TKTD growth model
+based on \emph{Schmitt et al.} (2013) – equation system and default parameters,
+implementation in R.
 Report of the working group \emph{Lemna} of the SETAC Europe Interest Group Effect
-Modeling. Version 1, uploaded on 22. Sept. 2021.
+Modeling. Version 1.1, uploaded on 09 May 2022.
 \url{https://www.setac.org/group/SEIGEffectModeling}
 }
diff --git a/src/lemna.c b/src/lemna.c
index be87d11..4d8a5c8 100644
--- a/src/lemna.c
+++ b/src/lemna.c
@@ -1,6 +1,6 @@
 /*******************************************************************
  *
- * Lemna model as described by Klein et al. (2021)
+ * Lemna model as described by Klein et al. (2022)
  * SETAC Europe Interest Group Effect Modeling. Version 1.1
  *
  * The model provides additional output on intermediary variables on
diff --git a/tests/testthat/test-verify-variable.R b/tests/testthat/test-verify-variable.R
index 8a8e84a..2d5f5e7 100644
--- a/tests/testthat/test-verify-variable.R
+++ b/tests/testthat/test-verify-variable.R
@@ -1,8 +1,6 @@
 #
 # Helper functions and setup
 #
-future::plan(future::multisession())
-
 # load Lemna model by Schmitt et al. (2013)
 source("mmc2.r", local=TRUE)
 source("mmc3.r", local=TRUE)
@@ -35,6 +33,9 @@ test_that("FOCUS D1 Ditch", {
   skip_on_cran()
   skip_if_not_installed("furrr")
 
+  # enable parallel processing for this and the following tests
+  future::plan(future::multisession())
+
   factors <- c(0, round(10^seq(0,2,0.2), 1))
   df_s <- furrr::future_map_dfr(factors, simulate_focus, model="schmitt", scenario=focusd1)
   df_kr <- furrr::future_map_dfr(factors, simulate_focus, model="klein", scenario=focusd1, ode_mode="r")
@@ -70,4 +71,7 @@ test_that("FOCUS R3 Stream", {
   df_kc <- furrr::future_map_dfr(factors, simulate_focus, model="klein", scenario=focusr3, ode_mode="c")
 
   expect_equal(df_kr, df_kc, tolerance=1e-4, ignore_attr=TRUE)
+
+  # disable parallel processing
+  future::plan(future::sequential())
 })
diff --git a/vignettes/lemna-introduction.Rmd b/vignettes/lemna-introduction.Rmd
index 56a800c..5c59ef8 100644
--- a/vignettes/lemna-introduction.Rmd
+++ b/vignettes/lemna-introduction.Rmd
@@ -32,7 +32,7 @@ The *lemna* package provides model equations and some useful helpers to
 simulate the growth of *Lemna* (duckweed) aquatic plant populations. 
 *Lemna* is a standard test macrophyte used in ecotox effect studies. The
 model was described and published by the SETAC Europe Interest Group Effect
-Modeling (Klein *et al.* 2021).
+Modeling (Klein *et al.* 2022).
 
 The model's main state variable is biomass, or `BM` for short, of the simulated
 *Lemna* population. Growth of *Lemna* is influenced by environmental variables
@@ -486,9 +486,10 @@ lemna(metsulfuron, ode_mode="c", nout=18)
 
 
 ## References
-- Klein J., Cedergreen N., Heine S., Reichenberger S., Rendal C.,
-  Schmitt W., Hommen U., 2021: Refined description of the *Lemna* TKTD growth model
-  based on *Schmitt et al.* (2013) – equation system and default parameters.
+- 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, uploaded on 22. Sept. 2021.
+  Modeling. Version 1.1, uploaded on 09 May 2022.
   https://www.setac.org/group/SEIGEffectModeling
diff --git a/vignettes/lemna-verification.Rmd b/vignettes/lemna-verification.Rmd
index 5cc72a6..39f7e4d 100644
--- a/vignettes/lemna-verification.Rmd
+++ b/vignettes/lemna-verification.Rmd
@@ -30,7 +30,7 @@ knitr::opts_chunk$set(eval = !is_check)
 ## Introduction
 
 This document aims at verifying the implementation of the *Lemna* 
-toxicokinetic-toxicodynamic (TKTD) model as described by *Klein et al.* (2021).
+toxicokinetic-toxicodynamic (TKTD) model as described by *Klein et al.* (2022).
 The model is a refined description of the *Lemna* TKTD model based on
 *Schmitt et al. (2013)*.
 
@@ -701,11 +701,12 @@ as well as J. Klein and J. Witt for their feedback and suggestions.
   Aquatic Macrophytes Could Support the Risk Assessment for Plant Protection
   Products. Integr Environ Assess Manag 12(1), pp. 82-95. DOI:
   [10.1002/ieam.1715](https://doi.org/10.1002/ieam.1715)
-- Klein J., Cedergreen N., Heine S., Reichenberger S., Rendal C.,
-  Schmitt W., Hommen U., 2021: Refined description of the *Lemna* TKTD growth model
-  based on *Schmitt et al.* (2013) – equation system and default parameters.
+- 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, uploaded on 22. Sept. 2021.
+  Modeling. Version 1.1, uploaded on 09 May 2022.
   https://www.setac.org/group/SEIGEffectModeling
 - Schmitt W., Bruns E., Dollinger M., Sowig P., 2013: Mechanistic TK/TD-model
   simulating the effect of growth inhibitors on *Lemna* populations. Ecol Model