simulate: Simulate an effect scenario

In cvasi: Calibration, Validation, and Simulation of TKTD Models

Simulate an effect scenario

Description

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.

Usage

simulate(x, ...)

## S4 method for signature 'EffectScenario'
simulate(x, ...)

## S4 method for signature 'Transferable'
simulate(x, ...)

## S4 method for signature 'ScenarioSequence'
simulate(x, ...)

Arguments

x	scenario to simulate
...	additional parameters passed on to ODE solver

Details

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.

Output times and windows

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.

Optional output variables

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)

Numerical precision and stability

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.

Value

A data.frame with the time-series of simulation results

Examples

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)

cvasi documentation built on Sept. 23, 2024, 9:08 a.m.