EcotoxSystems.jl
Quickstart
To get acquainted with using EcotoxSystems.jl, you can use the default model and default parameters. <br>
Let's start with simulating individual life-history based on the default parameters:
p = deepcopy(EcotoxSystems.debkiss_defaultparams)
p.glb.t_max = 21.
sim = EcotoxSystems.ODE_simulator(p)
first(sim, 5)Here, we first create a copy of the default parameters. These contain two components:
glbfor global parameters. This includes simulation settings like the maximum simulated time-span, and forcings such as the food input ratedX_in.spcfor species-specific parameters. For the default model, these are DEB-TKTD and some auxiliary parameters.
The second line,
p.glb.t_max = 21.adjusts the simulated time-span. <br> With
sim = EcotoxSystems.ODE_simulator(p), we then simulate the model. <br> ODE_simulator is the central function to execute models which are purely ODE-based. <br> The output `sim is a DataFrame containing all state variables over time. <br> Storing the solution in a DataFrame can be convenient, but, for simple models, consumes a considerable proportion of computation time. It is possible to skip the conversion and retrieve the ODEsolution instead by setting returntype=EcoxSystems.odesol. <br>
Instead of simulating the life-history of a single individual, we can also simulate population dynamics. To do so for the default model, we need to adjust the global parameters, set plausible value for the starvation parameters (which were irrelevant before), and call IBM_simulator:
using Distributions, ProgressMeter
p = deepcopy(EcotoxSystems.debkiss_defaultparams)
# adjusting global parameters
p.glb.t_max = 365. # simulated timespan [d]
p.glb.V_patch = 0.5 # simulated volume [L]
p.glb.dX_in = 10_000 # provide more food
# adjusting species-level parameters
p.spc.Z = Truncated(Normal(1, 0.1), 0, Inf) # induce individual variability
p.spc.W_S_rel_crit = 0.66
p.spc.h_S = -log(0.5)
sims = @replicates EcotoxSystems.IBM_simulator(p, showinfo = Inf) 10 # run replicated simulations
To deviate from the default model, there are a number of components to tweak in order to fully specify the system:
init_global_statevars: A function that initializes global state variables.global_ode!: ODE-portion of the dynamics of global states.global_rules!: Rule-based portion of the dynamics of global states.init_individual_statevars: A function that initializes individual-level state variables.individual_ode!: ODE-portion of the dynamics of individual-level states.individual_rules!: Rule-based portion of the dynamics of individual-level states.gen_ind_params: A function that translates species-specicif parametersspcinto individual-level parametersind. The default isgenerate_individual_paramsand is generic, apart from assuming thatspccontains an entry forZandpropagate_zoom.
API
EcotoxSystems.Euler! — Method
Euler!(u::CVOrParamStruct, du::CVOrParamStruct, dt::Real)::NothingApply Euler scheme to state variables.
EcotoxSystems.IBM_simulator — Method
IBM_simulator(
p::CVOrParamStruct;
global_ode! = default_global_ODE!,
global_rules! = default_global_rules!,
init_global_statevars = initialize_global_statevars,
individual_ode! = default_individual_ODE!,
individual_rules! = default_individual_rules!,
init_individual_statevars = initialize_individual_statevars,
dt = 1/24,
saveat = 1,
record_individuals = true,
showinfo::Number = Inf
)Simulate the individual-based version of the default model.
using EcotoxSystems
p = DEB.params()
sim = DEB.IBMsimulator(p)For explanation of arguments, see IndividualBasedModel.
EcotoxSystems.ODE_simulator — Method
ODE_simulator(
p::CVOrParamStruct;
alg = Tsit5(),
saveat = 1,
reltol = 1e-6,
model = default_ODE!,
statevars_init = initialize_statevars,
gen_ind_params = generate_individual_params,
param_links::Union{Nothing,NamedTuple} = nothing,
callbacks = DEBODE_callbacks,
returntype::ReturnType = dataframe,
kwargs...
)Run the model as ODE system. This function is essentially a wrapper around OrdinaryDiffEq.solve with additional input and output processing.
args:
p: Parameters, given as component vector.
At least two components are given: glb for global parameters, spc for species-level parameters. glb as to contain an entry t_max for the maximum simulation time. spc has to contain entries Z and propagate_zoom for the zoom factor and the parameters which are affected by the zoom factor.
kwargs:
The following kwargs are used internally by OrdinaryDiffEq.solve. See OrdinaryDiffEq documentation for more information.
alg: ODE solving algorithm to use..saveat: Interval or time-points at which to save the solution. Default is 1.reltol: Relative tolerance of ODE solution, default is 1e-6.model: ODE system to solve.callbacks: A callback set, i.e. pairs of conditions and events. Used to handle discontinuities.
In addition we have some kwargs that are used to further process inputs and outputs:
statevars_init: Function with signaturestatevars_init(p)that defines initial state variables as component vector. Components typically match those in the parameter vector.- `genindparams: Function that converts species-level parameters to individual-level parameters, for example by replacing parameters which are given as distributions with a random sample from the distribution. The inputs and outputs of this function should contain all components.
returntype: Indicating of how to return the result. Currently allowed aredataframe(complete solution converted to aDataFrame) andodesol(the ODE solution object as returned byOrdinaryDiffEq). Default isdataframe.
Run a model as purely ODE-based system:
using EcotoxSystems
sim = DEB.ODE_simulator(EcotoxSystems.debkiss_defaultparams(p))EcotoxSystems.calc_SL_max — Method
calc_SL_max(spc::CVOrParamStruct)::Float64Calculate maximum structural length slmax [m^(1/3)]
EcotoxSystems.calc_S_max — Method
calc_S_max(spc::CVOrParamStruct)::Float64Calculate maximum structural mass smax [m]
EcotoxSystems.debkiss! — Method
DEBkiss model without tox component.
EcotoxSystems.debkiss_global_statevars — Method
initialize_global_statevars(p::CVOrParamStruct)Function to initialize global state variables. In the default model, these are the resource abundance X, external chemical stressor concentration C_W and population size N.
Global state variables can be extended, modified or replaced in the same way as individual-level state variables.
EcotoxSystems.debkiss_individual_statevars — Method
debkissindividualstatevars( p::CVOrParamStruct; id = 1., cohort = 0.)::CVOrParamStruct
This function defines the individual-level state variables and their initial values for the default debkiss model.
p is a parameter vector containing all components which are relevant for this simulation (including global state varaiables). id is an individual's unique identifier in the IBM simulation.
cohortis the index of the cohort an indvidiual belongs to in the IBM simulation (i.e. initial generation is cohort 0, their offsrping is cohort 1, etc.)
To add state variables or overwrite the default initial value, one can define a wrapper around this function:
using ComponentArrays
custom_ind_statevars(p; kwargs...) = ComponentVector(
EcotoxSystems.initialize_individual_statevars();
new_statevar = 0.
)For IBM simulations, the new function custom_ind_statevars can be supplied as keyword-argument init_individual_statevars to IBM_simulator:
sim = IBM_simulator(p; init_individual_statevars = custom_ind_statevars)For the ODE, simulator, the state variables for all components are initialized simultaneously, and this function does not take any additional keyword-arguments:
custom_statevars(p) = ComponentVector( # state variables for some completely new model
glb = initialize_global_statevars(p),
spc = custom_ind_statevars(p)
)Alternatively, the custom ComponentVector can of course be defined from scratch in the new initialization functions for state variables:
custom_statevars(p) = ComponentVector( # state variables for some completely new model
glb = ComponentVector(t_max = 10., C_W = 0.),
spc = ComponentVector(N = 1.)
)EcotoxSystems.debkiss_mixture_IA! — Method
DEBkiss model with mixture toxicity, assuming independent action (IA) across all components. Supports an arbitrary number of stressors and stressor/PMoA combinations.
EcotoxSystems.default_global_rules! — Method
default_global_rules!(m)Global rule-based portion of the default model.
EcotoxSystems.default_individual_rules! — Method
default_individual_rules(a::AbstractIndividual, m::AbstractIBM)::NothingDefines the default rule-based portion for DEBIndividuals. <br>
The event functions which are used as callbacks during ODE solving are here re-used to apply rules for life stage transitions.
A crude rule for starvation mortality is implemented, applying a constant hazard rate of a certain relative amount of structural mass is lost.
Reproduction is assumed to occur in fixed time intervals, according to spc.tau_R.
EcotoxSystems.exposure — Method
exposure(
simulator::Function,
p::CVOrParamStruct,
Cmat::Union{Vector{R},Matrix{R}}
) where R <: RealSimulate constant chemical exposure with over an arbitrary number of treatments and chemical stressors, defined in Cmat.
The columns of C_W are stressors, the rows are treatments. <br> If C_W is supplied as a Vector, it is assumed that the elements represent different levels of single-stressor exposure. <br> That means, to simulate a single-stressor experiment with three treatments, do
Cmat = [0., 1., 2,], internally creating a 1 x 3 matrix with exposure concentrations 0, 1 and 2.
A single treatment with multiple stressors can be defined as
Cmat = [0 1 2;]Here, we would have three stressors with the simultaneous exposure concentrations 0,1,2. <br>
Defining four treatments for two stressors looks like this:
Cmat = [
0.0 0.0;
0.0 0.5;
1.0 0.0;
0.5 1.0
]This exposure matrix corresponds to a ray design with constant exposure ratios.
EcotoxSystems.filter_individuals! — Method
filter_individuals!(m::AbstractIBM)Remove individuals which have been flagged to die after the current time-step. Individuals for which the condition u.ind.cause_of_death == 0 applies are retained.
EcotoxSystems.generate_individual_params — Method
generate_individual_params(p::CVOrParamStruct; kwargs...)
Generate individual-specific parameter set from species-specific parameter set.If a parameter entry is a distribution, a random sample is taken.
This also works for Vectors of distributions. The kwargs need to be supplemented with additional components if there are more than just a global and an individual-level component.
EcotoxSystems.get_global_statevars! — Method
get_global_statevars!(a::AbstractIndividual, m::AbstractIBM)::NothingRetrieve global state variables and derivatives for use in the individual-level model step.
EcotoxSystems.groupedlineplot! — Method
groupedlineplot(x, y, g, q; estimator)Version of lineplot with additional grouping variable g.
EcotoxSystems.groupedlineplot! — Method
groupedlineplot(x, y, g, q; estimator)Version of lineplot with additional grouping variable g.
EcotoxSystems.groupedlineplot — Method
groupedlineplot(x, y, g, q; estimator)Version of lineplot with additional grouping variable g.
EcotoxSystems.individual_step! — Method
individual_step!(a::Agent, m::Model)::NothingThe individual-level model step follows a generic pattern:
First the ODE-portion of the model is executed, and the corresponding state variables are updated using the Euler scheme.
Then the rule-based portion of the model is executed. These are all the functions which cannot / should not be expressed as part of an ODE. At the minimum, this will include life stage transitions, reproduction and death of individuals.
For a spatially explicit model, movement should also most likely be part of the rule-based portion, as well as functions which require direct information exchange between individuals.
EcotoxSystems.initialize_statevars — Method
initialize_statevars(p::CVOrParamStruct)::CVOrParamStructFor initialization of ODE simulator, initialize the component vector of state variables, u, based on common oaraeter collection p.
EcotoxSystems.k_J! — Method
k_J!(spc::CVOrParamStruct)::NothingSet the maturity maintenance rate constant, assuming that the cumulative investment into maturity maintenance equals the cumulative investment into somatic maintenance (cf. DEBkiss book by Tjalling Jager).
EcotoxSystems.lineplot! — Method
lineplot(x, y, q; estimator)Plot aggregated values of y over x. Similar to seaborn.lineplot. By default, the arithmetic mean of y for every value of x is plotted, with 5th to 95th percentiles as ribbon. The plotted percentile range can be controlled via the positional argument q, which is a two-element Tuple. The arithmetic mean can be replaced with another function via the estimator argument.
EcotoxSystems.lineplot! — Method
lineplot(x, y, q; estimator)Plot aggregated values of y over x. Similar to seaborn.lineplot. By default, the arithmetic mean of y for every value of x is plotted, with 5th to 95th percentiles as ribbon. The plotted percentile range can be controlled via the positional argument q, which is a two-element Tuple. The arithmetic mean can be replaced with another function via the estimator argument.
EcotoxSystems.lineplot — Method
lineplot(x, y, q; estimator)Plot aggregated values of y over x. Similar to seaborn.lineplot. By default, the arithmetic mean of y for every value of x is plotted, with 5th to 95th percentiles as ribbon. The plotted percentile range can be controlled via the positional argument q, which is a two-element Tuple. The arithmetic mean can be replaced with another function via the estimator argument.
EcotoxSystems.link_params! — Function
link_params!(p::CVOrParamStruct, links::NamedTuple = (spc = linkfun,))::NothingApply functions to link parameters. <br>
p: ComponentVector containing parameterslinks: Component-specific functions to define links between parameters.
Examples
We can apply the link before running a simulation, in which case the link is applied to the spc component (species-level parameters).-
link_ind_params!(p) = begin
p.k_J_emb = (1-p.kappa_emb)/p.kappa_emb * p.k_M_emb
end
p = deepcopy(debkiss_defaultparams)
link_params!(p, (spc = link_ind_params!,) # apply link ahead of simulation Alternatively, we provide the links as a keyword argument to ODE_simulator, in which case the link is applied to the ind component. <br> This is useful if one of the parameters involved in the link is subject to individual variability, and we need to update the link for each simulated individual. <br>
sim = ODE_simulator(p, param_links = (ind = link_ind_params!,))EcotoxSystems.model_step! — Method
model_step!(m::AbstractIBM)::NothingGeneric definition of an individual-based model step.
EcotoxSystems.record_global! — Method
record_global!(m::AbstractIBM)::NothingStore global state variables in m.global_record.
EcotoxSystems.record_individual! — Method
record_individual!(a::AbstractIndividual, m::AbstractIBM)::NothingStore individual-level state variables in m.individual_record.
EcotoxSystems.relative_response — Method
relative_response(
sim::D,
response_vars::Vector{Symbol},
treatment_var::Symbol;
groupby_vars::Vector{Symbol} = Symbol[],
identify_control = minimum
) where D <: AbstractDataFrameCalculate relative responses.
args:
sim::AbstractDataFrame: resultsresponse_vars::Vector{Symbol}: response variables for which to calculate the relative responsestreatment_var::Symbol: Column indicating the treatment. Column values can be numerical or categorical, butidentify_controlkwarg has to be specified in the latter case
kwargs:
groupby_vars::Vector{Symbol}: relative response will be conditioned on these variables (e.g. time, separate experiments...). Empty by default.identify_control: function used to identify reference values fromtreatment_var. By default, this isminimum()(assuming that column values intreatment_varare numerical).
Example
using MechanisticEffectModels.EcotoxSystems, MechanisticEffectModels.Utils
simfunct(x) = @replicates EcotoxSystems.simulator(x) 10
sim = exposure(simfunct, Params(), [0., 100., 200.]) |>
x -> relative_response(x, [:S, :H, :R]) # -> data frame will contain columns y_S, y_H, y_R for control-normalized valuesEcotoxSystems.replicates — Method
replicates(simulator::Function, debkiss_defaultparams::CVOrParamStruct, nreps::Int64; kwargs...)Perform replicated runs of simulator with parameters debkiss_defaultparams (simulator(debkiss_defaultparams) has to be a valid function call). Analogous to @replicates, but a bit more flexible.
EcotoxSystems.set_global_statevars! — Method
set_global_statevars!(m::AbstractIBM, a::AbstractIndividual)::NothingUpdate global state variables and derivatives based on individual-level model step.
EcotoxSystems.texposure — Method
Threaded version of exposure().
EcotoxSystems.treplicates — Method
treplicates(
simulator::Function,
debkiss_defaultparams::CVOrParamStruct,
nreps::Int64;
kwargs...)Multi-threaded version of replicates.
Only useful if Julia has been started with multiple threads.
To check the number of threads, run using Base.Threads; Threads.nthreads().
In VSCode, you can use the entry "julia.NumThreads" in settings.json to set the default number of threads (searching for "julia threads" in the preferences will lead you there).
Check the Multi-threading documentation for more information.
EcotoxSystems.@replicates — Macro
@replicates(simcall::Expr, nreps::Int64)Perform replicated runs of simcall, where simcall is a call to a simulator function.
Example:
spc = SpeciesParams(Z = Truncated(Normal(1, 0.1), 0, Inf)) # initialize default parameters with variable zoom factor
sim = @replicates MechanisticEffectModels.simulator(Params(spc = spc))) 10 # execute replicated runs to simulatorIn this case, sim will contain the output of 10 replicated simulations. For each replicate, the zoom factor is sampled from a truncated Normal distribution. sim contains an additional column replicate.