BYOM function derivatives.m (the model in ODEs)

Syntax: dX = derivatives(t,X,par,c)

This function calculates the derivatives for the model system. It is linked to the script files in the directory 'simple_compound'. As input, it gets:

Time t and scenario name c are handed over as single numbers by call_deri.m (you do not have to use them in this function). Output dX (as vector) provides the differentials for each state at t.

Copyright (c) 2012-2019, Tjalling Jager, all rights reserved.
This source code is licensed under the MIT-style license found in the
LICENSE.txt file in the root directory of BYOM.



function dX = derivatives(t,X,par,c)
global glo   % allow for global parameters in structure glo (handy for switches)

Unpack states

The state variables enter this function in the vector X. Here, we give them a more handy name.

Dw = X(1); % state 1 is the scaled damage (referenced to water)
L  = X(2); % state 2 is body length
Rc = X(3); % state 3 is cumulative reproduction (not used)
S  = X(4); % state 4 is survival probability

Unpack parameters

The parameters enter this function in the structure par. The names in the structure are the same as those defined in the byom script file. The 1 between parentheses is needed as each parameter has 5 associated values.

% unpack globals
FBV  = glo.FBV;     % dry weight of egg as fraction of dry body weight (-)
KRV  = glo.KRV;     % part. coeff. repro buffer and structure (kg/kg)
kap  = glo.kap;     % approximation for kappa (-)
yP   = glo.yP;      % product of yVA and yAV (-)

% unpack model parameters for the basic life history
L0   = par.L0(1);   % body length at start (mm)
Lp   = par.Lp(1);   % body length at puberty (mm)
Lm   = par.Lm(1);   % maximum body length (mm)
rB   = par.rB(1);   % von Bertalanffy growth rate constant (1/d)
Rm   = par.Rm(1);   % maximum reproduction rate (#/d)
f    = par.f(1);    % scaled functional response (-)
hb   = par.hb(1);   % background hazard rate (d-1)

% unpack extra parameters for specific cases
Lf   = par.Lf(1);   % body length at half-saturation feeding (mm)
Tlag = par.Tlag(1); % lag time for start development (d)

% unpack model parameters for the response to toxicants
kd   = par.kd(1);   % dominant rate constant (d-1)
zb   = par.zb(1);   % effect threshold energy budget ([C])
bb   =;   % effect strength energy-budget effects (1/[C])
zs   = par.zs(1);   % effect threshold survival ([C])
bs   =;   % effect strength survival (1/([C] d))

Extract correct exposure for THIS time point

Allow for external concentrations to change over time, either continuously, or in steps, or as a static renewal with first-order disappearance. For constant exposure, the code in this section is skipped (and could also be removed).

if isfield(glo,'int_scen') % if it exists: use it to derive current external conc.
    if ismember(c,glo.int_scen) % is c in the scenario range global?
        c = make_scen(-1,c,t); % use make_scen again to derive actual exposure concentration
        % the -1 lets make_scen know we are calling from derivatives (so need one c)

Calculate the derivatives

This is the actual model, specified as a system of ODEs. This is the DEBkiss model, with toxicant effects and starvation module, expressed in terms of compound parameters, as presented in the DEBkiss e-book Section 5.5.5. (with updates that will be presented in a publication soon!).

L = max(1e-3*L0,L); % make sure that body length is not negative or almost zero (extreme shrinking may do that)

if Lf > 0 % to include feeding limitation for juveniles ...
    f  = f / (1+(Lf^3)/(L^3)); % hyperbolic relationship for f with body volume
    % kd = kd*f; % also reduce dominant rate by same factor? (open for discussion!)

% calculate stress factor and hazard rate
s = bb*max(0,Dw-zb); % stress level for metabolic effects
h = bs*max(0,Dw-zs); % hazard rate for effects on survival

% Define specific stress factors s*, depending on the mode of action as
% specified in the vector with switches glo.moa.
Si = glo.moa * s; % vector with specific stress factors from switches for mode of action
sA = min(1,Si(1)); % assimilation/feeding (maximise to 1 to avoid negative values for 1-sA)
sM = Si(2); % maintenance (somatic and maturity)
sG = Si(3); % growth costs
sR = Si(4); % reproduction costs

% calcululate the actual derivatives, with stress implemented
dL = rB * ((1+sM)/(1+sG)) * (f*Lm*((1-sA)/(1+sM)) - L); % ODE for body length

fR = f; % if there is no starvation, f for reproduction is the standard f
% starvation rules can modify the outputs here
if dL < 0 % then we are looking at starvation and need to correct things
    fR = (f - kap * (L/Lm) * ((1+sM)/(1-sA)))/(1-kap); % new f for reproduction alone
    if fR >= 0 % then we are in the first stage of starvation: 1-kappa branch can help pay maintenance
        dL = 0; % stop growth, but don't shrink
    else % we are in stage 2 of starvation and need to shrink to pay maintenance
        fR = 0; % nothing left for reproduction
        dL = (rB*(1+sM)/yP) * ((f*Lm/kap)*((1-sA)/(1+sM)) - L); % shrinking rate

R  = 0; % reproduction rate is zero, unless ...
if L >= Lp % if we are above the length at puberty, reproduce
    R = max(0,(Rm/(1+sR)) * (fR*Lm*(L^2)*(1-sA) - (Lp^3)*(1+sM))/(Lm^3 - Lp^3));

dRc = R; % cumulative reproduction rate
dS  = -(h + hb) * S; % change in survival probability (incl. background mort.)

% For the damage dynamics, there are four feedback factors x* that obtain a
% value based on the settings in the configuration vector glo.feedb: a
% vector with switches for various feedbacks: [surface:volume on uptake,
% surface:volume on elimination, growth dilution, losses with
% reproduction].
Xi = glo.feedb .* [Lm/L,Lm/L,(3/L)*dL,R*FBV*KRV]; % multiply switch factor with feedbacks
xu = max(1,Xi(1)); % if switch for surf:vol scaling is zero, the factor must be 1 and not 0!
xe = max(1,Xi(2)); % if switch for surf:vol scaling is zero, the factor must be 1 and not 0!
xG = Xi(3); % factor for growth dilution
xR = Xi(4); % factor for losses with repro

dDw = kd * (xu * c - xe * Dw) - (xG + xR) * Dw; % ODE for scaled damage

if L <= 1e-3*L0 % if an animal has size close to zero ...
    dL = 0; % don't let it grow or shrink any further (to avoid numerical issues)

if t < Tlag % when we use a lag time ...
    dX = [0;0;0;0;0]; % set all derivatives in vector dX to zero
    dX = [dDw;dL;dRc;dS]; % collect all derivatives in one vector dX