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:

• t is the time point, provided by the ODE solver
• X is a vector with the previous value of the states
• par is the parameter structure
• c is the external concentration (or scenario number)

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.

LICENSE.txt file in the root directory of BYOM.

Start

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   = par.bb(1);   % effect strength energy-budget effects (1/[C])
zs   = par.zs(1);   % effect threshold survival ([C])
bs   = par.bs(1);   % 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)
end
end

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!)
end

% 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
end
end

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));
end

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)
end

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