Learn R Programming

cvasi (version 1.4.0)

simulate: 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 EffectScenario simulate(x, ...)

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

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

# S4 method for SimulationBatch simulate(x, ...)

Value

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

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 define the simulated period. However, the simulation can also be limited to a subset of time points by enabling a moving exposure window, see set_window().

Results will be returned for each output 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.

Examples

Run this code
# base R syntax
simulate(minnow_sd)
# tidy syntax with the same result
minnow_sd %>% simulate()

# Extend the simulated time frame to the interval [0, 10]
minnow_sd %>%
  set_times(seq(0, 10)) %>%
  simulate()

# Use an alternative exposure profile, but keep the original output times
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 %>%
  set_times(seq(0, 10, 0.001)) %>%
  simulate() %>%
  tail()

# Defining only two output times allows the ODE solver to make larger steps
# in time during numerical integration. However, results can become
# imprecise.
minnow_long <- minnow_it %>% set_times(c(0, 10))
minnow_long %>% simulate()

# Numerical precision of results can be increased by limiting the solver's
# maximum step length in time using argument `hmax`.
minnow_long %>% simulate(hmax=0.005)

# A similar numerical precision can be achieved by switching to an alternative
# numerical integration scheme, such as the Radau scheme, without limiting
# the step length.
minnow_long %>% simulate(method="radau")

# Reducing the step length even further may increase numerical precision, but
# may exceed the solver's allowed number of integration steps per output interval.
# The following simulation will be aborted with a solver error:
try(
  minnow_long %>% simulate(hmax=0.001)
)

# However, the solver's maximum number of allowed steps can be increased,
# if needed, using the argument `maxsteps`:
minnow_long %>% simulate(hmax=0.001, maxsteps=10^5)

Run the code above in your browser using DataLab