Title: | Calibration, Validation, and Simulation of TKTD Models |
---|---|
Description: | Eases the use of ecotoxicological effect models. Can simulate common toxicokinetic-toxicodynamic (TK/TD) models such as General Unified Threshold models of Survival (GUTS) and Lemna. It can derive effects and effect profiles (EPx) from scenarios. It supports the use of 'tidyr' workflows employing the pipe symbol. Time-consuming tasks can be parallelized. |
Authors: | Nils Kehrein [aut, cre], Dirk Nickisch [aut], Peter Vermeiren [aut], Torben Wittwer [ctb], Johannes Witt [ctb], Andre Gergs [ctb] |
Maintainer: | Nils Kehrein <[email protected]> |
License: | GPL (>= 3) |
Version: | 1.3.1 |
Built: | 2024-11-25 14:56:11 UTC |
Source: | CRAN |
The model is a mechanistic combined toxicokinetic-toxicodynamic (TK/TD) and growth model for algae. It follows the concept of a simplified algae model described in Rendal et al. (2023). The model simulates the development of algal biomass. The growth of the algae population is simulated on the basis of growth rates, which are, in contrast to the Weber model, independent on environmental conditions which are usually optimal in laboratory effect studies. The toxicodynamic sub-model describes the effects of growth-inhibiting substances through a corresponding reduction in the photosynthesis rate on the basis of either external or internal concentrations (depending on user choice of 'scaled' parameter setting).
Algae_Simple()
Algae_Simple()
an S4 object of type AlgaeSimple
The model has two state variables:
A
, Biomass (ug fresh wt/mL, cells/mL *10^4)
Dw
, only used if scaled = 1
Growth model
mu_max
, Maximum growth rate (d-1)
Concentration response (Toxicodynamics)
EC_50
, Effect concentration of 50% inhibition of growth rate (ug L-1)
b
, slope of concentration effect curve at EC_50 (-)
dose_response
, shape of the dose response curve (0 = logit, 1 = probit)
External concentration (Toxicokinetics)
kD
, dominant rate constant of toxicant in aquatic environments (d-1)
scaled
, 0 = no internal scaled damage / 1 = yes (-)
Simplified model without additional forcings for e.g. irradiation or temperature
as implemented in Algae_Weber
. A constant growth over time is assumed.
In case that growth is time dependent, a forcing variable (f_growth) can be set.
Forcing time-series are represented by data.frame
objects consisting of two
columns. The first for time and the second for a scaling factor of mu_max.
The input format for all forcings is a list of the data frames. If f_growth is
not set, a default scaling factor of 1 is used.
Upper and lower parameter boundaries are set by default for each parameter. This, to avoid extreme values during calibration (particularly likelihood profiling)
Simulation results will contain the state variables biomass (A
) and
scaled damage concentration (Dw
).
It is possible to amend the output of simulate()
with additional model
quantities that are not state variables, for e.g. debugging purposes or to
analyze model behavior. To enable or disable additional outputs, use the
optional argument nout
of simulate()
. As an example, set nout=2
to
enable reporting of external concentration (Cw
) and growth scaling factor
(f_growth
). Set nout=0
to disable additional outputs (default).
The available output levels are as follows:
nout >= 1
: Cw
external concentration (ug L-1)
nout >= 2
: f_growth
growth scaling factor (-)
nout >= 3
: dA
, biomass derivative (µg)
nout >= 4
: dDw
, damage concentration derivative (ug L-1)
Weber D, Schaeffer D, Dorgerloh M, Bruns E, Goerlitz G, Hammel K, Preuss TG and Ratte HT, 2012. Combination of a higher-tier flow-through system and population modeling to assess the effects of time-variable exposure of isoproturon on the green algae Desmodesmus subspictatus and Pseudokirchneriella subcapitata. Environmental Toxicology and Chemistry, 31, 899-908. doi:10.1002/etc.1765
Other algae models:
Algae-models
,
Algae_TKTD()
,
Algae_Weber()
The model is a mechanistic combined toxicokinetic-toxicodynamic (TK/TD) and
growth model for algae. The model simulates the development of algal biomass
under laboratory and environmental conditions. The growth of the algae
population is simulated on the basis of growth rates, which are dependent on
environmental conditions (radiation, temperature and phosphorus).
The model is a variant of the Algae_Weber()
model (Weber 2012) as cited
in EFSA TKTD opinion (2018). This Algae model, Algae_TKTD()
, provides an
additional possibility (probit) to simulate the dose-response curve and
considers a scaled internal damage instead of the external concentration.
Algae_TKTD()
Algae_TKTD()
an S4 object of type AlgaeTKTD
The model has four state variables:
A
, Biomass (ug fresh wt/mL, cells/mL *10^4)
Q
, Mass of phosphorous internal (ug P/ug fresh wt)
P
, Mass of phosphorous external (ug P/L)
Dw
, Damage concentration (ug/L)
Growth model
mu_max
, Maximum growth rate (d-1)
Q_min
, Minimum intracellular P (ug P/ug fresh wt)
Q_max
, Maximum intracellular P (ug P/ug fresh wt)
v_max
, Maximum P-uptake rate at non-limited growth (ug P/ug fresh wt/d)
k_s
, Half-saturation constant for extracellular P (mg P/L)
m_max
, Natural mortality rate (1/d)
I_opt
, Optimum light intensity for growth (uE/m²/s)
T_opt
, Optimum temperature for growth (°C)
T_max
, Maximum temperature for growth (°C)
T_min
, Minimum temperature for growth (°C)
Concentration response (Toxicodynamics)
EC_50
, Effect concentration of 50% inhibition of growth rate (ug L-1)
b
, slope of concentration effect curve at EC_50 (-)
dose_resp
, shape of the dose response curve (0 = logit, 1 = probit)
External concentration (Toxicokinetics)
kD
, dominant rate constant (d-1)
Besides exposure events (Cw), the Algae model requires two environmental
properties as time-series input: Irradiance (I
, uE/m²/s) and
temperature (T_act
, deg C).
Forcings time-series are represented by data.frame
objects
consisting of two columns. The first for time and the second for the
environmental factor in question. The input format for all forcings is a
list of the data frames.
Simulation results will contain the state variables Biomass (A
), mass of
internal phosphorous (Q
), mass of external phosphorous (P
) and the damage
concentration (Dw
).
It is possible to amend the output of simulate()
with additional model
quantities that are not state variables, for e.g. debugging purposes or to
analyze model behavior. To enable or disable additional outputs, use the
optional argument nout
of simulate()
. As an example, set nout=2
to
enable reporting of model derivatives dA
and dQ
. Set nout=0
to disable
additional outputs (default).
The available output levels are as follows:
Derivatives
nout >= 1
: dA
, biomass derivative (µg)
nout >= 2
: dQ
, internal phosphorous derivative (mg P/ug fresh wt)
nout >= 3
: dP
, external phosphorous derivative (mg P L-1)
nout >= 4
: dDw
, damage concentration derivative (ug L-1)
Weber D, Schaeffer D, Dorgerloh M, Bruns E, Goerlitz G, Hammel K, Preuss TG and Ratte HT, 2012. Combination of a higher-tier flow-through system and population modeling to assess the effects of time-variable exposure of isoproturon on the green algae Desmodesmus subspictatus and Pseudokirchneriella subcapitata. Environmental Toxicology and Chemistry, 31, 899-908. doi:10.1002/etc.1765
Other algae models:
Algae-models
,
Algae_Simple()
,
Algae_Weber()
The model is a mechanistic combined toxicokinetic-toxicodynamic (TK/TD) and growth model for algae. The model simulates the development of algal biomass under laboratory and environmental conditions and was developed by Weber et al. (2012) as cited in EFSA TKTD opinion (2018). The growth of the algae population is simulated on the basis of growth rates, which are dependent on environmental conditions (radiation, temperature and phosphorus). The toxicodynamic sub-model describes the effects of growth-inhibiting substances through a corresponding reduction in the photosynthesis rate on the basis of internal concentrations. (the implementation of Weber et al. (2012) is followed where units differ with EFSA)
Algae_Weber()
Algae_Weber()
an S4 object of type AlgaeWeber
The model has four state variables:
A
, Biomass (ug fresh wt/mL, cells/mL *10^4)
Q
, Mass of phosphorous internal (mg P/L, or ug P/mL)
P
, Mass of phosphorous external (mg P/L, or ug P/mL)
C
, external substance concentration (ug/L)
Growth model
mu_max
, Maximum growth rate (d-1)
Q_min
, Minimum intracellular P (ug P/ug fresh wt)
Q_max
, Maximum intracellular P (ug P/ug fresh wt)
v_max
, Maximum P-uptake rate at non-limited growth (ug P/ug fresh wt/d)
k_s
, Half-saturation constant for extracellular P (mg P/L)
m_max
, Natural mortality rate (1/d)
I_opt
, Optimum light intensity for growth (uE/m²/s)
T_opt
, Optimum temperature for growth (°C)
T_max
, Maximum temperature for growth (°C)
T_min
, Minimum temperature for growth (°C)
D
, Dilution rate (1/d)
R_0
, Influx concentration of P (mg P/L)
Concentration response (Toxicodynamics)
EC_50
, Effect concentration of 50% inhibition of growth rate (ug/L)
b
, slope of concentration effect curve at EC_50 (-)
External concentration (Toxicokinetics)
k
, Degradation rate of toxicant in aquatic environments (d-1)
Besides exposure events (C_in), the Algae model requires three environmental
properties as time-series input: Irradiance (I
, uE/m²/s)
and temperature (T_act
, deg C).
Forcings time-series are represented by data.frame
objects
consisting of two columns. The first for time and the second for the
environmental factor in question. The input format for all forcings is a
list of the data frames.
Simulation results will contain the state variables Biomass (A
), mass of
internal phosphorous (Q
), mass of external phosphorous (P
) and the external
concentration (C
).
It is possible to amend the output of simulate()
with additional model
quantities that are not state variables, for e.g. debugging purposes or to
analyze model behavior. To enable or disable additional outputs, use the
optional argument nout
of simulate()
. As an example, set nout=2
to
enable reporting of model derivatives dA
and dQ
. Set nout=0
to disable
additional outputs (default).
The available output levels are as follows:
Derivatives
nout >= 1
: dA
, biomass derivative (µg)
nout >= 2
: dQ
, internal phosphorous derivative (mg P/ug fresh wt)
nout >= 3
: dP
, external phosphorous derivative (mg P L-1)
nout >= 4
: dC
, external concentration derivative (ug L-1)
Default values for parameter boundaries are set for all parameters by expert judgement, for calibration purposes. Values can be access from the object, and defaults overwritten.
Weber D, Schaeffer D, Dorgerloh M, Bruns E, Goerlitz G, Hammel K, Preuss TG and Ratte HT, 2012. Combination of a higher-tier flow-through system and population modeling to assess the effects of time-variable exposure of isoproturon on the green algae Desmodesmus subspictatus and Pseudokirchneriella subcapitata. Environmental Toxicology and Chemistry, 31, 899-908. doi:10.1002/etc.1765
EFSA PPR Panel (EFSA Panel on Plant Protection Products and their Residues), Ockleford C, Adriaanse P, Berny P, Brock T, Duquesne S, Grilli S, Hernandez-Jerez AF, Bennekou SH,Klein M, Kuhl T, Laskowski R, Machera K, Pelkonen O, Pieper S, Smith RH, Stemmer M, Sundh I, Tiktak A,Topping CJ, Wolterink G, Cedergreen N, Charles S, Focks A, Reed M, Arena M, Ippolito A, Byers H andTeodorovic I, 2018. Scientific Opinion on the state of the art of Toxicokinetic/Toxicodynamic (TKTD)effect models for regulatory risk assessment of pesticides for aquatic organisms. EFSA Journal, 16(8), 5377. doi:10.2903/j.efsa.2018.5377
Other algae models:
Algae-models
,
Algae_Simple()
,
Algae_TKTD()
Overview of supported Algae models
Algae_Weber()
by Weber et al. (2012)
Algae_TKTD()
based on Weber et al. (2012), but with scaled damage
Algae_Simple()
simplified growth model without additional forcing variables
Models supporting biomass transfer can be instructed to move a fixed amount of biomass to a new medium after a period of time. This feature replicates a procedure occurring in e.g. Lemna effect studies and may be necessary to recreate study results.
The biomass transfer feature assumes that always a fixed amount of
biomass is transferred. Transfers can occur at any fixed point in time or
in regular intervals. During a transfer, the biomass is reset to the
transferred amount and additional compartments can be scaled 1:1 accordingly,
to e.g. reflect the change in internal toxicant mass when biomass is modified.
Transfer settings can be modified using set_transfer()
.
Any transfer time point must also be an output time point. If a transfer occurs, simulation results of that time point will report the model state before the transfer. Be aware that in order to use transfers at regular intervals, the simulation must start at time point zero.
Weber D, Schaeffer D, Dorgerloh M, Bruns E, Goerlitz G, Hammel K, Preuss TG and Ratte HT, 2012. Combination of a higher-tier flow-through system and population modeling to assess the effects of time-variable exposure of isoproturon on the green algae Desmodesmus subspictatus and Pseudokirchneriella subcapitata. Environmental Toxicology and Chemistry, 31, 899-908. doi:10.1002/etc.1765
EFSA Panel on Plant Protection Products and their Residues, 2018. Scientific opinion on the state of the art of Toxicokinetic/Toxicodynamic (TKTD) effect models for regulatory risk assessment of pesticides for aquatic organisms. EFSA journal 16:5377 doi:10.2903/j.efsa.2018.5377
Other algae models:
Algae_Simple()
,
Algae_TKTD()
,
Algae_Weber()
Other models:
DEB-models
,
GUTS-RED-models
,
Lemna-models
,
Macrophyte-models
,
Myriophyllum-models
Species parameters were collected from the AddMyPet database entry on Americamysis bahia (Opossum shrimp). The exposure series consists of a constant exposure resulting in medium effects on length and reproduction.
americamysis
americamysis
An object of class DebAbj
of length 1.
Cache control simulations
cache_controls(x, windows, skipZeroExposure = FALSE, ...)
cache_controls(x, windows, skipZeroExposure = FALSE, ...)
x |
vector of scenario objects |
windows |
|
skipZeroExposure |
|
... |
additional parameters passed on to |
Modified scenario objects
The function calibrate()
performs the calibration (fitting) of model
parameters to observed data. The data can originate from one or more experiments or
trials. Experimental conditions, such as model parameters and exposure
level, can differ between trials; fitting can be performed on all datasets
at the same time.
calibrate(x, ...) ## S4 method for signature 'EffectScenario' calibrate( x, par, data, endpoint = deprecated(), output, by, metric_fun = deprecated(), err_fun, as_tibble = deprecated(), catch_errors = deprecated(), verbose = TRUE, ... ) ## S4 method for signature 'CalibrationSet' calibrate(x, par, output, err_fun, verbose = TRUE, ...) ## S4 method for signature 'list' calibrate( x, par, endpoint = deprecated(), output, metric_fun = deprecated(), metric_total = deprecated(), err_fun, as_tibble = deprecated(), catch_errors = deprecated(), verbose = TRUE, ... )
calibrate(x, ...) ## S4 method for signature 'EffectScenario' calibrate( x, par, data, endpoint = deprecated(), output, by, metric_fun = deprecated(), err_fun, as_tibble = deprecated(), catch_errors = deprecated(), verbose = TRUE, ... ) ## S4 method for signature 'CalibrationSet' calibrate(x, par, output, err_fun, verbose = TRUE, ...) ## S4 method for signature 'list' calibrate( x, par, endpoint = deprecated(), output, metric_fun = deprecated(), metric_total = deprecated(), err_fun, as_tibble = deprecated(), catch_errors = deprecated(), verbose = TRUE, ... )
x |
either a single scenario or a list of caliset objects to be fitted |
... |
additional parameters passed on to |
par |
named numeric vector with parameters to fit and their start values |
data |
|
endpoint |
deprecated |
output |
|
by |
optional |
metric_fun |
deprecated, please use |
err_fun |
vectorized error function to calculate an error term that is minimized during optimization, must accept exactly four vectorized arguments, defaults to sum of squared errors |
as_tibble |
deprecated, result can no longer be returned as a tibble |
catch_errors |
deprecated, simulation errors are always caught |
verbose |
|
metric_total |
deprecated |
Fitting of model parameters can be performed in two ways:
A single scenario is fitted to a single dataset. The dataset must represent a time-series of an output variable of the model, e.g. observed biomass over time (effect data). The dataset can represent results of one or more experimental replicates under identical conditions.
One or more datasets of observed data are fitted each to a scenario which describes the experimental conditions during observation, such as exposure level and environmental properties. Each combination of dataset and scenario is represented by a calibration set. During fitting, all calibration sets are evaluated and a total error term is calculated over all observed and predicted values.
Experimental, or effect, data must be supplied as a data.frame
in long format
with at least two columns: the first column contains numeric
timestamps and
the remaining columns must contain the observed quantity. The dataset must
contain a column that which matches with the contents of parameter output
.
As an example, the simulation result of Lemna_Schmitt model contains the
output column biomass (BM
), amongst others. To fit model parameters of said
Lemna_Schmitt scenario based on observed biomass, the observed data must
contain a column named BM
which represents the observed biomass.
A minimal observed dataset could look like this:
observed <- data.frame(time=c(0, 7, 14, 21), BM=c( 12, 23, 37, 56))
By default, the total sum of squared errors is used as the target
function which is minimized during fitting. A custom error function can be
supplied by the user: The function must accept four vectorized
arguments and return a numeric of length one, i.e. the total error value
which gets minimized by calibrate()
.
Example of a custom error function which returns the sum of absolute errors:
my_absolute_error <- function(observed, predicted, weights, tags) { sum(abs(observed - predicted)) }
The arguments to the error function will contain all observed and predicted
values, as well as any weights and tags that were defined by the calibration sets.
As tags are optional, the fourth argument may be a list containing NULL
values.
The fourth argument can be used to pass additional information to the error
function: For example, the tag may identify the study from where the data
originates from and the error function could group and evaluate the data
accordingly.
A list of fitted parameters (as produced by stats::optim()
)
is returned.
calibrate(EffectScenario)
: Fit single scenario using a dataset
calibrate(CalibrationSet)
: Fit using a CalibrationSet
calibrate(list)
: Fit using a list of caliset objects
library(dplyr) # Get observed biomass during control experiment by Schmitt et al. (2013) observed <- Schmitt2013 %>% filter(ID == "T0") %>% select(t, BM=obs) # Create a scenario that represents conditions during experiment scenario <- metsulfuron %>% set_param(c(k_phot_fix=TRUE, k_resp=0, Emax=1)) %>% set_init(c(BM=12)) %>% set_noexposure() # Fit parameter 'k_phot_max' to observed biomass growth from experiment calibrate( scenario, par=c(k_phot_max=1), data=observed, output="BM", method="Brent", # Brent is recommended for one-dimensional optimization lower=0, # lower parameter boundary upper=0.5 # upper parameter boundary ) -> fit fit$par
library(dplyr) # Get observed biomass during control experiment by Schmitt et al. (2013) observed <- Schmitt2013 %>% filter(ID == "T0") %>% select(t, BM=obs) # Create a scenario that represents conditions during experiment scenario <- metsulfuron %>% set_param(c(k_phot_fix=TRUE, k_resp=0, Emax=1)) %>% set_init(c(BM=12)) %>% set_noexposure() # Fit parameter 'k_phot_max' to observed biomass growth from experiment calibrate( scenario, par=c(k_phot_max=1), data=observed, output="BM", method="Brent", # Brent is recommended for one-dimensional optimization lower=0, # lower parameter boundary upper=0.5 # upper parameter boundary ) -> fit fit$par
A calibration set combines a scenario, observed data, and an
optional weighting factor into one object. The calibration set is used to fit
model parameters to observed data using calibrate()
.
caliset(scenario, data, weight = 1.0, tag = NULL)
caliset(scenario, data, weight = 1.0, tag = NULL)
scenario |
a scenario describing conditions during the experiment |
data |
a |
weight |
optional |
tag |
optional value to identify the data, e.g. a study number |
A calibration set usually represents a single experiment or trial.
Multiple experimental replicates can be combined into a single set, if model
parameters are identical between trials.
If model parameters were modified during a trial, e.g. a pump failure occurred
or flow rates changed, this can be represented by using a scenario sequence
instead of a basic scenario. Please refer to sequence()
for details.
An optional weighting factor can be used to scale the error term of a
whole set or of individual data points when fitting parameters using e.g.
calibrate()
.
The vector of weights must either be of length one or have the same length as the dataset. In the former case, the same weight will be applied to all values in the dataset. In the latter, individual weights are applied for each data point.
caliset()
returns a calibration set object
library(dplyr) # Get observed biomass during control experiment by Schmitt et al. (2013) observed <- Schmitt2013 %>% filter(ID == "T0") %>% select(t, BM=obs) # Create a scenario that represents conditions during experiment scenario <- metsulfuron %>% set_param(c(k_phot_fix=TRUE, k_resp=0, Emax=1)) %>% set_init(c(BM=12)) %>% set_noexposure() # Create a calibration set cs <- caliset(scenario, observed) # Fit parameter 'k_phot_max' to observed biomass growth from experiment calibrate( cs, par=c(k_phot_max=1), output="BM", method="Brent", # Brent is recommended for one-dimensional optimization lower=0, # lower parameter boundary upper=0.5 # upper parameter boundary ) -> fit fit$par
library(dplyr) # Get observed biomass during control experiment by Schmitt et al. (2013) observed <- Schmitt2013 %>% filter(ID == "T0") %>% select(t, BM=obs) # Create a scenario that represents conditions during experiment scenario <- metsulfuron %>% set_param(c(k_phot_fix=TRUE, k_resp=0, Emax=1)) %>% set_init(c(BM=12)) %>% set_noexposure() # Create a calibration set cs <- caliset(scenario, observed) # Fit parameter 'k_phot_max' to observed biomass growth from experiment calibrate( cs, par=c(k_phot_max=1), output="BM", method="Brent", # Brent is recommended for one-dimensional optimization lower=0, # lower parameter boundary upper=0.5 # upper parameter boundary ) -> fit fit$par
Creates a DEB abj scenario. The abj model with type M acceleration is
like model std, but acceleration occurs between birth and metamorphosis (V1-morph).
Isomorphy is assumed before and after acceleration. Metamorphosis is before
puberty and occurs at maturity E_Hj
, which might or might not correspond with
changes in morphology. The abj model is a one-parameter extension of model std
(DEB Wiki).
DEB_abj()
DEB_abj()
The following list describes the default names and standard units of the model's state variables:
L
, structural length (cm)
E
, energy reserve (J)
H
, energy invested in maturity (J)
R
, reproduction buffer (J)
cV
, internal concentration (C)
Lmax
, maximum structural length (cm)
All state variables are initialized with zero. See set_init()
on how to set
the initial state.
The following model parameters are required:
p_M
, vol-spec somatic maintenance (J/d.cm^3)
v
, energy conductance (cm/d)
k_J
, maturity maint rate coefficient (1/d)
p_Am
, surface-area specific maximum assimilation rate (J/d.cm^2)
kap
, allocation fraction to soma (-)
E_G
, spec cost for structure (J/cm^3)
f
, scaled functional response (-)
E_Hj
, maturity at metamorphosis (J)
E_Hp
, maturity at puberty (J)
kap_R
, reproduction efficiency (-)
L_b
, structural length at birth (cm)
L_j
, structural length at metamorphosis (cm)
ke
, elimination rate constant (d-1)
c0
, no-effect concentration sub-lethal (C)
cT
, tolerance concentration (C)
MoA
, mode of action switch (-)
Any combination of the following mode of actions (MoA) can be considered by the model:
MoA = 1
: effect on feeding
MoA = 2
: effect on maintenance costs
MoA = 4
: effect on overhead costs for making an egg
MoA = 8
: hazard during oogenesis
MoA = 16
: energy conductance
To activate more than one MoA, simply add up the corresponding
codes. To disable all MoAs, set the parameter to zero.
See also set_mode_of_action()
.
The state variables L (structural length) and R (reproduction buffer) are
set as effect endpoints by default. All state variables are available as
potential endpoints. The list of considered endpoints can be modified
by using set_endpoints()
.
To calculate effects, each DEB scenario is simulated twice: One simulation
which considers exposure to a toxicant and one simulation without exposure, i.e.
a control. See also effect()
.
an S4 object of type DebAbj
Simulation results will contain the state variables.
It is possible to amend the output of simulate()
with additional model
quantities that are not state variables, for e.g. debugging purposes or to
analyze model behavior. To enable or disable additional outputs, use the
optional argument nout
of simulate()
. As an example, set nout=2
to
enable reporting of the acceleration factor (MV
) and the mobilization flux
(pC
). Set nout=0
to disable additional outputs (default).
The available output levels are as follows:
nout
>= 1: MV
acceleration factor (-)
nout
>= 2: pC
mobilization flux (J/d)
nout
>= 3: pA
assimilation flux (J/d)
nout
>= 4: pJ
energy invested in maturity flux (J/d)
Other DEB models:
DEB-models
,
DEBtox()
# Create an abj scenario from scratch and simulate it DEB_abj() %>% set_init(c(L=0.02,E=0.1,H=0.01)) %>% set_param(c(p_M=3000,v=0.02,k_J=0.6,p_Am=300,kap=0.9,E_G=4000,f=1, E_Hj=0.05,E_Hp=0.3,kap_R=0.9,ke=1,c0=0,cT=1,L_b=0.02, L_j=0.04,MoA=0)) %>% set_exposure(no_exposure()) %>% set_times(0:10) %>% simulate() # Print information about sample scenario 'americamysis' americamysis # Simulate 'americamysis' scenario americamysis %>% simulate()
# Create an abj scenario from scratch and simulate it DEB_abj() %>% set_init(c(L=0.02,E=0.1,H=0.01)) %>% set_param(c(p_M=3000,v=0.02,k_J=0.6,p_Am=300,kap=0.9,E_G=4000,f=1, E_Hj=0.05,E_Hp=0.3,kap_R=0.9,ke=1,c0=0,cT=1,L_b=0.02, L_j=0.04,MoA=0)) %>% set_exposure(no_exposure()) %>% set_times(0:10) %>% simulate() # Print information about sample scenario 'americamysis' americamysis # Simulate 'americamysis' scenario americamysis %>% simulate()
Supported models:
Other DEB models:
DEB_abj()
,
DEBtox()
Other models:
Algae-models
,
GUTS-RED-models
,
Lemna-models
,
Macrophyte-models
,
Myriophyllum-models
Creates a DEBtox scenario as described by Jager (2020). It represents a simplified DEBtox model based on DEBkiss. In the BYOM application [link], this model is referred to as DEBtox 2019, version 4.7. It supports an optional feature of the ERA special model variant, which can consider a reference Lm parameter to compare results of multiple datasets.
DEBtox() DEB_Daphnia()
DEBtox() DEB_Daphnia()
The following list describes the names and units of the model's state variables:
D
, scaled damage ([C])
L
, body length (mm)
R
, cumulative reproduction (-)
S
, survival probability (-)
State variables D
, L
, and R
are initialized with zero. Variable S
is initialized with one (1.0
). See set_init()
on how to set
the initial state manually.
The following parameters are required:
General
L0
, body length at start (mm)
Lp
, body length at puberty (mm)
Lm
, maximum body length (mm)
rB
, von Bertalanffy growth rate constant (1/d)
Rm
, maximum reproduction rate (#/d)
f
, scaled functional response (-)
hb
, background hazard rate (d-1)
a
, Weibull background hazard coefficient (-). Set to 1
to disable.
Extra parameters
Lf
, body length at half-saturation feeding (mm)
Lj
, body length at which acceleration stops (mm)
Tlag
, lag time for start development (d)
TK/TD parameters
kd
, dominant rate constant (d-1)
zb
, effect threshold energy budget ([C])
bb
, effect strength energy-budget effects (1/[C])
zs
, effect threshold survival ([C])
bs
, effect strength survival (1/([C] d))
Other parameters (formerly globals in BYOM)
FBV
, dry weight egg as fraction of structural body weight (-)
KRV
, part. coeff. repro buffer and structure (kg/kg) (for losses with reproduction)
kap
, approximation for kappa (for starvation response)
yP
, product of yVA and yAV (for starvation response)
Lm_ref
, optional reference max length for scaling rate constants (mm).
Set to zero to disable the reference length. Disabled by default.
len
, a switch to control body length dynamics: 1
organism can shrink,
2
organism cannot shrink. Default value is 1
.
Tbp
, optional brood-pouch delay (d). Set to NA
or zero to disable.
Default value is 0
.
MoA
, mode of action switches (-). Default value is 0
.
FB
, feedback on damage dynamics switches (-). Default value is 0
.
A reference Lm_ref
is needed to properly compare different data sets,
or when calibrating on more than one data set. If Lm
differs, one would not
want to have different rate constants at the same length.
Any combination of the following mode of actions (MoA) can be considered by the model:
MoA = 1
: assimilation/feeding
MoA = 2
: costs for maintenance
MoA = 4
: costs for growth and reproduction
MoA = 8
: costs for reproduction
MoA = 16
: hazard for reproduction
To activate more than one mode of action, simply add up the corresponding
codes and set parameter MoA
to the desired value. To disable all mode of
actions, set parameter MoA
to zero. See also set_moa()
.
As an example, to consider effects on feeding and maintenance, set the
mode of action to three (3
):
DEBtox() %>% set_param(c(MoA=3))
Any combination of the following damage feedbacks can be considered by the model:
1
: surf:vol scaling uptake rate
2
: surf:vol scaling elimination rate
4
: growth dilution
8
: losses with reproduction
To activate more than one feedback, simply add up the corresponding codes. To disable all feedbacks, set the parameter to zero.
The state variables L (body length), R (cumulative reproduction), and
S (survival probability) are set as effect endpoints by default. All state
variables are available as potential endpoints. The list of considered
endpoints can be modified by using set_endpoints()
.
To calculate effects, each DEBtox scenario is simulated twice: One simulation
which considers exposure to a toxicant and one simulation without exposure, i.e.
a control. See also effect()
.
The following intermediary model variables can be added to the model
output on demand. Simply set the optional parameter nout
to the
required output level and pass it to simulate()
.
nout
>= 1: f
, actual scaled response
nout
>= 2: fR
, actual f considering starvation
nout
>= 3: kd
, actual kd
nout
>= 4: s
, stress level
nout
>= 5: h
, hazard rate
nout
>= 6: sA
, stress factor on assimilation/feeding
nout
>= 7: sM
, stress factor on maintenance
nout
>= 8: sG
, stress factor on growth costs
nout
>= 9: sR
, stress factor on reproduction costs
nout
>= 10: sH
, stress factor on hazard to reproduction
nout
>= 11: xu
, damage feedback factor for surf:vol scaling uptake rate
nout
>= 12: xe
, damage feedback factor for surf:vol scaling elimination rate
nout
>= 13: xG
, damage feedback factor for growth dilution
nout
>= 14: xR
, damage feedback factor for losses with repro
cvasi v1.0.0
The DEB_Daphnia()
model implemented BYOM's DEBtox 2019
model version 4.5
cvasi v1.2.0
The model equations were updated to conform with BYOM's DEBtox 2019
version 4.7. This introduced a new model parameter a
, the Weibull
background hazard coefficient, and limited the maximum hazard rate to
99% per hour.
The scenario constructor was renamed to DEBtox()
.
Additional intermediary model variables available as optional simulation output
an S4 object of type DebTox
DEB_Daphnia()
: Deprecated model variant of DEBtox()
Jager T, 2020: Revisiting simplified DEBtox models for analysing ecotoxicity data. Ecol Model 416. doi:10.1016/j.ecolmodel.2019.108904
Romoli et al., 2024: Environmental risk assessment with energy budget models: a comparison between two models of different complexity. Environ Toxicol Chem 43(2):440-449. doi:10.1002/etc.5795
Other DEB models:
DEB-models
,
DEB_abj()
Species and substance parameters were collected from test runs of the original DEBtox Daphnia model.
dmagna
dmagna
An object of class DebTox
of length 1.
Returns a data.frame
with points on the dose response curve for the given
effect scenario.
dose_response( scenario, range = c(1, 99), n = 20, strategy = c("exponential", "decadic", "vanilla"), verbose = FALSE, ... )
dose_response( scenario, range = c(1, 99), n = 20, strategy = c("exponential", "decadic", "vanilla"), verbose = FALSE, ... )
scenario |
used for calculation |
range |
numeric vector specifying the required range of effect levels in
percent (%), defaults to |
n |
minimum number of points on the dose response curve |
strategy |
controls how multiplication factors are chosen, |
verbose |
logical, set to |
... |
additional arguments passed on to |
Derives a dose response curve from a scenario. The result will
cover the requested range of effect levels. The tested multiplication factors
can be chosen by different strategies, i.e. a vanilla
approach using a
fixed set of factors, or decadic
and exponential
approaches
employing logarithmic and exponential factor scaling, respectively.
data.frame
with two columns, i.e. mf
and effect
# basic dose response curve minnow_sd %>% dose_response() # modify the minimum number of points on the curve minnow_sd %>% dose_response(n=10) # select a subset of the effect range minnow_sd %>% dose_response(range=c(10,20)) # use an alternative strategy for the selection of multiplication factors minnow_sd %>% dose_response(strategy="decadic") # provide additional output how multiplication factors were selected minnow_sd %>% dose_response(verbose=TRUE)
# basic dose response curve minnow_sd %>% dose_response() # modify the minimum number of points on the curve minnow_sd %>% dose_response(n=10) # select a subset of the effect range minnow_sd %>% dose_response(range=c(10,20)) # use an alternative strategy for the selection of multiplication factors minnow_sd %>% dose_response(strategy="decadic") # provide additional output how multiplication factors were selected minnow_sd %>% dose_response(verbose=TRUE)
Derives the effect level due to toxicant exposure in the supplied scenarios. Either relative to a control scenario or derived directly from model endpoints, depending on model type. For scenarios with moving exposure windows, the maximum effect is returned.
effect(x, factor = 1, max_only = TRUE, ep_only = FALSE, marginal_effect, ...)
effect(x, factor = 1, max_only = TRUE, ep_only = FALSE, marginal_effect, ...)
x |
vector of |
factor |
optional numeric value which scales the exposure time-series |
max_only |
|
ep_only |
logical, if TRUE only effect endpoints are returned as a vector |
marginal_effect |
|
... |
additional parameters passed on to |
By default, only the maximum effect in all moving exposure windows will
be returned. If argument max_only=FALSE
is set, the returned table will
be converted to long-format and will contain effect levels for each
assessed exposure window.
Start and end time of exposure windows can be disabled by setting ep_only=TRUE
.
Effect levels smaller than a certain threshold can be automatically set to
zero (0.0
) to avoid spurious effect levels introduced by numerical errors.
Set marginal_effect
to an adequate value less than 1%.
Calculations can be sped up by providing a data.frame
of pre-calculated
control scenarios for each assessed time window. As control scenarios are
by definition independent of any exposure multiplication factor, they can
be reused for repeated calculations, e.g. to derive effect profiles or
dose-response relationships.
a tibble
, by default containing scenarios, effect levels, and the
exposure window where the maximum effect level occurred. The number of columns
depends on the enabled effect endpoints and function arguments.
By default, the first column, named scenarios
, contains the original scenario
objects that were the basis of the calculation. For each effect endpoint, it
will be followed by one column with the maximum effect level and two columns
containing start and end time of the associated exposure window. If exposure
windows are disabled, the columns will just contain the start and end time of
the simulation. The effect level column will have the name of the effect
endpoint, start and end time will additionally have the suffixes .dat.start
and .dat.end
, respectively.
Derives one or more EPx/LPx values for the supplied effect scenarios, i.e. it calculates the multiplication factors of an exposure profile that cause x% of effect. Scenarios are processed in parallel, if possible.
epx( scenarios, level = c(10, 50), effect_tolerance = 0.001, factor_cutoff = NA, min_factor = 1e-30, max_factor = 1e+30, verbose = FALSE, ep_only = FALSE, long_format = FALSE, ... )
epx( scenarios, level = c(10, 50), effect_tolerance = 0.001, factor_cutoff = NA, min_factor = 1e-30, max_factor = 1e+30, verbose = FALSE, ep_only = FALSE, long_format = FALSE, ... )
scenarios |
table or vector of |
level |
effect levels in percent (%), defaults to |
effect_tolerance |
|
factor_cutoff |
optional |
min_factor |
|
max_factor |
|
verbose |
|
ep_only |
|
long_format |
|
... |
additional arguments passed on to |
To estimate EPx values, a binary search on multiplication factors is conducted.
The algorithm can achieve arbitrary precision in terms of effects. The
same approach is implemented in the morse
package in the MFx()
function.
Convergence is often achieved in less than 10 iterations per effect level and
endpoint.
Internally, a knowledge base of all tried factors and resulting effect levels is
kept to speed up convergence if more than one endpoint or effect level was
requested. The algorithm will automatically sweep the range of multiplication
factors as needed but hard cutoff values are implemented to avoid infinite loops;
the algorithm will halt with an error message if tried factors are
smaller than 1e-30
or greater than 1e30
.
The precision of reported EPx values is controlled by the argument
effect_tolerance
and is given as the upper absolute error threshold of
effects that is deemed acceptable. The default value of 0.001
ensures that
a derived EPx will result in an effect of x% ± 0.1. Decreasing the
effect_tolerance
will result in additional model iterations and longer
runtime. Setting an extremely small tolerance value may lead to a breakdown
of the algorithm due to the occurrence of extremely small, quasi-random
numerical errors in simulation results.
The original tibble
with additional columns named after the request effect levels, e.g. L.EP10.
If no tibble was used as argument, then a new one is created. The first column scenario
will contain
the supplied EffectScenario
objects.
minnow_sd %>% epx() minnow_sd %>% epx(level=c(10,23,42)) # displays infos about tested multiplication factors minnow_sd %>% epx(verbose=TRUE) # return results as a table in wide format minnow_sd %>% epx(long_format=TRUE)
minnow_sd %>% epx() minnow_sd %>% epx(level=c(10,23,42)) # displays infos about tested multiplication factors minnow_sd %>% epx(verbose=TRUE) # return results as a table in wide format minnow_sd %>% epx(long_format=TRUE)
Calls epx()
to calculate the EPx value (i.e. the multiplication factors of
an exposure profile that cause x% of effect) for moving windows with length
window_length
that move timesteps defined by window_interval
.
epx_mtw( x, level = c(10, 50), factor_cutoff = 1000, window_length = 7, window_interval = 1, ... )
epx_mtw( x, level = c(10, 50), factor_cutoff = 1000, window_length = 7, window_interval = 1, ... )
x |
a scenario |
level |
The target effect level of the effect, ie. the x of EPx. |
factor_cutoff |
above which cutoff is the EPx is not relevant |
window_length |
the length of the moving time window |
window_interval |
the interval that the moving time window moves |
... |
arguments passed to |
a tibble with five columns
window.start
window.end
endpoint
level
EPx
metsulfuron %>% set_window(length=7, interval=1) %>% epx_mtw()
metsulfuron %>% set_window(length=7, interval=1) %>% epx_mtw()
The function is aimed at getting an idea of how the parameter space
of a model behaves, so that parameter identifiability problems and correlations
between parameters can be explored. Therefore, the function samples a large
number of parameter sets by randomly drawing from each parameter's 95%
confidence interval (generated by lik_profile()
). It then
checks how many of the parameter sets are within acceptable limits by comparing
the likelihood ratio of a parameter set vs. the original parameter set against
a chi-square distribution as degrees of freedom (df) the total number of profile
parameters (outer rim) or one df (inner rim). If needed, the function resamples
until at least nr_accept
parameters sets are within the inner rim
explore_space( x, par, res, output, sample_size = 1000, max_runs = 30, nr_accept = 100, sample_factor = 1.2 )
explore_space( x, par, res, output, sample_size = 1000, max_runs = 30, nr_accept = 100, sample_factor = 1.2 )
x |
a list of caliset objects |
par |
best fit parameters from joined calibration |
res |
output of 'lik_profile()' function |
output |
character vector, name of output column of |
sample_size |
number of samples to draw from each parameter interval |
max_runs |
max number of times to redraw samples (within a smaller space), and repeat the process |
nr_accept |
threshold for number of points sampled within the inner circle |
sample_factor |
multiplication factor for sampling (95% interval * sample factor) |
a list containing a plot to explore the parameter space, and the data.frame
supporting it
library(dplyr) # Example with Lemna model - physiological params # Before applying the function, a model needs to be calibrated and its parameters profiled # Inputs for likelihood profiling # exposure - control run exp <- Schmitt2013 %>% filter(ID == "T0") %>% select(time=t, conc) # observations - control run obs <- Schmitt2013 %>% filter(ID == "T0") %>% select(t, BM=obs) # parameters after calibration params <- c( k_phot_max = 5.663571, k_resp = 1.938689, Topt = 26.7 ) # set parameter boundaries (if different from defaults) bounds <- list( k_resp = list(0, 10), k_phot_max = list(0, 30), Topt = list(20, 30) ) # update metsulfuron myscenario <- metsulfuron %>% set_init(c(BM = 5, E = 1, M_int = 0)) %>% set_param(list( k_0 = 5E-5, a_k = 0.25, BM50 = 17600, mass_per_frond = 0.1 )) %>% set_exposure(exp) %>% set_param(params) %>% set_bounds(bounds) # Likelihood profiling res <- lik_profile( x = myscenario, data = obs, output = "BM", par = params, refit = FALSE, type = "fine", method = "Brent" ) # plot plot_lik_profile(res) # parameter space explorer set.seed(1) # for reproducibility res_space <- explore_space( x = list(caliset(myscenario, obs)), par = params, res = res, output = "BM", sample_size = 1000, max_runs = 20, nr_accept = 100) plot_param_space(res_space)
library(dplyr) # Example with Lemna model - physiological params # Before applying the function, a model needs to be calibrated and its parameters profiled # Inputs for likelihood profiling # exposure - control run exp <- Schmitt2013 %>% filter(ID == "T0") %>% select(time=t, conc) # observations - control run obs <- Schmitt2013 %>% filter(ID == "T0") %>% select(t, BM=obs) # parameters after calibration params <- c( k_phot_max = 5.663571, k_resp = 1.938689, Topt = 26.7 ) # set parameter boundaries (if different from defaults) bounds <- list( k_resp = list(0, 10), k_phot_max = list(0, 30), Topt = list(20, 30) ) # update metsulfuron myscenario <- metsulfuron %>% set_init(c(BM = 5, E = 1, M_int = 0)) %>% set_param(list( k_0 = 5E-5, a_k = 0.25, BM50 = 17600, mass_per_frond = 0.1 )) %>% set_exposure(exp) %>% set_param(params) %>% set_bounds(bounds) # Likelihood profiling res <- lik_profile( x = myscenario, data = obs, output = "BM", par = params, refit = FALSE, type = "fine", method = "Brent" ) # plot plot_lik_profile(res) # parameter space explorer set.seed(1) # for reproducibility res_space <- explore_space( x = list(caliset(myscenario, obs)), par = params, res = res, output = "BM", sample_size = 1000, max_runs = 20, nr_accept = 100) plot_param_space(res_space)
Creates an object that encapsulates an exposure time-series with its
metadata, such as formatted datetime strings and file name where the
series was loaded from. no_exposure()
is shorthand to create a time-series
of constant zero exposure.
ExposureSeries(series, dates, file, meta, context)
ExposureSeries(series, dates, file, meta, context)
series |
|
dates |
|
file |
|
meta |
|
context |
|
an S4 object of type ExposureSeries
dates
original time points of time-series, e.g. time stamps of the form 2000-01-01 12:00
file
character
, file name where data originates from, may be empty
meta
list
, contains metadata
context
list
, contains contextual metadata, such as project ids
series
data.frame
containing the actual time-series
A mechanistic combined toxicokinetic-toxicodynamic (TK/TD) and growth model for the aquatic macrophytes Lemna spp. as published by Klein et al. (2021).
focusd1
focusd1
An object of class LemnaSetac
of length 1.
The scenario will simulate a period of 365 days, a start population of 80 g/m² dry weight, variable environmental conditions, and a complex, time-varying exposure pattern.
The scenario setup was published by Hommen et al. (2015). Exposure pattern and substance specific parameters are of exemplary character and represent the herbicide metsulfuron-methyl. The parameters were derived by Schmitt et al. (2013) based on literature data.
Hommen U., Schmitt W., Heine S., Brock Theo CM., 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/effect-modeling.html
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
# Simulate the example scenario focusd1 %>% simulate()
# Simulate the example scenario focusd1 %>% simulate()
Generic to calculate effects for a particular scenario
fx(scenario, ...) ## S4 method for signature 'ANY' fx(scenario, ...) ## S4 method for signature 'Algae' fx(scenario, ...) ## S4 method for signature 'GutsRedSd' fx(scenario, ...) ## S4 method for signature 'GutsRedIt' fx(scenario, ...) ## S4 method for signature 'Lemna' fx(scenario, ...) ## S4 method for signature 'Myriophyllum' fx(scenario, ...)
fx(scenario, ...) ## S4 method for signature 'ANY' fx(scenario, ...) ## S4 method for signature 'Algae' fx(scenario, ...) ## S4 method for signature 'GutsRedSd' fx(scenario, ...) ## S4 method for signature 'GutsRedIt' fx(scenario, ...) ## S4 method for signature 'Lemna' fx(scenario, ...) ## S4 method for signature 'Myriophyllum' fx(scenario, ...)
scenario |
scenario object |
... |
additional parameters |
numeric named vector
fx(ANY)
: Use state variables at end of simulation
fx(Algae)
: Effect at end of simulation of Algae-models
fx(GutsRedSd)
: Survival and lethality in GUTS-RED-models
fx(GutsRedIt)
: Survival and lethality in GUTS-RED-models
fx(Lemna)
: Effect at end of simulation of Lemna-models
fx(Myriophyllum)
: Effect at end of simulation of Myriophyllum-models
Returns the unique model name that is associated with a scenario, e.g.
GUTS-RED-IT
. The function supports vectorized arguments.
get_model(x)
get_model(x)
x |
(vector of) scenarios or |
vector of character
# returns `GUTS-RED-IT` get_model(minnow_it)
# returns `GUTS-RED-IT` get_model(minnow_it)
Returns the user-defined, custom tag of a scenario, if available. Tags can be helpful to quickly distinguish scenarios of the same model type. The function supports vectorized inputs.
get_tag(x)
get_tag(x)
x |
(vector of) scenarios or |
vector of character
# returns `fathead minnow` get_tag(minnow_it) # update or set a tag myscenario <- minnow_it %>% set_tag("My Custom Tag") # returns `My Custom Tag` get_tag(myscenario)
# returns `fathead minnow` get_tag(minnow_it) # update or set a tag myscenario <- minnow_it %>% set_tag("My Custom Tag") # returns `My Custom Tag` get_tag(myscenario)
Reduced General Unified Threshold models of Survival (GUTS) with individual tolerance (IT).
GUTS_RED_IT(param, init)
GUTS_RED_IT(param, init)
param |
optional named |
init |
optional named numeric |
an S4 object of type GutsRedIt
The return value of simulate()
will contain values for the state variables,
as well as an additional column S
which represents the survival probability
for each time point. S
is calculated as described in EFSA (2018) as
S = (1- F(t)). The background hazard rate hb
is already considered in state
variable H
and therefore does not occur as an additional term to derive S
.
The following list describes the default names and standard units of GUTS-RED state variables:
D
, scaled damage (conc)
H
, cumulative hazard (-)
The state variables are initialized with zero by default.
kd
, dominant rate constant (time^-1)
hb
, background hazard rate (time^-1)
alpha
, median of thresholds (conc)
beta
, shape parameter (-)
The effect endpoint L
(lethality) is available for GUTS-RED models.
A value of zero (0.0
) denotes no effect on organism survival. A value of
one (1.0
) denotes a lethality rate of 100%, i.e. no survivors.
The survival probability S
is available in the return value of simulate()
.
EFSA PPR Panel (EFSA Panel on Plant Protection Products and their Residues), Ockleford C, Adriaanse P, Berny P, et al., 2018: Scientific Opinion on the state of the art of Toxicokinetic/Toxicodynamic (TKTD) effect models for regulatory risk assessment of pesticides for aquatic organisms. EFSA Journal 2018; 16(8):5377, 188 pp. doi:10.2903/j.efsa.2018.5377
Other GUTS-RED models:
GUTS-RED-models
,
GUTS_RED_SD()
Reduced General Unified Threshold models of Survival (GUTS) with stochastic death (SD).
GUTS_RED_SD(param, init)
GUTS_RED_SD(param, init)
param |
optional named |
init |
optional named numeric |
an S4 object of type GutsRedSd
The return value of simulate()
will contain values for the state variables,
as well as an additional column S
which represents the survival probability
for each time point. S
is calculated as described in EFSA (2018) as
S = exp(-H). The background hazard rate hb
is already considered in state
variable H
and therefore does not occur as an additional term to derive S
.
The following list describes the default names and standard units of GUTS-RED state variables:
D
, scaled damage (conc)
H
, cumulative hazard (-)
The state variables are initialized with zero by default.
kd
, dominant rate constant (time^-1)
hb
, background hazard rate (time^-1)
z
, threshold for effects (conc)
kk
, killing rate constant (time^-1)
The effect endpoint L
(lethality) is available for GUTS-RED models.
A value of zero (0.0
) denotes no effect on organism survival. A value of
one (1.0
) denotes a lethality rate of 100%, i.e. no survivors.
The survival probability S
is available in the return value of simulate()
.
EFSA PPR Panel (EFSA Panel on Plant Protection Products and their Residues), Ockleford C, Adriaanse P, Berny P, et al., 2018: Scientific Opinion on the state of the art of Toxicokinetic/Toxicodynamic (TKTD) effect models for regulatory risk assessment of pesticides for aquatic organisms. EFSA Journal 2018; 16(8):5377, 188 pp. doi:10.2903/j.efsa.2018.5377
Other GUTS-RED models:
GUTS-RED-models
,
GUTS_RED_IT()
Reduced General Unified Threshold models of Survival (GUTS) with stochastic death (SD) and individual tolerance (IT)
The TKTD models GUTS-RED-SD and GUTS-RED-IT were described by EFSA (2018). GUTS-RED models assume a one-compartment model which directly links external concentration to the scaled damage. The scaled damage is given in units of concentration, equal to the units of measurement in the external medium, e.g. ug/L. The damage dynamics is connected to an individual hazard state variable, resulting in simulated mortality when an internal damage threshold is exceeded. The death mechanisms stochastic death (SD) and individual threshold (IT) are extreme cases of the GUTS theory.
For SD models, the threshold parameter for lethal effects is fixed and identical for all individuals of a group, meaning that the variance of the threshold values is zero. Hence, the killing rate relates the probability of a mortality event in proportion to the scaled damage. For IT models, the thresholds for effects are distributed among individuals of a group. Mortality of an individual follows immediately once the individual's tolerance is exceeded. Meaning in model terms that the killing rate is set to infinity (EFSA 2018).
The following list describes the default names and standard units of GUTS-RED state variables:
D
, scaled damage (conc)
H
, cumulative hazard (-)
The state variables are initialized with zero by default.
kd
, dominant rate constant (time^-1)
hb
, background hazard rate (time^-1)
z
, threshold for effects (conc)
kk
, killing rate constant (time^-1)
kd
, dominant rate constant (time^-1)
hb
, background hazard rate (time^-1)
alpha
, median of thresholds (conc)
beta
, shape parameter (-)
The effect endpoint L
(lethality) is available for GUTS-RED models.
A value of zero (0.0
) denotes no effect on organism survival. A value of
one (1.0
) denotes a lethality rate of 100%, i.e. no survivors.
The survival probability S
is available in the return value of simulate()
.
EFSA PPR Panel (EFSA Panel on Plant Protection Products and their Residues), Ockleford C, Adriaanse P, Berny P, et al., 2018: Scientific Opinion on the state of the art of Toxicokinetic/Toxicodynamic (TKTD) effect models for regulatory risk assessment of pesticides for aquatic organisms. EFSA Journal 2018; 16(8):5377, 188 pp. doi:10.2903/j.efsa.2018.5377
Other GUTS-RED models:
GUTS_RED_IT()
,
GUTS_RED_SD()
Other models:
Algae-models
,
DEB-models
,
Lemna-models
,
Macrophyte-models
,
Myriophyllum-models
morse
model parametersLoads GUTS model parameters which were fitted by the morse package.
import_morse( fit, find_sd = TRUE, find_it = TRUE, reset_hb = TRUE, params = c("estim", "all"), mcmc_size, find.SD = deprecated(), find.IT = deprecated(), reset.hb = deprecated(), mcmc.size = deprecated(), file = deprecated() ) morse(...)
import_morse( fit, find_sd = TRUE, find_it = TRUE, reset_hb = TRUE, params = c("estim", "all"), mcmc_size, find.SD = deprecated(), find.IT = deprecated(), reset.hb = deprecated(), mcmc.size = deprecated(), file = deprecated() ) morse(...)
fit |
Either a string with a file path to an .Rdata or .RDS file containing a morse fit, or a morse fit object itself |
find_sd |
a logical value. If |
find_it |
a logical value. If |
reset_hb |
a logical value. If |
params |
|
mcmc_size |
optional |
find.SD |
deprecated |
find.IT |
deprecated |
reset.hb |
deprecated |
mcmc.size |
deprecated |
file |
deprecated |
... |
Arguments passed on to |
list of parameter_set objects
morse()
: deprecated alias
# import all parameter fits try(import_morse("path/to/morse_fit.RData")) # import parameters for a specific model try(import_morse("path/to/morse_fit.RData", find_it=TRUE, find_sd=FALSE)) # modify model objects try(models %>% set_param(import_morse("path/to/morse_fit.RData")))
# import all parameter fits try(import_morse("path/to/morse_fit.RData")) # import parameters for a specific model try(import_morse("path/to/morse_fit.RData", find_it=TRUE, find_sd=FALSE)) # modify model objects try(models %>% set_param(import_morse("path/to/morse_fit.RData")))
Read all TOXSWA files within a SWASH project directory.
import_swash(swash_dir, ...)
import_swash(swash_dir, ...)
swash_dir |
path to the SWASH project directory |
... |
arguments passed on to |
a list of imported exposure series, see import_toxswa()
for details
Read one or more TOXSWA exposure series from TOXSWA's .out
files. By default, the
concentration dissolved in water (ConLiqWatLay) at the end of the simulated waterbody
(i.e. at the maximum of the x dimension) is returned. The unit of the time
scale as well as of the imported model output variable can be scaled as needed.
import_toxswa( files, alias = NA, output_var = "ConLiqWatLay", output_unit = "ug/L", time_unit = "days", substance = NULL, split = TRUE )
import_toxswa( files, alias = NA, output_var = "ConLiqWatLay", output_unit = "ug/L", time_unit = "days", substance = NULL, split = TRUE )
files |
vector of strings with absolute or relative paths to files |
alias |
optional vector with strings, will be used as an alias to identify a TOXSWA series instead of its filename |
output_var |
character, single output variable from TOXSWA that is imported, defaults to ConLiqWatLay |
output_unit |
character, target unit of the imported output variable,
defaults to ug/L, syntax must be compatible with |
time_unit |
character, target unit of the imported time scale,
defaults to days, syntax must be compatible with |
substance |
optional vector of characters, if set, only the substance codes defined in this vector are imported |
split |
logical, if |
The numerical time scale is shifted to always start at time zero (0.0).
Numerical columns of the returned data.frame objects will be of type
units::units. Please be aware that the use of units
objects may not be
supported by all functions in this package. However, set_times()
and
set_exposure()
can handle units
objects safely.
Incomplete list of alternative TOXSWA v5.5.3 output variables:
ConLiqWatLay: Concentration dissolved in water (g/m3)
ConLiqSed: Concentration in pore water sediment (g/m3)
ConSysWatLay: Total concentration in water (g/m3)
CntSorSusSol: Content sorbed to suspended solids (g/kg)
CntSorSed: Content sorbed to sediment (g/kg)
list of data.frame objects with exposure series. Each data.frame has at least three columns:
time
: numerical time scale, always starts at zero
timestamp
: time as datetime objects such as POSIXct
one or more additional columns for each imported substance
Test if argument is a DEB model
is_DEB(x)
is_DEB(x)
x |
vector of |
vector of logical
values
Test if argument is a GUTS model
is_GUTS(x) is_GUTS_IT(x) is_GUTS_SD(x)
is_GUTS(x) is_GUTS_IT(x) is_GUTS_SD(x)
x |
vector of |
vector of logical
values
is_GUTS_IT()
: Test if argument is a GUTS-IT model
is_GUTS_SD()
: Test if argument is a GUTS-IT model
# returns `TRUE` is_GUTS(minnow_it) is_GUTS(GUTS_RED_IT()) # returns `c(TRUE,TRUE,TRUE)` is_GUTS(c(minnow_it, minnow_it, minnow_it)) # returns `FALSE` is_GUTS_SD(minnow_it)
# returns `TRUE` is_GUTS(minnow_it) is_GUTS(GUTS_RED_IT()) # returns `c(TRUE,TRUE,TRUE)` is_GUTS(c(minnow_it, minnow_it, minnow_it)) # returns `FALSE` is_GUTS_SD(minnow_it)
Also returns TRUE
for LemnaThreshold
models
is_Lemna(x)
is_Lemna(x)
x |
vector of scenarios objects |
vector of logical values
Test if argument is a LemnaThreshold model
is_LemnaThreshold(x)
is_LemnaThreshold(x)
x |
vector of scenarios objects |
vector of logical
values
Supports vectorized arguments.
is_scenario(x)
is_scenario(x)
x |
Some value or object |
vector of logical
values
# returns `TRUE` is_scenario(minnow_it) # returns `FALSE` is_scenario(list())
# returns `TRUE` is_scenario(minnow_it) # returns `FALSE` is_scenario(list())
The model is a mechanistic combined toxicokinetic-toxicodynamic (TK/TD) and growth model for the aquatic macrophytes Lemna spp. The model simulates the development of Lemna biomass under laboratory and environmental conditions and was developed by Schmitt et al. (2013). Growth of the Lemna population is simulated on basis of photosynthesis and respiration rates which are functions of environmental conditions. The toxicodynamic sub-model describes the effects of growth-inhibiting substances by a respective reduction in the photosynthesis rate based on internal concentrations. This is the historical version of the Lemna model. For current uses, we recommend the Lemna (SETAC) model, which is a more recent version of the Schmitt model.
Lemna_Schmitt(param, init) Lemna_SchmittThold(param, init)
Lemna_Schmitt(param, init) Lemna_SchmittThold(param, init)
param |
optional named |
init |
optional named numeric |
Constructors to ease creation of scenarios based on the Lemna model by
Schmitt et al. (2013).
A variant of this Lemna model, Lemna_SchmittThold()
, provides an additional
cumulative exposure threshold parameter. The Lemna biomass stops growing
if the integral of exposure over time exceeds the threshold. The integral
of exposure is internally accounted for by an additional state variable
AUC
(Area Under Curve).
an S4 object of type LemnaSchmitt
Lemna_SchmittThold()
: model variant with cumulative exposure threshold
The following list describes the default names and standard units of the model's state variables:
BM, g_dw/m2, dry weight biomass per square meter
E, -, effect [0,1]
M_int, ug, internal toxicant mass
AUC, ug/L, cumulative exposure (only for LemnaThreshold
model)
Biomass (BM) and internal toxicant mass (M_int) are initialized to zero by
default. See set_init()
on how to set the initial states.
The following model parameters are required:
Fate and biomass
k_phot_fix
, logical, TRUE then k_phot_max is not changed by environmental factors, else FALSE
k_phot_max
, 1/d, maximum photosynthesis rate
k_resp
, 1/d, respiration rate
k_loss
, 1/d, rate of loss (e.g. flow rate)
mass_per_frond
, g_dw/frond, dry weight per frond
BMw2BMd
, g_fw/g_dw, Fresh weight/dry weight
Effect
Emax
, -, maximum effect [0,1]
EC50
, ug/L, midpoint of effect curve
b
, -, slope of effect curve
Toxicokinetics
P_up
, cm/d, Permeability for uptake
AperBM
, cm2/g_dw, A_leaf / d_leaf = 1/d_leaf (for circular disc, d=0.05 cm)
Kbm
, -, Biomass(fw) : water partition coefficient
P_Temp
, logical, TRUE to enable temperature dependence of cuticle permeability, else FALSE
MolWeight
, g/mol, Molmass of molecule (determines Q10_permeability)
Temperature dependence
Tmin
, deg C, minimum temperature for growth
Tmax
, deg C, maximum temperature for growth
Topt
, deg C, optimal temperature for growth
t_ref
, deg C, reference temperature for respiration rate
Q10
, -, temperature dependence factor for respiration rate
Light dependence
k_0
, 1/d, light dependence: intercept of linear part
a_k
, (1/d)/(kJ/m2.d), light dependence: slope of linear part
Phosphorus dependence (Hill like dep.)
C_P
, mg/L, phosphorus concentration in water
CP50
, mg/L, phosphorus conc. where growth rate is halfed
a_p
, -, Hill coefficient
KiP
, mg/L, p-inhibition constant for very high p-conc.
Nitrogen dependence (Hill like dep.)
C_N
, mg/L, nitrogen concentration in water
CN50
, mg/L, n-conc. where growth rate is halfed
a_N
, -, Hill coefficient
KiN
, mg/L, n-inhibition constant for very high p-conc.
Density dependence
BM50
, g_dw/m2, cut off BM
The Lemna_SchmittThold
model requires the following additional parameter:
threshold
, ug/L, cumulative exposure threshold
Besides exposure, the Lemna model requires two environmental properties as
time-series input: global radiation (rad
, kJ/m2.d) and temperature (temp
, deg C).
Forcings time-series are represented by data.frame
objects consisting of two
columns. The first for time and the second for the environmental factor in question.
Entries of the data.frame
need to be ordered chronologically. A time-series
can consist of only a single row; in this case it will represent constant
environmental conditions. See scenarios for more details.
Supported effect endpoints include BM (biomass) and r (average growth rate during simulation). The effect on biomass is calculated from the last state of a simulation. Be aware that endpoint r is incompatible with frond transfers.
Default values for parameter boundaries are set for all parameters by expert judgement, for calibration purposes. Values can be access from the object, and defaults overwritten.
Simulation results will contain two additional columns besides state variables:
C_int
, ug/L, internal concentration of toxicant
FrondNo
, -, number of fronds
It is possible to amend the output of simulate()
with additional model
quantities that are not state variables, for e.g. debugging purposes or to
analyze model behavior. To enable or disable additional outputs, use the
optional argument nout
of simulate()
, see examples below. nout=1
enables reporting of internal concentration (C_int), nout=14
enables all
additional outputs, and nout=0
will disable additional outputs.
The available output levels are as follows:
nout >= 1
: C_int
, internal concentration (ug/L)
nout >= 2
: FrondNo
, number of fronds (-)
nout >= 3
: C_int_u
, unbound internal concentration (ug/l)
Growth and TK/TD
nout >= 4
: BM_fresh
, fresh weight biomass (g_fw/m2)
nout >= 5
: k_photo_eff
, current photosynthesis rate (1/d)
nout >= 6
: k_resp_eff
, current respiration rate (1/d)
nout >= 7
: f_Eff
, toxic effect factor (-)
nout >= 8
: P_up_eff
, current permeability for uptake (cm/d)
Environmental variables
nout >= 9
: actConc
, current toxicant concentration in surrounding medium (ug/L)
nout >= 10
: actTemp
, current environmental temperature (deg C)
nout >= 11
: actRad
, current environmental radiation (kJ/m2.d)
Derivatives
nout >= 12
: d BM/dt
, current change in state variable BM
nout >= 13
: d E/dt
, current change in effect
nout >= 14
: d M_int/dt
, current change in internal toxicant mass
Models supporting biomass transfer can be instructed to move a fixed amount of biomass to a new medium after a period of time. This feature replicates a procedure occurring in e.g. Lemna effect studies and may be necessary to recreate study results.
The biomass transfer feature assumes that always a fixed amount of
biomass is transferred. Transfers can occur at any fixed point in time or
in regular intervals. During a transfer, the biomass is reset to the
transferred amount and additional compartments can be scaled 1:1 accordingly,
to e.g. reflect the change in internal toxicant mass when biomass is modified.
Transfer settings can be modified using set_transfer()
.
Any transfer time point must also be an output time point. If a transfer occurs, simulation results of that time point will report the model state before the transfer. Be aware that in order to use transfers at regular intervals, the simulation must start at time point zero.
Schmitt W., Bruns E., Dollinger M., and 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
Lemna-models, Macrophyte-models, Transferable, Scenarios
Other Lemna models:
Lemna-models
,
Lemna_SETAC()
Other macrophyte models:
Lemna_SETAC()
,
Macrophyte-models
,
Myrio()
,
Myrio_log()
The model was described and published by the SETAC Europe Interest Group Effect Modeling (Klein et al. 2022). It is based on the Lemna model by Schmitt (2013). The model is a mechanistic combined toxicokinetic-toxicodynamic (TK/TD) and growth model for the aquatic macrophytes Lemna spp.. The model simulates the development of Lemna biomass under laboratory and environmental conditions. Growth of the Lemna population is simulated on basis of photosynthesis and respiration rates which are functions of environmental conditions. The toxicodynamic sub-model describes the effects of growth-inhibiting substances by a respective reduction in the photosynthesis rate based on internal concentrations.
Lemna_SETAC()
Lemna_SETAC()
an S4 object of type LemnaSetac
The model has two state variables:
BM
, Biomass (g dw m-2)
M_int
, Mass of toxicant in plant population (mass per m2, e.g. ug m-2)
Growth model
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)
Temperature response of photosynthesis
T_opt
, Optimum growth temperature (°C)
T_min
, Minimum growth temperature (°C)
T_max
, Maximum growth temperature (°C)
Temperature response of biomass loss rate
Q10
, Temperature coefficient (-)
T_ref
, Reference temperature for response=1 (°C)
Irradiance reponse of photosynthesis
alpha
, Slope of irradiance response (m2 d kJ-1)
beta
, Intercept of irradiance response (-)
Nutrient response of photosynthesis
N_50
, Half-saturation constant of Nitrogen (mg N L-1)
P_50
, Half-saturation constant of Phosphorus (mg P L-1)
Density dependence of photosynthesis
BM_L
, Carrying capacity (g dw m-2)
Concentration response (Toxicodynamics)
EC50_int
, Internal concentration resulting in 50% effect (ug L-1)
E_max
, Maximum inhibition (-)
b
, Slope parameter (-)
Internal concentration (Toxicokinetics)
P
, Permeability (cm d-1)
r_A_DW
, Area per dry-weight ratio (cm2 g-1)
r_FW_DW
, Fresh weight per dry weight ratio (-)
r_FW_V
, Fresh weight density (g cm-3)
r_DW_FN
, Dry weight per frond ratio (g dw)
K_pw
, Partitioning coefficient plant:water (-)
k_met
, Metabolisation rate (d-1)
Besides exposure, the model requires four environmental properties as time-series input:
tmp
, temperature (°C)
irr
, irradiance (kJ m-2 d-1)
P
, Phosphorus concentration (mg P L-1)
N
, Nitrogen concentration (mg N L-1)
Forcings time-series are represented by data.frame
objects consisting of two
columns. The first for time and the second for the environmental factor in question.
Entries of the data.frame
need to be ordered chronologically. A time-series
can consist of only a single row; in this case it will represent constant
environmental conditions. See scenarios for more details.
Supported effect endpoints include BM (biomass) and r (average growth rate during simulation). The effect on biomass is calculated from the last state of a simulation. Be aware that endpoint r is incompatible with biomass transfers.
For reasons of convenience, the return value contains by default two additional
variables derived from simulation results: the internal concentration C_int
as well as the number of fronds FrondNo
. These can be disabled by setting
the argument nout = 0
.
The available output levels are as follows:
nout >= 1
: C_int
, internal concentration (mass per volume)
nout >= 2
: FrondNo
, frond number (-)
Response functions
nout >= 3
: f_loss
, respiration dependency function (-)
nout >= 4
: f_photo
, photosynthesis dependency function (-)
nout >= 5
: fT_photo
, temperature response of photosynthesis (-)
nout >= 6
: fI_photo
, irradiance response of photosynthesis (-)
nout >= 7
: fP_photo
, phosphorus response of photosynthesis (-)
nout >= 8
: fN_photo
, nitrogen response of photosynthesis (-)
nout >= 9
: fBM_photo
, density response of photosynthesis (-)
nout >= 10
: fCint_photo
, concentration response of photosynthesis (-)
Environmental variables
nout >= 11
: C_int_unb
, unbound internal concentration (mass per volume)
nout >= 12
: C_ext
, external concentration (mass per volume)
nout >= 13
: Tmp
, temperature (deg C)
nout >= 14
: Irr
, irradiance (kJ m-2 d-1)
nout >= 15
: Phs
, Phosphorus concentration (mg P L-1)
nout >= 16
: Ntr
, Nitrogen concentration (mg N L-1)
Derivatives
nout >= 17
: dBM
, biomass derivative (g dw m-2 d-1)
nout >= 18
: dM_int
, mass of toxicant in plants derivative (mass per m2 d-1)
Models supporting biomass transfer can be instructed to move a fixed amount of biomass to a new medium after a period of time. This feature replicates a procedure occurring in e.g. Lemna effect studies and may be necessary to recreate study results.
The biomass transfer feature assumes that always a fixed amount of
biomass is transferred. Transfers can occur at any fixed point in time or
in regular intervals. During a transfer, the biomass is reset to the
transferred amount and additional compartments can be scaled 1:1 accordingly,
to e.g. reflect the change in internal toxicant mass when biomass is modified.
Transfer settings can be modified using set_transfer()
.
Any transfer time point must also be an output time point. If a transfer occurs, simulation results of that time point will report the model state before the transfer. Be aware that in order to use transfers at regular intervals, the simulation must start at time point zero.
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/effect-modeling.html
Schmitt W., Bruns E., Dollinger M., and 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
Lemna-models, Macrophyte-models, Transferable, Scenarios
Other Lemna models:
Lemna-models
,
Lemna_Schmitt()
Other macrophyte models:
Lemna_Schmitt()
,
Macrophyte-models
,
Myrio()
,
Myrio_log()
Overview of supported Lemna models
Lemna_Schmitt()
by Schmitt et al. (2013)
Lemna_SETAC()
by Klein et al. (2021)
Models supporting biomass transfer can be instructed to move a fixed amount of biomass to a new medium after a period of time. This feature replicates a procedure occurring in e.g. Lemna effect studies and may be necessary to recreate study results.
The biomass transfer feature assumes that always a fixed amount of
biomass is transferred. Transfers can occur at any fixed point in time or
in regular intervals. During a transfer, the biomass is reset to the
transferred amount and additional compartments can be scaled 1:1 accordingly,
to e.g. reflect the change in internal toxicant mass when biomass is modified.
Transfer settings can be modified using set_transfer()
.
Any transfer time point must also be an output time point. If a transfer occurs, simulation results of that time point will report the model state before the transfer. Be aware that in order to use transfers at regular intervals, the simulation must start at time point zero.
Other Lemna models:
Lemna_SETAC()
,
Lemna_Schmitt()
Other models:
Algae-models
,
DEB-models
,
GUTS-RED-models
,
Macrophyte-models
,
Myriophyllum-models
The aim of the function is two-fold: 1) estimate a 95% confidence
around each parameter of a calibrated model, and 2) see if perhaps a local minimum was found rather
than a global minimum. To achieve this, the likelihood profiling goes through
every parameter one by one. For each parameter,
the model is sequentially refit with the parameter value set to
increasingly lower and higher values, and the likelihood of the model given the
data calculated (using log_lik()
). The likelihood is then compared
to the likelihood of the original model (using a likelihood ratio). This leads
to the development of a likelihood profile, from which a plot a 95%
confidence interval for the parameter is derived.
The idea of the function is a variable stepwise algorithm:
When the likelihood ratio changes very little (less than l_crit_min
), the stepsize is
increased (up to a maximum, specified by f_step_max
). When the lik.
ratio changes too much (more than l_crit_max
), the algorithm tries again
with a smaller stepsize (also bound to a minimum: f_step_min
). Note that
the stepsize is used as a fraction of the parameter value that is tried.
To prevent very small stepsizes when the value goes towards zero (as can
be the case for effect thresholds), an absolute minimum
stepsize (f_step_abs
), which is specified as a fraction of the best
parameter value (Xhat
) (unless it is zero, then algoritm takes
something small).
The function was inspired by a MatLab BYOM v.6.8 procedure, created by Tjalling Jager. For details, please refer to BYOM (http://debtox.info/byom.html) as well as Jager (2021).
lik_profile( x, par, output, data = NULL, bounds = NULL, refit = TRUE, type = c("coarse", "fine"), break_prof = FALSE, ... )
lik_profile( x, par, output, data = NULL, bounds = NULL, refit = TRUE, type = c("coarse", "fine"), break_prof = FALSE, ... )
x |
|
par |
named vector - parameters (names and values) to be profiled |
output |
character vector, name of output column of |
data |
only needed if |
bounds |
optional list of lists (including lower and upper bound): uses defaults in |
refit |
if 'TRUE' (default), refit if a better minimum is found |
type |
"fine" or "coarse" (default) likelihood profiling |
break_prof |
if 'TRUE' (default), stop the profiling if a better optimum is located |
... |
additional parameters passed on to |
A list containing, for each parameter profiled, the likelihood profiling results as a dataframe; the 95% confidence interval; the original parameter value; the likelihood plot object; and the recalibrated parameter values (in case a lower optimum was found)
Jager T, 2021. Robust Likelihood-Based Optimization and Uncertainty Analysis of Toxicokinetic-Toxicodynamic Models. Integrated Environmental Assessment and Management 17:388-397. doi:10.1002/ieam.4333
# Example with Lemna model - physiological params library(dplyr) # exposure - control run exp <- Schmitt2013 %>% filter(ID == "T0") %>% select(time=t, conc) # observations - control run obs <- Schmitt2013 %>% filter(ID == "T0") %>% select(t, BM=obs) # update metsulfuron myscenario <- metsulfuron %>% set_param(c(k_phot_fix = TRUE,Emax = 1)) %>% set_init(c(BM = 12)) %>% set_exposure(exp) fit <- calibrate( x = myscenario, par = c(k_phot_max = 1), data = obs, output = "BM", lower=0, upper=1, method="Brent" ) # Likelihood profiling res <- lik_profile( x = myscenario, data = obs, output = "BM", par = fit$par, bounds = list( k_phot_max = list(0, 30) ), refit = FALSE, type = "fine", method = "Brent" ) # plot plot_lik_profile(res)
# Example with Lemna model - physiological params library(dplyr) # exposure - control run exp <- Schmitt2013 %>% filter(ID == "T0") %>% select(time=t, conc) # observations - control run obs <- Schmitt2013 %>% filter(ID == "T0") %>% select(t, BM=obs) # update metsulfuron myscenario <- metsulfuron %>% set_param(c(k_phot_fix = TRUE,Emax = 1)) %>% set_init(c(BM = 12)) %>% set_exposure(exp) fit <- calibrate( x = myscenario, par = c(k_phot_max = 1), data = obs, output = "BM", lower=0, upper=1, method="Brent" ) # Likelihood profiling res <- lik_profile( x = myscenario, data = obs, output = "BM", par = fit$par, bounds = list( k_phot_max = list(0, 30) ), refit = FALSE, type = "fine", method = "Brent" ) # plot plot_lik_profile(res)
Start and stop logging
log_enable(file = NULL, append = TRUE, envir = parent.frame()) log_disable()
log_enable(file = NULL, append = TRUE, envir = parent.frame()) log_disable()
file |
|
append |
|
envir |
log will be automatically disabled if |
no return value
Log R environment properties
log_envir()
log_envir()
no return value
Calculates the sum of log likelihoods of each observation given the model parameterization (considering a normal distribution around the prediction for each datapoint)
log_lik(npars, obs, pred)
log_lik(npars, obs, pred)
npars |
named numeric vector of parameters that the model was calibrated on |
obs |
numeric vector of observed values |
pred |
numeric vector of predicted values |
the log likelihood value
# observations obs <- c(12, 38, 92, 176, 176, 627, 1283, 2640) # intercept, a, and slope, b, of a Poisson regression fitted through obs pars <- c(a = 2, b = 0.73) # predictions with the Poisson regression pred <- c(15.43, 32.15, 66.99, 139.57, 290.82, 605.94, 1262.52, 2630.58) # example plot plot(seq(1:length(obs)), obs) lines(seq(1:length(obs)), pred) log_lik( npars = length(pars), obs = obs, pred = pred )
# observations obs <- c(12, 38, 92, 176, 176, 627, 1283, 2640) # intercept, a, and slope, b, of a Poisson regression fitted through obs pars <- c(a = 2, b = 0.73) # predictions with the Poisson regression pred <- c(15.43, 32.15, 66.99, 139.57, 290.82, 605.94, 1262.52, 2630.58) # example plot plot(seq(1:length(obs)), obs) lines(seq(1:length(obs)), pred) log_lik( npars = length(pars), obs = obs, pred = pred )
Message will only appear in the console or in log file if logging was
enabled using log_enable()
.
log_msg(...)
log_msg(...)
... |
elements will be concatenated using |
no return value
log_msg("this message will not appear") log_enable() log_msg("this message will appear") log_msg("a number of ","elements to ",42," concatenate")
log_msg("this message will not appear") log_enable() log_msg("this message will appear") log_msg("a number of ","elements to ",42," concatenate")
Log scenario properties
log_scenarios(x, header = TRUE)
log_scenarios(x, header = TRUE)
x |
vector of |
header |
|
unmodified argument x
Population models of standard test macrophytes, such as Lemna spp.
Available macrophyte models:
Models supporting biomass transfer can be instructed to move a fixed amount of biomass to a new medium after a period of time. This feature replicates a procedure occurring in e.g. Lemna effect studies and may be necessary to recreate study results.
The biomass transfer feature assumes that always a fixed amount of
biomass is transferred. Transfers can occur at any fixed point in time or
in regular intervals. During a transfer, the biomass is reset to the
transferred amount and additional compartments can be scaled 1:1 accordingly,
to e.g. reflect the change in internal toxicant mass when biomass is modified.
Transfer settings can be modified using set_transfer()
.
Any transfer time point must also be an output time point. If a transfer occurs, simulation results of that time point will report the model state before the transfer. Be aware that in order to use transfers at regular intervals, the simulation must start at time point zero.
Other macrophyte models:
Lemna_SETAC()
,
Lemna_Schmitt()
,
Myrio()
,
Myrio_log()
Other models:
Algae-models
,
DEB-models
,
GUTS-RED-models
,
Lemna-models
,
Myriophyllum-models
Data set for the parametrisation of a mechanistic combined toxicokinetic-toxicodynamic (TK/TD) and growth model for the aquatic macrophytes Lemna spp. as published by Schmitt et al. (2013). The growth model was parameterised by Schmitt et al. based on these data while toxicokinetic and toxicodynamic parameters were determined by calibrating the model using substance specific effect data of metsulfuron-methyl.
metsulfuron
metsulfuron
An object of class LemnaSchmittScenario
of length 1.
Schmitt W., Bruns E., Dollinger M., and 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
The scenario consists of a parameterized GUTS-RED-IT model and a constant exposure series. Model parameters were derived from an acute fish toxicity study of the fathead minnow and chlorpyrifos by Geiger et al. (1988). The dataset is also referred to as GUTS Ring-test dataset C by EFSA (2018). Fitted parameters were estimated using the morse package.
minnow_it
minnow_it
An object of class GutsRedIt
of length 1.
The background mortality rate (hb
) was set to zero.
https://mosaic.univ-lyon1.fr/guts
Geiger D.L., Call D.J., and Brooke L.T., 1988: Acute toxicities of organic chemicals to fathead minnows (Pimephales promelas): Volume IV, pp. 195-197. University of Wisconsin-Superior, Center for Lake Superior Environmental Studies. ISBN 9780961496838.
EFSA PPR Panel (EFSA Panel on Plant Protection Products and their Residues), Ockleford C, Adriaanse P, Berny P, et al., 2018: Scientific Opinion on the state of the art of Toxicokinetic/Toxicodynamic (TKTD) effect models for regulatory risk assessment of pesticides for aquatic organisms. EFSA Journal 2018; 16(8):5377, 188 pp. doi:10.2903/j.efsa.2018.5377
The scenario consists of a parameterized GUTS-RED-SD model and a constant exposure series. Model parameters were derived from an acute fish toxicity study of the fathead minnow and chlorpyrifos by Geiger et al. (1988). The dataset is also referred to as GUTS Ring-test dataset C by EFSA (2018). Fitted parameters were estimated using the morse package.
minnow_sd
minnow_sd
An object of class GutsRedSd
of length 1.
The background mortality rate (hb
) was set to zero.
https://mosaic.univ-lyon1.fr/guts
Geiger D.L., Call D.J., and Brooke L.T., 1988: Acute toxicities of organic chemicals to fathead minnows (Pimephales promelas): Volume IV, pp. 195-197. University of Wisconsin-Superior, Center for Lake Superior Environmental Studies. ISBN 9780961496838.
EFSA PPR Panel (EFSA Panel on Plant Protection Products and their Residues), Ockleford C, Adriaanse P, Berny P, et al., 2018: Scientific Opinion on the state of the art of Toxicokinetic/Toxicodynamic (TKTD) effect models for regulatory risk assessment of pesticides for aquatic organisms. EFSA Journal 2018; 16(8):5377, 188 pp. doi:10.2903/j.efsa.2018.5377
The Myriophyllum model is derived from the Lemna TKTD model by
Klein et al. (2021). The Myriophyllum model is mathematically equivalent
to the Tier 2C version of the Lemna model by Klein et al. (2021),
cf. Lemna_SETAC()
. Recommended settings for Tier 2C are k_photo_fixed=TRUE
and k_resp=0
(Klein et al. 2021).
In particular, the growth model is a simple exponential growth model,
which is considered to be the typical situation for a laboratory macrophyte
study. Instead of frond numbers as for Lemna, the biomass is also returned as
total shoot length (TSL
) in simulation results.
Consequently, the model has the additional parameter r_DW_TSL
(dry weight per total shoot length ratio) instead of r_DW_FN
(dry weight
per frond number ratio).
Myrio()
Myrio()
an S4 object of type MyrioExpScenario
The model has two state variables:
BM
, Biomass (g dw m-2 for field studies or mg dw for lab)
M_int
, Mass of toxicant in plant population (mass per m2, e.g. ug m-2)
Growth model
k_photo_max
, Maximum photosynthesis rate (d-1), default: 0.47
Concentration response (Toxicodynamics)
EC50_int
, Internal concentration resulting in 50% effect (ug L-1)
E_max
, Maximum inhibition (-), default: 1
b
, Slope parameter (-)
Internal concentration (Toxicokinetics)
P
, Permeability (cm d-1)
r_A_DW
, Area per dry-weight ratio (cm2 g-1), default: 1000
r_FW_DW
, Fresh weight per dry weight ratio (-), default: 16.7
r_FW_V
, Fresh weight density (g cm-3), default: 1
r_DW_TSL
, Dry weight per total shoot length ratio (g (field) or mg (lab) dw cm-1)
K_pw
, Partitioning coefficient plant:water (-), default: 1
k_met
, Metabolisation rate (d-1), default: 0
None.
Default values for parameter boundaries are set for all parameters by expert
judgement, for calibration purposes. Values can be modified using set_bounds()
.
Simulation results will contain the state variables
It is possible to amend the output of simulate()
with additional model
quantities that are not state variables, for e.g. debugging purposes or to
analyze model behavior. To enable or disable additional outputs, use the
optional argument nout
of simulate()
. As an example, set nout=2
to
enable reporting of the acceleration factor (MV
) and the mobilization flux
(pC
). Set nout=0
to disable additional outputs (default).
The available output levels are as follows:
nout
>= 1: C_int
, internal concentration (mass per volume)
nout
>= 2: TSL
, total shoot length (?)
nout
>= 3: f_photo
, photosynthesis dependency function (-)
Growth and TK/TD
nout
>= 4: C_int_unb
, unbound internal concentration (mass per volume)
nout
>= 5: C_ext
, external concentration (mass per volume)
Derivatives
nout
>= 6: dBM
, biomass derivative (g dw m-2 d-1)
nout
>= 7: dM_int
, mass of toxicant in plants derivative (mass per m2 d-1)
Supported effect endpoints include BM (biomass) and r (average growth rate during simulation). The effect on biomass is calculated from the last state of a simulation. Be aware that endpoint r is incompatible with biomass transfers.
Models supporting biomass transfer can be instructed to move a fixed amount of biomass to a new medium after a period of time. This feature replicates a procedure occurring in e.g. Lemna effect studies and may be necessary to recreate study results.
The biomass transfer feature assumes that always a fixed amount of
biomass is transferred. Transfers can occur at any fixed point in time or
in regular intervals. During a transfer, the biomass is reset to the
transferred amount and additional compartments can be scaled 1:1 accordingly,
to e.g. reflect the change in internal toxicant mass when biomass is modified.
Transfer settings can be modified using set_transfer()
.
Any transfer time point must also be an output time point. If a transfer occurs, simulation results of that time point will report the model state before the transfer. Be aware that in order to use transfers at regular intervals, the simulation must start at time point zero.
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/effect-modeling.html
Macrophyte-models, Transferable, Scenarios
Other Myriophyllum models:
Myrio_log()
,
Myriophyllum-models
Other macrophyte models:
Lemna_SETAC()
,
Lemna_Schmitt()
,
Macrophyte-models
,
Myrio_log()
The Myriophyllum model is derived from the Lemna TKTD model
by Klein et al. (2021).
Myrio_log()
modifies the Myrio()
model to feature logistic growth, i.e.
control growth is described by the differential equation
d BM/dt = k_photo_max*BM*(1 - BM/BM_L)
where BM_L
is the carrying capacity.
Myrio_log()
Myrio_log()
an S4 object of type MyrioLogScenario
Growth model
k_photo_max
, Maximum photosynthesis rate (d-1), default: 0.47
BM_L
, Carrying capacity (g dw m-2)
Concentration response (Toxicodynamics)
EC50_int
, Internal concentration resulting in 50% effect (ug L-1)
E_max
, Maximum inhibition (-), default: 1
b
, Slope parameter (-)
Internal concentration (Toxicokinetics)
P
, Permeability (cm d-1)
r_A_DW
, Area per dry-weight ratio (cm2 g-1), default: 1000
r_FW_DW
, Fresh weight per dry weight ratio (-), default: 16.7
r_FW_V
, Fresh weight density (g cm-3), default: 1
r_DW_TSL
, Dry weight per total shoot length ratio (?)
K_pw
, Partitioning coefficient plant:water (-), default: 1
k_met
, Metabolisation rate (d-1), default: 0
The model has two state variables:
BM
, Biomass (g dw m-2 for field studies or mg dw for lab)
M_int
, Mass of toxicant in plant population (mass per m2, e.g. ug m-2)
None.
Simulation results will contain the state variables
It is possible to amend the output of simulate()
with additional model
quantities that are not state variables, for e.g. debugging purposes or to
analyze model behavior. To enable or disable additional outputs, use the
optional argument nout
of simulate()
. As an example, set nout=2
to
enable reporting of the acceleration factor (MV
) and the mobilization flux
(pC
). Set nout=0
to disable additional outputs (default).
The available output levels are as follows:
nout
>= 1: C_int
, internal concentration (mass per volume)
nout
>= 2: TSL
, total shoot length (?)
nout
>= 3: f_photo
, photosynthesis dependency function (-)
Growth and TK/TD
nout
>= 4: C_int_unb
, unbound internal concentration (mass per volume)
nout
>= 5: C_ext
, external concentration (mass per volume)
Derivatives
nout
>= 6: dBM
, biomass derivative (g dw m-2 d-1)
nout
>= 7: dM_int
, mass of toxicant in plants derivative (mass per m2 d-1)
Supported effect endpoints include BM (biomass) and r (average growth rate during simulation). The effect on biomass is calculated from the last state of a simulation. Be aware that endpoint r is incompatible with biomass transfers.
Models supporting biomass transfer can be instructed to move a fixed amount of biomass to a new medium after a period of time. This feature replicates a procedure occurring in e.g. Lemna effect studies and may be necessary to recreate study results.
The biomass transfer feature assumes that always a fixed amount of
biomass is transferred. Transfers can occur at any fixed point in time or
in regular intervals. During a transfer, the biomass is reset to the
transferred amount and additional compartments can be scaled 1:1 accordingly,
to e.g. reflect the change in internal toxicant mass when biomass is modified.
Transfer settings can be modified using set_transfer()
.
Any transfer time point must also be an output time point. If a transfer occurs, simulation results of that time point will report the model state before the transfer. Be aware that in order to use transfers at regular intervals, the simulation must start at time point zero.
Default values for parameter boundaries are set for all parameters by expert
judgement, for calibration purposes. Values can be modified using set_bounds()
.
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/effect-modeling.html
Other Myriophyllum models:
Myrio()
,
Myriophyllum-models
Other macrophyte models:
Lemna_SETAC()
,
Lemna_Schmitt()
,
Macrophyte-models
,
Myrio()
Supported models:
Myrio()
, with exponential growth
Myrio_log()
, with logistic growth
Other Myriophyllum models:
Myrio()
,
Myrio_log()
Other models:
Algae-models
,
DEB-models
,
GUTS-RED-models
,
Lemna-models
,
Macrophyte-models
Creates an ExposureSeries with zero concentration. When setting the zero exposure, pay attention not to accidentally reset the output times of your scenario as the zero exposure series contains only a single time point. See the examples.
no_exposure()
no_exposure()
an S4 object of type ExposureSeries
# this will reset the output times of the sample scenario, # simulate() will quit with an error try( minnow_it %>% set_exposure(no_exposure()) %>% simulate() ) # set zero exposure, but keep original output times minnow_it %>% set_exposure(no_exposure(), reset_times=FALSE) %>% simulate()
# this will reset the output times of the sample scenario, # simulate() will quit with an error try( minnow_it %>% set_exposure(no_exposure()) %>% simulate() ) # set zero exposure, but keep original output times minnow_it %>% set_exposure(no_exposure(), reset_times=FALSE) %>% simulate()
Set of model parameters
parameter_set(model, param = list(), tag = NA_character_)
parameter_set(model, param = list(), tag = NA_character_)
model |
|
param |
named |
tag |
|
an S4 object of type ParameterSet
model
character
, a string containing a model name, e.g. "GUTS-RED-IT"
tag
character
, an optional identifier
param
named list
of model parameters
# create a parameter set and assign it ps <- parameter_set("GUTS-RED-IT", list(kd=0.12, hb=0.3)) GUTS_RED_IT() %>% set_param(ps) # multiple scenarios can be modified at once c(GUTS_RED_IT(), GUTS_RED_IT()) %>% set_param(ps) # model names must match, otherwise an error will be raised try(GUTS_RED_SD() %>% set_param(ps))
# create a parameter set and assign it ps <- parameter_set("GUTS-RED-IT", list(kd=0.12, hb=0.3)) GUTS_RED_IT() %>% set_param(ps) # multiple scenarios can be modified at once c(GUTS_RED_IT(), GUTS_RED_IT()) %>% set_param(ps) # model names must match, otherwise an error will be raised try(GUTS_RED_SD() %>% set_param(ps))
In certain cases it might be beneficial to disable parallel execution of e.g. effect profile calculations. By disabling, all processes run sequentially and instantly pass messages to the console which would be delayed during parallel processing. This makes it easier to pinpoint problems within the data or algorithm.
pll_debug(state = TRUE)
pll_debug(state = TRUE)
state |
|
no return value
These functions overload base::plot()
to provide simple plotting
routines to display various time-series and scenario objects.
## S3 method for class 'cvasi.drc' plot(x, y, scale_x = c("auto", "log10", "none"), ...) ## S3 method for class 'cvasi.simulate' plot(x, y, ...)
## S3 method for class 'cvasi.drc' plot(x, y, scale_x = c("auto", "log10", "none"), ...) ## S3 method for class 'cvasi.simulate' plot(x, y, ...)
x |
object to plot |
y |
unused parameter |
scale_x |
character, controls how the x-axis is scaled. |
... |
unused parameters |
ggplot2 plot object
plot(cvasi.drc)
: Plot dose response curves
plot(cvasi.simulate)
: Plot return value of simulate()
plot_epx( EPx_ts, exposure_ts, draw = TRUE, time_col = "time", conc_col = "conc", epx_x_title = "Start time", conc_y_title = "Exposure conc." )
plot_epx( EPx_ts, exposure_ts, draw = TRUE, time_col = "time", conc_col = "conc", epx_x_title = "Start time", conc_y_title = "Exposure conc." )
EPx_ts |
the result of |
exposure_ts |
an exposure time series with columns for time 't' and concentration 'conc' |
draw |
Should the whole plot be drawn? If FALSE the exposure plot and the EPx plot are returned as a list for later modification |
time_col |
the name of the time column in the exposure dataset |
conc_col |
the name of the concentration column in the exposure dataset |
epx_x_title |
title of the x-axis of the epx panel |
conc_y_title |
title of the y-axis of the concentration panel |
a grid of ggplots
ti <- 0:21 expo <- abs(0.01*ti + rnorm(length(ti), 0, 0.05)) exposure <- data.frame(time = ti, conc = expo) metsulfuron_epx_mtw <- metsulfuron %>% set_exposure(exposure) %>% epx_mtw(level = 10, factor_cutoff = 1000) metsulfuron_epx_mtw plot_epx(EPx_ts = metsulfuron_epx_mtw, exposure_ts = exposure, conc_y_title = "env. concentration [µg/L]")
ti <- 0:21 expo <- abs(0.01*ti + rnorm(length(ti), 0, 0.05)) exposure <- data.frame(time = ti, conc = expo) metsulfuron_epx_mtw <- metsulfuron %>% set_exposure(exposure) %>% epx_mtw(level = 10, factor_cutoff = 1000) metsulfuron_epx_mtw plot_epx(EPx_ts = metsulfuron_epx_mtw, exposure_ts = exposure, conc_y_title = "env. concentration [µg/L]")
The function provides a combined plot of the likelihood profiles of all parameters profiled.
plot_lik_profile(x)
plot_lik_profile(x)
x |
object of class lik_profile |
plots
The function provides bivariate parameter space plots indicating parameter draws (from the 95%confidence intervals per parameter obtained through likelihood profiling) that fall within the inner rim (in green, i.e. parameter sets which are not significantly different from the original, based on a chi-square test). The original parameter set is also indicated (in orange), and, if different from the original set, the best fit parameter set is indicated (in red)
plot_param_space(x)
plot_param_space(x)
x |
object of class param_space |
plots
plot_ppc( rs_mean, rs_range, col_number = 2, obs_mean = NULL, obs_full = NULL, xy_lim = NULL, study = NULL )
plot_ppc( rs_mean, rs_range, col_number = 2, obs_mean = NULL, obs_full = NULL, xy_lim = NULL, study = NULL )
rs_mean |
|
rs_range |
|
col_number |
column to plot, default = 2 |
obs_mean |
|
obs_full |
|
xy_lim |
optional |
study |
optional |
A sample of parameters representing the uncertainty within the dataset is passed to the function. All parameter combinations and exposure patterns are simulated and the range of predicted frond numbers is derived for a single study. The uncertainty is displayed by a Posterior Predictive Plot (PPC). The data (rs_mean, obs_mean and obs_full) must have the following format (col1 = time, col2 = data of interest, col3 = trial name). Data for uncertainties (rs_range) must have the format: col1 = time, col2 = lower boundaries, col3 = upper boundaries, col4 = trial. The user should take care of the input data and consider whether control data and data at time zero should be included in the model check.
a ggplot2 plot object
plot_ppc_combi(table, xy_lim = NULL)
plot_ppc_combi(table, xy_lim = NULL)
table |
|
xy_lim |
optional |
The function expects a data.frame with five mandatory and one optional column. The mandatory columns are as follows:
pred
: mean of predictions e.g. frond number for lemna
max
: maximum of predictions
min
: minimum of predictions
obs
: observations
PPC
: color code
The optional column is to be named study
and contains a study identifier.
If more than one study identifier is present in the table, individual
studies will be plotted in different colors and a legend will be displayed.
The function is called by plot_ppc where the column names are defined
(see rs_ppc object).
a ggplot2 plot object
plot_scenario(model_base, plot_col = 2, trial_number = NULL)
plot_scenario(model_base, plot_col = 2, trial_number = NULL)
model_base |
effect scenario object with mean parameters |
plot_col |
output column which should be plotted, default = 2 |
trial_number |
name for model run (if available tag is used) |
This function has been deprecated and replaced by the generic plot()
.
Sometimes it is helpful if the user can plot results of one effect scenario. This is for instance the case for test simulations or predictions for one profile. This function runs the simulation for one effect scenario and plots the results. Function plots the time (column 1) and the predictions (column 2, can be changed by the user plot_col)
plot of the results for one effect scenario
# Please use `plot()` instead metsulfuron %>% simulate() %>% plot()
# Please use `plot()` instead metsulfuron %>% simulate() %>% plot()
plot_sd( model_base, treatments, rs_mean, rs_range = NULL, obs_mean = NULL, obs_full = NULL, x_breaks = NULL, y_lim = NULL, grid_labels = NULL, grid_ncol = 2, plot_col = 2, y_title = NULL, ... )
plot_sd( model_base, treatments, rs_mean, rs_range = NULL, obs_mean = NULL, obs_full = NULL, x_breaks = NULL, y_lim = NULL, grid_labels = NULL, grid_ncol = 2, plot_col = 2, y_title = NULL, ... )
model_base |
effect scenario object with mean parameters |
treatments |
treatments exposure levels as data frame |
rs_mean |
|
rs_range |
|
obs_mean |
|
obs_full |
|
x_breaks |
optional vector of breaks of x-axis |
y_lim |
optional vector containing limits of y-axis |
grid_labels |
optional labels of grid headers |
grid_ncol |
optional number of grid columns |
plot_col |
output column which should be plotted |
y_title |
optional title of y-axis |
... |
any additional parameters |
All parameter combinations and exposure patterns are simulated and the mean of predictions is derived for a single study. The uncertainty is passed to the function due to computation time. Results are displayed by plotting the time series including the uncertainty interval. Observation data can be optionally displayed. Data should be provided in long format. Function plots the time (column 1) and the predictions (column 2, can be changed by the user plot_col)
a ggplot2 plot object
set.seed(124) exposure <- data.frame( time = 0:21, conc = rnorm(n = 22, mean = 0.1, sd = 0.06), trial = "T1" ) forcings <- list(temp = 12, rad = 15000) param <- list(EC50 = 0.3, b = 4.16, P_up = 0.0054) inits <- list(BM = 0.0012, E = 1, M_int = 0) scenario <- Lemna_Schmitt() %>% set_forcings(forcings) %>% set_param(param) %>% set_init(inits) sim_result <- simulate_batch( model_base = scenario, treatments = exposure, param_sample = NULL ) plot_sd( model_base = scenario, treatments = exposure, rs_mean = sim_result )
set.seed(124) exposure <- data.frame( time = 0:21, conc = rnorm(n = 22, mean = 0.1, sd = 0.06), trial = "T1" ) forcings <- list(temp = 12, rad = 15000) param <- list(EC50 = 0.3, b = 4.16, P_up = 0.0054) inits <- list(BM = 0.0012, E = 1, M_int = 0) scenario <- Lemna_Schmitt() %>% set_forcings(forcings) %>% set_param(param) %>% set_init(inits) sim_result <- simulate_batch( model_base = scenario, treatments = exposure, param_sample = NULL ) plot_sd( model_base = scenario, treatments = exposure, rs_mean = sim_result )
The method pulls available metadata from scenario objects and returns a
table with additional columns. If the argument already was a data.frame
object, the columns are appended. May overwrite existing columns of the
same name.
pull_metadata(x, model = TRUE, exposure = TRUE)
pull_metadata(x, model = TRUE, exposure = TRUE)
x |
vector of scenarios or a |
model |
|
exposure |
|
a data.frame
metsulfuron %>% pull_metadata()
metsulfuron %>% pull_metadata()
Data are from Weber 2012 publication.
Rsubcapitata
Rsubcapitata
An object of class AlgaeTKTD
of length 1.
Weber D, Schaeffer D, Dorgerloh M, Bruns E, Goerlitz G, Hammel K, Preuss TG and Ratte HT, 2012. Combination of a higher-tier flow-through system and population modeling to assess the effects of time-variable exposure of isoproturon on the green algae Desmodesmus subspictatus and Pseudokirchneriella subcapitata. Environmental Toxicology and Chemistry, 31, 899-908. doi:10.1002/etc.1765
The EffectScenario
class is the base for all of the basic scenario
types and models. It contains slots for data and settings that are
required by most models such as a vector of model parameters and a vector
of initial states. For each particular model, the class's slots are
filled with certain default or fixed values. Some models derive from this
class and add slots to store additional data.
Certain behaviors that are required to model complex processes cannot be
represented by a single EffectScenario. As an example, the parameters of a scenario
are generally fixed during the simulated time period. In order to represent a
change in parameter values, the original scenario would
need to split into two scenarios A and B which differ by parameter values
and simulated time period. By combining these scenarios to a scenario sequence,
the sequence would be treated as a single, complex scenario. See
sequence()
for more information.
Most parameters are represented by numerical types but other types are possible depending on model. Please refer to the model description which parameters are required and in which unit. Some or all parameters may be required to start a simulation. If required parameters are missing, simulation will fail with an error message.
The initial state represents the starting values of state variables when
starting a simulation. A scenario's default initial state may be insufficient
to get sensible results. It is advisable to set an initial state explicitly
when creating a new scenario, see set_init()
.
In theory, a scenario's state variables can be renamed by modifying the names
of the initial state vector. However, this is strongly
discouraged as this will affect other routines such as effect()
and epx()
and may render results useless.
Exposure refers to the concentration of toxicant an organism is exposed to. In case of aquatic organisms, this would commonly be the concentration of a toxicant in water. Other interpretations are possible depending on model assumptions.
Exposure time-series are generally represented by a data.frame
containing two
columns. The first column representing time, the second representing the
exposure level. The ordering of columns is mandatory. The column names
are essentially irrelevant but sensible names may help documenting the
scenario and its data. The rows must be ordered chronologically. A time-series
can consist of only a single row; in this case it will represent constant
exposure. Exposure time-series are set to a scenario using set_exposure()
.
Handling time-series is a costly task for the ODE solver due to consistency
checks and interpolation between time steps. How the solver interpolates
the time-series can be controlled by certain arguments to functions
such as simulate()
and effect()
. Please refer to simulate()
for a brief
overview and deSolve::forcings for a detailed description.
Exposure time-series should be kept as short as possible and as complex as needed for optimal computational efficiency.
Forcings generally refer to model parameters that change over time as part of an external function such as environmental temperature and exposure levels. Due to the importance of exposure in regulatory assessments, this R package explicitly distinguishes between environmental forcings and exposure. However, the same restrictions and features apply to both of them.
Forcing time-series are handled the same way as exposure time-series, i.e.
they are represented by a data.frame
containing two columns. The first column
representing time, the second representing the parameter that is a function of
time. The ordering of columns is mandatory. The rows must be ordered
chronologically. Forcings time-series are set using set_forcings()
.
Please refer to the Exposure section for more information on how time-series
are handled.
A scenario's simulated time period is defined by its minimum and maximum output
time. Simulation results will only be returned for the defined output times
even though the ODE solver may use smaller time steps between output times.
Output times can be explicitly set using set_times()
. The number and
distance of output times may have influence on the precision of simulation
results and numerical stability, cf. simulate()
.
Be aware that set_exposure()
will overwrite previously defined output times
if not requested otherwise.
Generally, all state variables can be used as effect endpoints but models may
provide additional endpoints. Use set_endpoints()
to enable or disable
endpoints for a scenario.
Some scenarios or models require control runs to calculate effects under
exposure. Generally, control simulations will run automatically
where needed. However, when conducting a large number of repeated simulations, e.g.
when calculating effect profiles (EPx values) or simulating moving exposure
windows, it may be computational efficient to run control simulations only
once and cache their results within the scenario. Please refer to cache_controls()
for details.
The time frame relevant for effects may be much shorter than the assessed
exposure time-series for certain organisms. This fact can be represented
by moving exposure windows which divide a long time period in a number of
consecutive windows of the same length. Each window is simulated individually and effects are
calculated. By default, methods such as effect()
will only return the
maximum effect of all considered windows but detailed results can be presented
on demand.
To use moving exposure windows, the exposure time-series must be regular,
i.e. must have an equidistant step length in time. The length of the window
is defined as the number of time steps of the exposure time-series. As an
example, assume the time-series has daily granularity and a moving window
of seven days length is required. In this case, the moving window must have
a length of seven (7) time steps. If the exposure time-series had hourly granularity,
the same window would need to have a length of 168 (=7*24) time steps. Please
refer to set_window()
for details.
name
character
, unique model name
tag
character
, an optional identifier
param
list
of parameter key-value pairs
param.bounds
named list
of parameter boundaries
param.req
character
vector of required parameters
forcings
list
of data.frame
s representing forcing time-series
forcings.req
character
vector or required model forcings data, e.g. temperature
init
list
of initial model states
times
numeric
vector of output times, beginning and end also define
the simulated period
endpoints
character
vector of endpoints to calculate results for
exposure
data.frame
with two columns representing an exposure time-series
control
list
of named numerical vectors, contains the control values
for all relevant moving windows
control.req
logical
, if TRUE
then control values are required to
calculate effects
window.length
numeric
, maximum length of the simulated period, if
window.length
is shorter than the exposure pattern, then all possible
exposure sub-patterns are evaluated for effect calculation. This is also
referred to as a moving window approach.
window.interval
numeric
, interval determining distance between moving
windows during effect calculation. First window starts at first time point
in exposure pattern.
Other scenario-related:
Transferable
Data are from Schmitt 2013 publication.
Schmitt2013
Schmitt2013
An object of class data.frame
with 56 rows and 4 columns.
Schmitt W., Bruns E., Dollinger M., and 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
A sequence of scenarios is treated as a single scenario and each scenario
is simulated one after the other. If scenario n
in a sequence was simulated,
scenario n+1
will start off in the model state where n
had ended.
Scenario sequences can be used to e.g. implement changes in model parameters
over time.
All scenarios in a sequence must fulfill the following requirements:
All scenarios must have identical state variables
The output times of all scenarios must represent a continuous time series without gaps or overlaps
Only simulation of sequences are supported, at the moment. Effects and effect profiles (EPx values) cannot be derived, yet.
an S4 object of type ScenarioSequence
# create two scenarios that need to be simulated one after the other scen1 <- minnow_it %>% set_times(0:3) scen2 <- minnow_it %>% set_times(3:6) %>% set_param(c(kd=0)) # create a sequence and assign scenarios sq <- sequence(list(scen1, scen2)) # simulate the sequence simulate(sq)
# create two scenarios that need to be simulated one after the other scen1 <- minnow_it %>% set_times(0:3) scen2 <- minnow_it %>% set_times(3:6) %>% set_param(c(kd=0)) # create a sequence and assign scenarios sq <- sequence(list(scen1, scen2)) # simulate the sequence simulate(sq)
Modifies the boundaries of model parameters for one or more scenario or caliset objects.
set_bounds(x, bounds) ## S4 method for signature 'EffectScenario,list' set_bounds(x, bounds) ## S4 method for signature 'CalibrationSet,list' set_bounds(x, bounds) ## S4 method for signature 'list,list' set_bounds(x, bounds)
set_bounds(x, bounds) ## S4 method for signature 'EffectScenario,list' set_bounds(x, bounds) ## S4 method for signature 'CalibrationSet,list' set_bounds(x, bounds) ## S4 method for signature 'list,list' set_bounds(x, bounds)
x |
|
bounds |
named list of numerical vectors, where the first level lists the parameters by name, and the second level lists the lower and upper boundary |
scenario or caliset with modified parameter boundaries
metsulfuron %>% set_bounds(list(k_phot_max = c(0, 30), k_resp = c(0, 10)))
metsulfuron %>% set_bounds(list(k_phot_max = c(0, 30), k_resp = c(0, 10)))
Effect endpoints calculated by functions such as effect()
and epx()
can be enabled and disabled. If an endpoint is not required for an assessment,
it should be disabled for reasons of computational efficiency. Please refer
to the model description for a list of available endpoints.
set_endpoints(x, endpoints)
set_endpoints(x, endpoints)
x |
vector of |
endpoints |
|
Modified EffectScenario
objects
# Only enable reproduction (R) endpoint for americamysis scenario americamysis %>% set_endpoints("R") %>% effect() # Enable endpoints length (L) and reproduction (R) americamysis %>% set_endpoints(c("L","R")) %>% effect()
# Only enable reproduction (R) endpoint for americamysis scenario americamysis %>% set_endpoints("R") %>% effect() # Enable endpoints length (L) and reproduction (R) americamysis %>% set_endpoints(c("L","R")) %>% effect()
Exposure refers to the toxicant concentration an organism is exposed to. In case of aquatic organisms, this would commonly be the concentration of a toxicant in water. Other interpretations are possible depending on model assumptions.
set_exposure(scenarios, series, ...) ## S4 method for signature 'ANY,ANY' set_exposure(scenarios, series) ## S4 method for signature 'EffectScenario,data.frame' set_exposure(scenarios, series, ...) ## S4 method for signature 'EffectScenario,ExposureSeries' set_exposure(scenarios, series, reset_times = TRUE) ## S4 method for signature 'EffectScenario,list' set_exposure(scenarios, series, ...) ## S4 method for signature 'list,list' set_exposure(scenarios, series, ...) ## S4 method for signature 'list,ANY' set_exposure(scenarios, series, ...)
set_exposure(scenarios, series, ...) ## S4 method for signature 'ANY,ANY' set_exposure(scenarios, series) ## S4 method for signature 'EffectScenario,data.frame' set_exposure(scenarios, series, ...) ## S4 method for signature 'EffectScenario,ExposureSeries' set_exposure(scenarios, series, reset_times = TRUE) ## S4 method for signature 'EffectScenario,list' set_exposure(scenarios, series, ...) ## S4 method for signature 'list,list' set_exposure(scenarios, series, ...) ## S4 method for signature 'list,ANY' set_exposure(scenarios, series, ...)
scenarios |
|
series |
|
... |
additional arguments |
reset_times |
|
Exposure time-series are generally represented by a data.frame
containing two
columns. The first column for time, the second representing the
exposure level. The ordering of columns is mandatory. The column names
are non-relevant but sensible names may help documenting the
scenario and its data. The data.frame
's rows must be ordered chronologically.
A time-series can consist of only a single row; in this case it will represent constant
exposure.
For convenience, a time-series with zero exposure can be set using
set_noexposure()
.
Handling time-series is a costly task for the ODE solver due to consistency
checks and interpolation between time steps. How the solver interpolates
the time-series can be controlled by optional arguments to functions
such as simulate()
and effect()
. Please refer to simulate()
for a brief
overview and deSolve::forcings for a detailed description.
Exposure time-series should be kept as short as possible and as complex as needed for optimal computational efficiency.
By default, the exposure time-series' time points will also be used as output
times of the scenario. Any output times previously set by set_times()
will
be lost. If this behavior is undesired, set the function argument
reset_times=FALSE
.
The functions supports modifying multiple scenarios at once: by
calling it with lists of scenario and ExposureSeries
objects. The cartesian product of all scenarios and exposure series will
be returned, iff the parameter expand = TRUE
is set.
As an example for the expand mode, two scenarios A
and B
and one
exposure series g
will result in two scenarios Ag
and Bg
, both
using exposure series g
. Two scenarios A
and B
as well as two
exposure seres g
and h
will result in four scenarios Ag
,Ah
,Bg
,
and Bh
.
list
of EffectScenario
objects
# set a data.frame as exposure series Lemna_Schmitt() %>% set_exposure(data.frame(time=c(0, 1, 2, 3), conc=c(1, 1, 0, 0))) # set one ExposureSeries es1 <- ExposureSeries(data.frame(time=0, conc=0)) Lemna_Schmitt() %>% set_exposure(es1) # set two ExposureSeries to create two scenarios es2 <- ExposureSeries(data.frame(time=5:10, conc=1)) Lemna_Schmitt() %>% set_exposure(c(es1, es2)) # set one ExposureSeries without resetting existing output times Lemna_Schmitt() %>% set_times(0:5) %>% set_exposure(es1, reset_times=FALSE)
# set a data.frame as exposure series Lemna_Schmitt() %>% set_exposure(data.frame(time=c(0, 1, 2, 3), conc=c(1, 1, 0, 0))) # set one ExposureSeries es1 <- ExposureSeries(data.frame(time=0, conc=0)) Lemna_Schmitt() %>% set_exposure(es1) # set two ExposureSeries to create two scenarios es2 <- ExposureSeries(data.frame(time=5:10, conc=1)) Lemna_Schmitt() %>% set_exposure(c(es1, es2)) # set one ExposureSeries without resetting existing output times Lemna_Schmitt() %>% set_times(0:5) %>% set_exposure(es1, reset_times=FALSE)
Parameters which change their value over time are referred to as forcings. If and what parameters can vary over time depends on the model in question. In many cases, forcings represent time-series of environmental properties.
set_forcings(x, ...) ## S4 method for signature 'EffectScenario' set_forcings(x, ...) ## S4 method for signature 'list' set_forcings(x, ...)
set_forcings(x, ...) ## S4 method for signature 'EffectScenario' set_forcings(x, ...) ## S4 method for signature 'list' set_forcings(x, ...)
x |
(vector of) scenario objects |
... |
named argument list to set as forcings |
Forcing time-series are always represented by a
data.frame
containing two columns. The first column representing time,
the second representing the parameter that is a function of time. The
ordering of columns is mandatory. The column names are essentially irrelevant
but may help documenting the scenario and its data. The rows must be
ordered chronologically. A time-series can consist of only a single row; in
this case it will represent constant conditions.
Handling forcing time-series is a costly task for the ODE solver due to consistency
checks and interpolation between timesteps. How the solver interpolates
the forcing time-series can be controlled by certain arguments to functions
such as simulate()
and effect()
. Please refer to simulate()
for a brief
overview and deSolve::forcings for a detailed description.
Forcing time-series should be kept as short as possible and as complex as needed for optimal computational efficiency.
Modified scenarios
# constant values will be automatically converted to a data.frame Lemna_Schmitt() %>% set_forcings(temp=20) -> lemna lemna@forcings # setting multiple forcings at once df <- data.frame(t=0:14, temp=rnorm(15, mean=20)) # random temperature series Lemna_Schmitt() %>% set_forcings(temp=df, rad=15000) -> lemna lemna@forcings # forcings can also be supplied as a named list Lemna_Schmitt() %>% set_forcings(list(temp=20, rad=15000)) -> lemna lemna@forcings
# constant values will be automatically converted to a data.frame Lemna_Schmitt() %>% set_forcings(temp=20) -> lemna lemna@forcings # setting multiple forcings at once df <- data.frame(t=0:14, temp=rnorm(15, mean=20)) # random temperature series Lemna_Schmitt() %>% set_forcings(temp=df, rad=15000) -> lemna lemna@forcings # forcings can also be supplied as a named list Lemna_Schmitt() %>% set_forcings(list(temp=20, rad=15000)) -> lemna lemna@forcings
The initial state represents the starting values of a scenario's state variables when starting a simulation. A scenario's default initial state may be insufficient to get sensible results.
set_init(x, init) ## S4 method for signature 'vector' set_init(x, init) ## S4 method for signature 'EffectScenario' set_init(x, init)
set_init(x, init) ## S4 method for signature 'vector' set_init(x, init) ## S4 method for signature 'EffectScenario' set_init(x, init)
x |
vector of |
init |
named numeric vector |
In theory, a scenarios's state variables can be renamed by modifying the names
of the initial state vector. However, this is strongly
discouraged as this will affect other routines such as effect()
and epx()
and may render results useless.
modified EffectScenario
objects
# Set initial biomass to 1.0 metsulfuron %>% set_init(c(BM=1.0)) %>% simulate()
# Set initial biomass to 1.0 metsulfuron %>% set_init(c(BM=1.0)) %>% simulate()
Updates the model parameter MoA
to a certain value
set_mode_of_action(x, code) set_moa(x, code)
set_mode_of_action(x, code) set_moa(x, code)
x |
vector of scenarios |
code |
a code for a mode of action, refer to model description for details |
modified scenarios
set_moa()
: Shorthand version
# Set MoA=8, i.e. hazard during oogenesis americamysis %>% set_mode_of_action(8) %>% effect(method="ode45") # alternative approach using the parameter directly americamysis %>% set_param(c(MoA=8)) %>% effect(method="ode45")
# Set MoA=8, i.e. hazard during oogenesis americamysis %>% set_mode_of_action(8) %>% effect(method="ode45") # alternative approach using the parameter directly americamysis %>% set_param(c(MoA=8)) %>% effect(method="ode45")
The scenarios current exposure is replaced by a constant exposure time-series
of value zero(0.0
). Output times are unaffected.
set_noexposure(x)
set_noexposure(x)
x |
vector of scenarios |
vector of scenarios
# Derive effect size in sample scenario without toxicant exposure minnow_it %>% set_noexposure() %>% effect()
# Derive effect size in sample scenario without toxicant exposure minnow_it %>% set_noexposure() %>% effect()
Modifies the parameters of one or more EffectScenario
objects.
set_param(x, param) ## S4 method for signature 'EffectScenario,vector' set_param(x, param) ## S4 method for signature 'EffectScenario,ParameterSet' set_param(x, param) ## S4 method for signature 'list,ParameterSet' set_param(x, param) ## S4 method for signature 'list,vector' set_param(x, param) ## S4 method for signature 'ScenarioSequence,vector' set_param(x, param) ## S4 method for signature 'ScenarioSequence,ParameterSet' set_param(x, param)
set_param(x, param) ## S4 method for signature 'EffectScenario,vector' set_param(x, param) ## S4 method for signature 'EffectScenario,ParameterSet' set_param(x, param) ## S4 method for signature 'list,ParameterSet' set_param(x, param) ## S4 method for signature 'list,vector' set_param(x, param) ## S4 method for signature 'ScenarioSequence,vector' set_param(x, param) ## S4 method for signature 'ScenarioSequence,ParameterSet' set_param(x, param)
x |
object(s) to modify |
param |
named numeric |
Most parameters are represented by numerical types but other types are possible depending on model. Please refer to the model description which parameters are required and in which unit. Some or all parameters may be required to start a simulation. If required parameters are missing, simulation will fail with an error message.
Vector of modified objects
Lemna_Schmitt() %>% set_param(c(Emax=1,EC50=0.12))
Lemna_Schmitt() %>% set_param(c(Emax=1,EC50=0.12))
Sets the user-defined, custom tag of a scenario. Tags can be helpful to quickly distinguish scenarios of the same model type.
set_tag(x, tag)
set_tag(x, tag)
x |
(vector of) |
tag |
vector of |
(vector of) modified EffectScenario
objects
# set a custom tag myscenario <- GUTS_RED_SD() %>% set_tag("My Custom Tag") # returns `My Custom Tag` get_tag(myscenario) # the tag also appears in the scenario overview myscenario
# set a custom tag myscenario <- GUTS_RED_SD() %>% set_tag("My Custom Tag") # returns `My Custom Tag` get_tag(myscenario) # the tag also appears in the scenario overview myscenario
Minimum and maximum output times define the simulated period for a scenario.
Simulation results will be returned for each output time, see simulate()
.
set_times(x, times)
set_times(x, times)
x |
vector of scenarios |
times |
|
Be aware that output times may be modified by set_exposure()
. Precision of
simulation results may be influenced by chosen output times, see simulate()
for more information.
Vector of modified scenarios
# Set simulated period to [2,4] with output intervals of length 1 minnow_it %>% set_times(c(2,3,4)) # Decrease output interval length to 0.1 minnow_it %>% set_times(seq(2, 4, 0.1))
# Set simulated period to [2,4] with output intervals of length 1 minnow_it %>% set_times(c(2,3,4)) # Decrease output interval length to 0.1 minnow_it %>% set_times(seq(2, 4, 0.1))
A transfer refers to an event where a certain amount of biomass (BM) is moved to a new medium after a period of time. This feature replicates a procedure occurring e.g. in Lemna effect studies and may be necessary to recreate study results.
set_transfer(x, interval, times, biomass, scaled_comp) ## S4 method for signature 'ANY' set_transfer(x, interval, times, biomass, scaled_comp) ## S4 method for signature 'Transferable' set_transfer(x, interval, times, biomass, scaled_comp)
set_transfer(x, interval, times, biomass, scaled_comp) ## S4 method for signature 'ANY' set_transfer(x, interval, times, biomass, scaled_comp) ## S4 method for signature 'Transferable' set_transfer(x, interval, times, biomass, scaled_comp)
x |
vector of |
interval |
optional |
times |
optional |
biomass |
optional |
scaled_comp |
optional |
Any transfer time point must also be an output time point. If a transfer occurs, simulation results of that time point will report the model state before the transfer. Be aware that in order to use transfers at regular intervals, the simulation must start at time point zero.
At each transfer, a defined amount of biomass is transferred to a new medium.
This is modeled by interrupting the simulation at a transfer time point, modifying
the biomass level BM
, and scaling affected compartments according to new
biomass levels. Scaling of compartments depending on biomass, such as
internal toxicant mass, is necessary to correctly reflect mass balances and
concentrations over time.
Transferred biomass is set using the biomass
parameter. Is is either a
single numerical value in which case the same biomass level is set at each
transfer. Or it is a vector of numerical values with the same length as the
times
parameter in which case a custom biomass level can be set for each
transfer. Multiple biomass levels can only be set in conjunction with
custom transfer time points.
Some scenario types define default values for transferred biomass based on common study set ups.
Transfers can occur either in regular intervals of time or at selected, custom
time points. For regular intervals, the parameter interval
is set to a single
numeric value which has the same unit as the scenario's time dimension. As an
example: if a scenario uses the unit of days for time, the transfer interval
is also specified in days:
Transfers occurring at custom time points are set by passing a numerical vector
to the parameter times
. The time points' units must match with the unit of
time in the scenario. A custom transfer time point must not occur at the
starting time point of a simulation.
Some compartments depend on biomass to correctly reflect mass balances and
concentrations over time, such as internal toxicant mass. These compartments
need to be scaled linearly to reflect the change in biomass levels.
The parameter scaled_comp
accepts a character vector of compartment names
which are scaled at each transfer. This parameter should only be used with
custom, user-defined models. If no compartment needs to be scaled, set or
use the default value of character(0)
.
Modified EffectScenario
objects
# Simulate biomass transfer of 50 *g dw/m²* at a regular interval of 7 *days* metsulfuron %>% set_transfer(interval=7, biomass=50) %>% simulate() # Simulate irregular biomass transfers occuring at days 5, 10, and 12 metsulfuron %>% set_transfer(times=c(5,10,12), biomass=50) %>% simulate() # Simulate irregular transfers with changing amounts of transferred biomass metsulfuron %>% set_transfer(times=c(5,10,12), biomass=c(50,20,10)) %>% simulate()
# Simulate biomass transfer of 50 *g dw/m²* at a regular interval of 7 *days* metsulfuron %>% set_transfer(interval=7, biomass=50) %>% simulate() # Simulate irregular biomass transfers occuring at days 5, 10, and 12 metsulfuron %>% set_transfer(times=c(5,10,12), biomass=50) %>% simulate() # Simulate irregular transfers with changing amounts of transferred biomass metsulfuron %>% set_transfer(times=c(5,10,12), biomass=c(50,20,10)) %>% simulate()
Exposure windows are defined as a period of time at the scale of the exposure series.
As an example: if an exposure series has an hourly time step, a window length
of 24
will consider the exposure within 24 hours intervals for effect
calculation. The same applies for the window interval, i.e. the period between
considered exposure windows. Set length=-1
to disable moving windows.
set_window(x, length, interval)
set_window(x, length, interval)
x |
vector of |
length |
|
interval |
|
modified EffectScenario
objects
# calculate the maximum effect for all windows of 10 days length metsulfuron %>% set_window(length=10, interval=1) %>% effect()
# calculate the maximum effect for all windows of 10 days length metsulfuron %>% set_window(length=10, interval=1) %>% effect()
The supplied EffectScenario
is passed on to the ODE solver for numerical
integration. Internally, simulate()
is split up into several functions
dedicated to particular models, e.g. one for GUTS and one for Lemna type
models. The package will take care of using the correct function
for each model when simulate()
is called.
simulate(x, ...) ## S4 method for signature 'EffectScenario' simulate(x, ...) ## S4 method for signature 'Transferable' simulate(x, ...) ## S4 method for signature 'ScenarioSequence' simulate(x, ...)
simulate(x, ...) ## S4 method for signature 'EffectScenario' simulate(x, ...) ## S4 method for signature 'Transferable' simulate(x, ...) ## S4 method for signature 'ScenarioSequence' simulate(x, ...)
x |
scenario to simulate |
... |
additional parameters passed on to ODE solver |
Simulation results are returned as a time-series for each state variable. Some models provide additional properties describing the model state, e.g. the internal concentration of a toxicant within the organism. Refer to the respective scenario for more information.
Additional arguments to simulate()
will be passed on to deSolve::ode()
which enables control of the numerical integration parameters.
The minimum and maximum of given time points generally define the simulated
period. Function argument times
overrides settings of the scenario, i.e.
time points set in scenario@times
. However, the simulation can be limited to
a subset of time points by enabling a moving exposure window, see set_window()
.
Results will be returned for each time point. Precision of the numeric solver
may be affected by chosen output times in certain cases. Hence, small deviations
in results should be expected if different output times are set. This effect
can be mitigated by either defining are sufficiently small time step for the solver
using argument hmax
or by decreasing the error tolerances atol
and rtol
.
These arguments are passed to the solver, see e.g. deSolve::lsoda()
for details.
Some models support adding intermediary model variables to the return value
of simulate()
. Analyzing the additional outputs may be helpful to understand
model behavior and support finding potential issues in model parameterization.
Optional outputs are enabled by setting the parameter nout
to a value greater
than zero. If nout
is set to n
, then the first n
optional output columns
will be returned along the normal simulation result.
Which optional outputs are available depends on the model/scenario at hand.
Please refer to the model documentation for details.
As an example, the GUTS-RED-IT model supports adding the
external toxicant concentration to the output by setting nout=1
:
minnow_it %>% simulate(nout=1)
Each model was assigned a default ODE solver which handles most of the
occurring inputs well. In most cases, this will be an explicit numerical
scheme of the Runge-Kutta family with variable step width. For certain extreme
parameters settings, such as very high uptake/permeability of the contaminant
or exposure series which represent step functions, the numerical approximation
might deteriorate and yield invalid results. In this case try to decrease the
allowed max step width by setting the argument hmax
with various values.
Start with hmax=1
and decrease the value by orders of 10. It is not
possible or desirable to reduce hmax
to extremely small values, as the
ODE solver will require more CPU time and simulation will become inefficient.
Oftentimes, it will be computational more efficient to adapt the solver's
error tolerances atol
and rtol
than reducing the step width hmax
to achieve
stable numerics. Start by decreasing deSolve's default values by
orders of ten until the simulation yields acceptable results, see e.g.
deSolve::lsoda()
for more information on error tolerances.
As an alternative to adapting solver parameters, it might be worthwhile to
try other numerical schemes which can handle stiff ODEs, such as Radau, LSODA,
or LSODES. To change solvers, set the method
argument.
To select e.g. the Radau scheme, set method="radau"
. For LSODA, set method="lsoda"
.
Radau performs better than LSODA in some cases, as the latter method can return
biologically nonsensical results without raising an error. See deSolve::ode()
for details on available ODE solvers.
A data.frame
with the time-series of simulation results
minnow_sd %>% simulate() # tidy syntax simulate(minnow_sd) # base R syntax # Set new output times minnow_sd %>% simulate(times=c(0,4)) # Modify simulated time frame minnow_sd %>% simulate(times=c(0,10)) # Use an alternative exposure profile than defined in the scenario minnow_sd %>% set_exposure(data.frame(t=0,c=10), reset_times=FALSE) %>% simulate() ## ## Precision of results # A large number of output times forces smaller solver time steps minnow_it %>% simulate(times=seq(0,1,0.001)) %>% tail() # Defining only two output times allows for larger solver time steps minnow_it %>% simulate(times=c(0,1)) # The difference between results can be reduced by limiting the solver's # maximum step length minnow_it %>% simulate(times=c(0,1), hmax=0.001) # The same numerical precision can be achieved by switching to # the Radau scheme minnow_it %>% simulate(times=c(0,1), method="radau") # A very small step length may reduce the difference even further but may # also exceed the allowed number of steps per output interval. The following # simulation will be aborted with a solver error: try( minnow_it %>% simulate(times=c(0,1), hmax=0.0001) ) # However, the solver's max number of allowed steps can be increased: minnow_it %>% simulate(times=c(0,1), hmax=0.0001, maxsteps=10^5)
minnow_sd %>% simulate() # tidy syntax simulate(minnow_sd) # base R syntax # Set new output times minnow_sd %>% simulate(times=c(0,4)) # Modify simulated time frame minnow_sd %>% simulate(times=c(0,10)) # Use an alternative exposure profile than defined in the scenario minnow_sd %>% set_exposure(data.frame(t=0,c=10), reset_times=FALSE) %>% simulate() ## ## Precision of results # A large number of output times forces smaller solver time steps minnow_it %>% simulate(times=seq(0,1,0.001)) %>% tail() # Defining only two output times allows for larger solver time steps minnow_it %>% simulate(times=c(0,1)) # The difference between results can be reduced by limiting the solver's # maximum step length minnow_it %>% simulate(times=c(0,1), hmax=0.001) # The same numerical precision can be achieved by switching to # the Radau scheme minnow_it %>% simulate(times=c(0,1), method="radau") # A very small step length may reduce the difference even further but may # also exceed the allowed number of steps per output interval. The following # simulation will be aborted with a solver error: try( minnow_it %>% simulate(times=c(0,1), hmax=0.0001) ) # However, the solver's max number of allowed steps can be increased: minnow_it %>% simulate(times=c(0,1), hmax=0.0001, maxsteps=10^5)
A convenience function to simulate a single base scenario with one or more exposure series. This aims at reproducing the setup and results of common effect studies.
simulate_batch(model_base, treatments, param_sample = deprecated())
simulate_batch(model_base, treatments, param_sample = deprecated())
model_base |
effect scenario object with mean parameters |
treatments |
treatments exposure levels as data frame (time, conc, trial) |
param_sample |
deprecated parameter, no longer in use |
A scenario contains only one exposure series. However, laboratory experiments commonly examine the effects of multiple exposure levels on a biological system. A batch simulation approach would involve running multiple simulations with varying exposure or treatment conditions. To illustrate, if the objective is to examine the impact of a substance on cell growth, the simulation model could be designed to replicate the cell growth dynamics under varying concentrations of the substance. Each simulation run would represent a specific exposure level, ranging from low to high concentrations of the chemical. To simulate such a laboratory experiment, the simulate_batch function can be used. All exposure series are saved in the treatment argument. The first column contains the time, the second column the concentration, and the third column the trial name (exposure level, e.g. 'T1', 'T2', 'T3').
a data.frame
with simulation results
t1 <- data.frame(time=0:10, conc=0, trial="control") # 1st treatment level t2 <- data.frame(time=0:10, conc=1, trial="T1") # 2nd treatment level treatments <- rbind(t1, t2) metsulfuron %>% simulate_batch(treatments)
t1 <- data.frame(time=0:10, conc=0, trial="control") # 1st treatment level t2 <- data.frame(time=0:10, conc=1, trial="T1") # 2nd treatment level treatments <- rbind(t1, t2) metsulfuron %>% simulate_batch(treatments)
Please refer to the Modeling Howto vignette on how to implement custom
models by overloading the solver
function.
solver(scenario, times, ...) ## S4 method for signature 'ANY' solver(scenario, times, ...) ## S4 method for signature 'AlgaeWeber' solver( scenario, times, approx = c("linear", "constant"), f = 1, method = "lsoda", hmax = 0.1, ... ) ## S4 method for signature 'AlgaeTKTD' solver( scenario, times, approx = c("linear", "constant"), f = 1, method = "lsoda", hmax = 0.1, ... ) ## S4 method for signature 'AlgaeSimple' solver( scenario, times, approx = c("linear", "constant"), f = 1, method = "ode45", hmax = 0.01, ... ) ## S4 method for signature 'DebAbj' solver(scenario, times, ...) ## S4 method for signature 'DebTox' solver( scenario, times, approx = c("linear", "constant"), f = 1, rule = 2, method = "ode45", ... ) ## S4 method for signature 'DebDaphnia' solver(scenario, times, ...) ## S4 method for signature 'GutsRedSd' solver(scenario, times, ...) ## S4 method for signature 'GutsRedIt' solver(scenario, times, ...) ## S4 method for signature 'LemnaSetac' solver(scenario, times, ...) ## S4 method for signature 'LemnaSchmitt' solver(scenario, times, ...) ## S4 method for signature 'MyrioExp' solver(scenario, times, ...) ## S4 method for signature 'MyrioLog' solver(scenario, times, ...)
solver(scenario, times, ...) ## S4 method for signature 'ANY' solver(scenario, times, ...) ## S4 method for signature 'AlgaeWeber' solver( scenario, times, approx = c("linear", "constant"), f = 1, method = "lsoda", hmax = 0.1, ... ) ## S4 method for signature 'AlgaeTKTD' solver( scenario, times, approx = c("linear", "constant"), f = 1, method = "lsoda", hmax = 0.1, ... ) ## S4 method for signature 'AlgaeSimple' solver( scenario, times, approx = c("linear", "constant"), f = 1, method = "ode45", hmax = 0.01, ... ) ## S4 method for signature 'DebAbj' solver(scenario, times, ...) ## S4 method for signature 'DebTox' solver( scenario, times, approx = c("linear", "constant"), f = 1, rule = 2, method = "ode45", ... ) ## S4 method for signature 'DebDaphnia' solver(scenario, times, ...) ## S4 method for signature 'GutsRedSd' solver(scenario, times, ...) ## S4 method for signature 'GutsRedIt' solver(scenario, times, ...) ## S4 method for signature 'LemnaSetac' solver(scenario, times, ...) ## S4 method for signature 'LemnaSchmitt' solver(scenario, times, ...) ## S4 method for signature 'MyrioExp' solver(scenario, times, ...) ## S4 method for signature 'MyrioLog' solver(scenario, times, ...)
scenario |
scenario object |
times |
numeric vector of output times, overrides any scenario setting |
... |
additional parameters passed on to |
approx |
string, interpolation method of exposure series, see deSolve::forcings |
f |
if |
method |
string, numerical solver used by |
hmax |
numeric, maximum step length in time, see |
rule |
if |
Please note that not all solvers support setting the parameters listed here. In addition, some solvers may set reasonable default values for e.g. maximum step length in time, but not all do. Please check the model documentation for details.
data.frame
with simulation results
solver(ANY)
: Default solver, raises an error
solver(AlgaeWeber)
: numerically integrates Algae_Weber models
solver(AlgaeTKTD)
: numerically integrates Algae_TKTD models
solver(AlgaeSimple)
: numerically integrates Algae_Simple models
solver(DebAbj)
: Numerically integrates DEB_abj models
solver(DebTox)
: Numerically integrates DEBtox scenarios
solver(DebDaphnia)
: (deprecated) Numerically integrates DEBtox_Daphnia scenarios
solver(GutsRedSd)
: Numerically integrates GUTS-RED-SD models
solver(GutsRedIt)
: Numerically integrates GUTS-RED-IT models
solver(LemnaSetac)
: Numerically integrates Lemna_SETAC models
solver(LemnaSchmitt)
: Numerically integrates Lemna_Schmitt models
solver(MyrioExp)
: Numerically integrates MyrioExp
models
solver(MyrioLog)
: Numerically integrates MyrioLog
models
Derives the survival rate of individuals for
Reduced GUTS models. Function was replaced by output of simulate()
and
will be removed in a later version.
survival(scenario, ...)
survival(scenario, ...)
scenario |
an |
... |
additional parameters passed on to |
The survival rate describes the survival probability at each
time point. The function simulates the GUTS scenario and appends a column
survival
to the simulation result. A value of one (1.0
) denotes that
all individuals survive. A value of zero (0.0
) denotes that no individuals
survived.
Only available for Reduced GUTS models, see GUTS-RED-models. The equations were described by EFSA (2018).
a data.frame
containing simulation results
EFSA PPR Panel (EFSA Panel on Plant Protection Products and their Residues), Ockleford C, Adriaanse P, Berny P, et al., 2018: Scientific Opinion on the state of the art of Toxicokinetic/Toxicodynamic (TKTD) effect models for regulatory risk assessment of pesticides for aquatic organisms. EFSA Journal 2018; 16(8):5377, 188 pp. doi:10.2903/j.efsa.2018.5377
# calculate survival rate minnow_it %>% survival() # plot survival over time based on a random exposure profile minnow_sd %>% set_exposure(data.frame(t=1:100, c=runif(100)*10)) %>% survival() -> df plot(df$time, df$survival, "l")
# calculate survival rate minnow_it %>% survival() # plot survival over time based on a random exposure profile minnow_sd %>% set_exposure(data.frame(t=1:100, c=runif(100)*10)) %>% survival() -> df plot(df$time, df$survival, "l")
By inheriting from class Transferable
, a scenario's behavior can be
extended to support transfer and reset of biomass at dedicated points during
simulation.
transfer.times
numeric
, vector of time points at which transfers occur,
e.g. c(7,10,14)
transfer.interval
numeric
, interval length until frond transfer to new
medium
transfer.biomass
numeric
, amount of biomass transferred to new medium
transfer.comp.biomass
character
state variable which describes
biomass
transfer.comp.scaled
character
vector of state variable which will
be scaled 1:1 when biomass is modified, e.g. internal toxicant mass
Models supporting biomass transfer can be instructed to move a fixed amount of biomass to a new medium after a period of time. This feature replicates a procedure occurring in e.g. Lemna effect studies and may be necessary to recreate study results.
The biomass transfer feature assumes that always a fixed amount of
biomass is transferred. Transfers can occur at any fixed point in time or
in regular intervals. During a transfer, the biomass is reset to the
transferred amount and additional compartments can be scaled 1:1 accordingly,
to e.g. reflect the change in internal toxicant mass when biomass is modified.
Transfer settings can be modified using set_transfer()
.
Any transfer time point must also be an output time point. If a transfer occurs, simulation results of that time point will report the model state before the transfer. Be aware that in order to use transfers at regular intervals, the simulation must start at time point zero.
Other scenario-related:
Scenarios