EcotoxSystems.jl

Default parameter object

DEB-ODE model with arbitrary number of stressors, assuming IA to compute combined effects.

Individual-level part of the DEB-ODE model with arbitrary number of stressors, assuming IA to compute combined effects.

DEBkiss!(du, u, p, t)::Nothing

Dynamics of DEBkiss model with maturity and explicit simulation of resource dynamics.

The density of structure is ignored, and instead S^(2/3) is applied for surface-area scaling. This affects the dimension of dI_max, but has no effect on the model dynamics.

If model output is to be compared to length data, a statistical weight-lenght relationship has to be applied to the model output.

Euler!(u::ComponentVector, du::ComponentVector, dt::Real)::Nothing

Apply Euler scheme to state variables.

IBM_simulator(
    p::ComponentVector; 
    global_ode! = DEBODE_global!,
    global_rules! = default_global_rules!,
    init_global_statevars = initialize_global_statevars,
    individual_ode! = DEBkiss_individual!,
    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.

ODE_simulator(
    p::ComponentVector; 
    alg = Tsit5(),
    saveat = 1,
    reltol = 1e-6,
    model = DEBODE!,
    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 signature statevars_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 are dataframe (complete solution converted to a DataFrame) and odesol (the ODE solution object as returned by OrdinaryDiffEq). Default is dataframe.

Run a model as purely ODE-based system:

using EcotoxSystems
sim = DEB.ODE_simulator(EcotoxSystems.defaultparams(p))

function TKTD_mix_IA!(du, u, p, t)::Nothing

Mixture-TKTD for an arbitrary number of stressors, assuming Independent Action.

calc_SL_max(spc::ComponentVector)::Float64

Calculate maximum structural length slmax [m^(1/3)]

calc_S_max(spc::ComponentVector)::Float64

Calculate maximum structural mass smax [m]

Clip negative values at 0 as a continuous function, using sig.

dI_embryo(
    embryo::Float64,
    S::Float64,
    dI_max_emb::Float64,
    y_T::Float64
    )::Float64

Resource ingestion by embryos, i.e. uptake of yolk/vitellus.

Arguments

• embryo: State variable indicating whether current state is embryonic
• S: Structural mass
• dI_max_emb: Maxmimum size-specific ingestion rate of embryos
• y_T: Temperature coefficient

dS(
    kappa::Float64,
    dA::Float64, 
    dM::Float64, 
    eta_SA::Float64,
    y_G::Float64,
    eta_AS::Float64
    )::Float64

Structural growth rate.

Arguments

• kappa: Allocation fraction to somatic growth and maintenance
• dA: Assimilation rate
• dM: Somatic maintenance rate
• eta_SA: Growth efficiency (yield of somatic mass on assimilates)
• y_G: Relative response of growth efficiency
• eta_SA: Shrinking efficiency

default_global_rules!(m)

Global rule-based portion of the default model.

default_individual_rules(a::AbstractDEBIndividual, m::AbstractDEBIBM)::Nothing

Defines 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.

exposure(
    simulator::Function, 
    p::ComponentVector, 
    C_Wmat::Union{Vector{R},Matrix{R}}
) where R <: Real

Simulate constant chemical exposure with over an arbitrary number of treatments and chemical stressors, defined in C_Wmat.

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

C_Wmat = [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

C_Wmat = [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:

C_Wmat = [
    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.

filter_individuals!(m::AbstractDEBIBM)

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.

generate_individual_params(p::ComponentVector; 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.

get_global_statevars!(a::AbstractDEBIndividual, m::AbstractDEBIBM)::Nothing

Retrieve global state variables and derivatives for use in the individual-level model step.

groupedlineplot(x, y, g, q; estimator)

Version of lineplot with additional grouping variable g.

groupedlineplot(x, y, g, q; estimator)

Version of lineplot with additional grouping variable g.

groupedlineplot(x, y, g, q; estimator)

Version of lineplot with additional grouping variable g.

individual_step!(a::Agent, m::Model)::Nothing

The 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.

initialize_global_statevars(p::ComponentVector)

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.

initialize_individual_statevars(
    p::ComponentVector; 
    id = 1., 
    cohort = 0.)::ComponentVector

This function defines the individual-level state variables and their initial values for the default 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.)
)

initialize_statevars(p::ComponentVector)::ComponentVector

For initialization of ODE simulator, initialize the component vector of state variables, u, based on common oaraeter collection p.

k_J!(spc::ComponentVector)::Nothing

Set 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).

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.

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.

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.

link_params!(p::ComponentVector, links::NamedTuple = (spc = linkfun,))::Nothing

Apply functions to link parameters. <br>

• p: ComponentVector containing parameters
• links: 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(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!,))

model_step!(m::AbstractDEBIBM)::Nothing

Generic definition of an individual-based model step.

params(;kwargs...) = ComponentVector(defaultparams; kwargs...)

Initialize the default parameter set, modifying and/or adding parameter with kwargs.

Examples

# simulate the default parameters
p = params()
sim = ODE_simulator(p)  

# simulate the default parameters, but modify kappa
p = params(kappa = 0.5) 
sim = ODE_simulator

# simulate the default parameters with individual variabilty in kappa
p = params(kappa = truncated(Normal(0.5, 0.1), 0, 1))
sim = @replicates ODE_simulator(p) 10

record_global!(m::AbstractDEBIBM)::Nothing

Store global state variables in m.global_record.

record_individual!(a::AbstractDEBIndividual, m::AbstractDEBIBM)::Nothing

Store individual-level state variables in m.individual_record.

relative_response(
    sim::D, 
    response_vars::Vector{Symbol},
    treatment_var::Symbol; 
    groupby_vars::Vector{Symbol} = Symbol[],
    identify_control = minimum
    ) where D <: AbstractDataFrame

Calculate relative responses.

args:

• sim::AbstractDataFrame: results
• response_vars::Vector{Symbol}: response variables for which to calculate the relative responses
• treatment_var::Symbol: Column indicating the treatment. Column values can be numerical or categorical, but identify_control kwarg 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 from treatment_var. By default, this is minimum() (assuming that column values in treatment_var are 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 values

replicates(simulator::Function, defaultparams::ComponentVector, nreps::Int64; kwargs...)

Perform replicated runs of simulator with parameters defaultparams (simulator(defaultparams) has to be a valid function call). Analogous to @replicates, but a bit more flexible.

set_global_statevars!(m::AbstractDEBIBM, a::AbstractDEBIndividual)::Nothing

Update global state variables and derivatives based on individual-level model step.

sig(
    x::Real, 
    x_thr::Real,
    y_left::Real, 
    y_right::Real; 
    beta::Real = 30
    )::Real

Sigmoid switch function. Used to replace simple if-statements with a continuous function in ODE models.

y_left and y_right are the function values left and right of the threshold x_thr.

Threaded version of exposure().

treplicates(
    simulator::Function, 
    defaultparams::ComponentVector, 
    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.

@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 simulator

In 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.