ODEs and parameters adapted to conform to the draft SETAC report vers…

…ion 2

* The original parameter BM_min was removed
* The original parameter BM_threshold was renamed to BM_min
* The biomass ODE was modified to handle the lower biomass density
  threshold (i.e if BM(t) < BM_min) in a consistent manner, see report
* The sample scenarios were modified to reflect the changes in
  input parameters

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: lemna
 Title: Lemna Ecotox Effect Model
-Version: 0.9.1
+Version: 0.9.2
 Authors@R: c(
     person("Nils", "Kehrein", , "[email protected]", c("aut", "cre")),
     person("SETAC Europe IG Effect Modeling", role="ccp")

R/lemna_ode.R

@@ -118,11 +118,6 @@ lemna_ode <- function(t, state, param) {
     ## Population growth
     ##
 
-    # If population size falls below threshold value, set BM=BM_min (Box 1)
-    if(BM < BM_threshold) {
-      BM <- BM_min
-    }
-
     # Growth model ODE (Box 1)
     dBM <- (k_photo_max * f_photo - k_loss * f_loss) * BM
     # avoid biomass decrease below BM_min

R/param.R

@@ -8,8 +8,7 @@
 #' - `k_photo_fixed`, Model switch for unlimited growth conditions (TRUE/FALSE)
 #' - `k_photo_max`, Maximum photosynthesis rate (d-1)
 #' - `k_loss`, Reference loss rate (d-1)
-#' - `BM_threshold`, Lower biomass abundance threshold (g dw m-2)
-#' - `BM_min`, Reservoir for biomass recovery (g dw m-2)
+#' - `BM_min`, Threshold density for setting dBM/dt to zero (g dw m-2)
 #'
 #' ### Temperature response of photosynthesis
 #' - `T_opt`, Optimum growth temperature (deg C)
@@ -80,8 +79,7 @@ param_defaults <- function(values) {
     k_photo_fixed = FALSE, # model switch, if TRUE then f_loss = 1 and f_photo = fCint_photo()
     k_photo_max = 0.47,    # max photosynthesis rate (d-1)
     k_loss = 0.05,         # reference loss rate (d-1)
-    BM_threshold = 5e-4,   # lower biomass abundance threshold (g dw m-2)
-    BM_min = 0,            # reservoir for biomass recovery (g dw m-2)
+    BM_min = 5e-4,         # threshold density for setting dBM/dt to zero (g dw m-2)
 
     # response parameters
     T_opt = 26.7,   # optimum growth temperature (deg C)
@@ -132,7 +130,6 @@ param_new <- function(values) {
     k_photo_fixed = NA,
     k_photo_max = NA,
     k_loss = NA,
-    BM_threshold = NA,
     BM_min = NA,
 
     # response parameters

README.Rmd

@@ -33,11 +33,12 @@ published by Schmitt *et al.* (2013).
 ## Installation
 
 ```{r, eval=FALSE}
-# install directly from CRAN
+## install directly from CRAN
 install.packages("lemna")
 
-# install latest development version from GitHub
-devtools::install_github("nkehrein/lemna")
+## install latest development version from GitHub
+#install.packages("remotes")
+remotes::install_github("nkehrein/lemna", dependencies=TRUE)
 ```
 
 ## Usage

README.md

@@ -20,11 +20,12 @@ published by Schmitt *et al.* (2013).
 ## Installation
 
 ``` r
-# install directly from CRAN
+## install directly from CRAN
 install.packages("lemna")
 
-# install latest development version from GitHub
-devtools::install_github("nkehrein/lemna")
+## install latest development version from GitHub
+#install.packages("remotes")
+remotes::install_github("nkehrein/lemna", dependencies=TRUE)
 ```
 
 ## Usage

data-raw/focus.R

@@ -7,7 +7,7 @@ focus_param <- within(param_defaults(), {
   # variable environmental conditions
   k_photo_fixed <- FALSE
   # growth model
-  BM_threshold <- 0
+  BM_min <- 0
   BM_L <- 176
   # toxicodynamics
   b <- 4.16

data-raw/metsulfuron.R

@@ -3,7 +3,7 @@ metsulfuron_param <- within(param_defaults(), {
   # constant photosynthesis rates for backwards compatibility
   k_photo_fixed <- TRUE
   # growth model
-  BM_threshold <- 0
+  BM_min <- 0
   BM_L <- 176
   # toxicodynamics
   b <- 4.16

lemna.Rproj

@@ -1,7 +1,7 @@
 Version: 1.0
 
-RestoreWorkspace: No
-SaveWorkspace: No
+RestoreWorkspace: Default
+SaveWorkspace: Default
 AlwaysSaveHistory: Default
 
 EnableCodeIndexing: Yes

src/lemna.c

@@ -13,7 +13,7 @@
 /**
  * Allocate memory for global parameter array
  */
-static double parms[25];
+static double parms[24];
 /**
  * Allocate memory for forcing function data
  *
@@ -37,31 +37,30 @@ static double forc[5];
 #define k_photo_fixed parms[0]
 #define k_photo_max   parms[1]
 #define k_loss        parms[2]
-#define BM_threshold  parms[3]
-#define BM_min        parms[4]
+#define BM_min        parms[3]
 // response parameters
-#define T_opt         parms[5]
-#define T_min         parms[6]
-#define T_max         parms[7]
-#define Q10           parms[8]
-#define T_ref         parms[9]
-#define alpha         parms[10]
-#define beta          parms[11]
-#define N_50          parms[12]
-#define P_50          parms[13]
-#define BM_L          parms[14]
+#define T_opt         parms[4]
+#define T_min         parms[5]
+#define T_max         parms[6]
+#define Q10           parms[7]
+#define T_ref         parms[8]
+#define alpha         parms[9]
+#define beta          parms[10]
+#define N_50          parms[11]
+#define P_50          parms[12]
+#define BM_L          parms[13]
 // toxicodynamic parameters
-#define E_max         parms[15]
-#define EC50_int      parms[16]
-#define b             parms[17]
+#define E_max         parms[14]
+#define EC50_int      parms[15]
+#define b             parms[16]
 // toxicokinetic parameters
-#define P             parms[18]
-#define r_A_DW        parms[19]
-#define r_FW_DW       parms[20]
-#define r_FW_V        parms[21]
-#define r_DW_FN       parms[22]
-#define K_pw          parms[23]
-#define k_met         parms[24]
+#define P             parms[17]
+#define r_A_DW        parms[18]
+#define r_FW_DW       parms[19]
+#define r_FW_V        parms[20]
+#define r_DW_FN       parms[21]
+#define K_pw          parms[22]
+#define k_met         parms[23]
 // forcings by environmental variables
 #define C_ext forc[0]
 #define Tmp forc[1]
@@ -75,7 +74,7 @@ static double forc[5];
  */
 void lemna_init(void (* odeparms)(int *, double *))
 {
-  int N=25;
+  int N=24;
   odeparms(&N, parms);
 
 }
@@ -213,7 +212,7 @@ void lemna_func(int *neq, double *t, double *y, double *ydot, double *yout, int*
   //
 
   // Growth model ODE (Box 1)
-  dBM = (k_photo_max * f_photo - k_loss * f_loss) * (BM >= BM_threshold ? BM : BM_min);
+  dBM = (k_photo_max * f_photo - k_loss * f_loss) * BM;
   // avoid biomass decrease below BM_min
   if(BM <= BM_min && dBM < 0) {
     dBM = 0;