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 byom_debkisstox_daphnia.m in the directory 'compound_model'. 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.

Contents

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)
WV = X(2); % state 2 is dry structural body mass
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.

dV   = glo.dV;      % dry weight density of structure (mg/mm3)
delM = glo.delM;    % shape correction coefficient (-)
yAV  = glo.yAV;     % yield of assimilates on structure (starvation) (-)
yBA  = glo.yBA;     % yield of egg buffer on assimilates (-)
yVA  = glo.yVA;     % yield of structure on assimilates (growth) (-)

LpM  = par.LpM(1);  % actual body length at puberty (mm)
LmM  = par.LmM(1);  % actual maximum body length (mm)
rB   = par.rB(1);   % von Bertalanffy growth rate constant (1/d)
Rm   = par.Rm(1);   % maximum reproduction rate (#/d)
WB0  = par.WB0(1);  % initial dry weight of egg (mg)
f    = par.f(1);    % scaled functional response (-)
hb   = par.hb(1);   % background hazard rate (d-1)

LfM  = par.LfM(1);  % actual body length at half-saturation feeding (mm)
Tlag = par.Tlag(1); % lag time before everything starts (d)

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))
KRV  = 1;           % part. coeff. repro buffer and structure (kg/kg)

WVp = dV * (LpM * delM)^3; % translate puberty length to dry weight
WVf = dV * (LfM * delM)^3; % translate feeding-saturation length to dry weight
Lm  = delM * LmM; % volumetric maximum length
Lp  = delM * LpM; % volumetric length at puberty

% recalculate compound to primary parameters (see DEBkiss e-book, Section 5.5.2)
sJM  = rB * (3*dV)/yVA; % specific maintenance rate (mg/mm^3/d)
sJAm = Rm * (WB0/yBA) * Lm/(Lm^3-Lp^3) + sJM * Lm; % specific assimilation rate (mg/mm^2/d)
kap  = Lm * sJM/sJAm; % allocation fraction to soma (-)

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 full DEBkiss model, expressed in primary parameters, extended with damage dynamics and effects. See DEBkiss e-book Chapter 5. Even though the input parameters are now compound parameters, they have been translated to primary ones. Note: some unused fluxes are calculated because they may be needed for future modules (e.g., respiration and ageing).

L  = (WV/dV)^(1/3); % calculate current volumetric length from dry weight

if WVf > 0 % to include feeding limitation for juveniles ...
    sf   = 1 / (1+WVf/WV); % hyperbolic relationship for f with body weight
    sJAm = sJAm * sf; % reduce the max assimilation rate
    kd   = kd*sf; % also reduce elimination rate by same factor (open for discussion!)
end

% select what to do with maturity maintenance
if glo.mat == 1
    sJJ = sJM * (1-kap)/kap; % add specific maturity maintenance with the suggested value
else
    sJJ = 0; % or ignore it completely (simplest DEBkiss model)
end
WVb = WB0 * yVA * kap; % body mass at birth (to put a stop on shrinking)

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

% for the mode of action, there is a choice ...
switch glo.moa % switch for mode of action (feel free to add more options, such as an effect on kappa)
    case 1 % assimilation/feeding
        sJAm = sJAm * max(0,1-s); % effect on assimilation, alternative: /(1+s) increase in handling time
    case 2 % maintenance costs (somatic and maturity)
        sJM = sJM * (1+s);   % effect on somatic maintenance
        sJJ = sJJ * (1+s);   % effect on maturity maintenance
    case 3 % growth and repro costs
        yVA = yVA / (1+s);   % effect on growth costs, alternative: *(1-s)
        yBA = yBA / (1+s);   % add effect on repro costs
    case 4 % repro costs
        yBA = yBA / (1+s);   % effect on repro costs
end

% calculate the main fluxes
JA = f * sJAm * L^2;          % assimilation flux
JM = sJM * L^3;               % somatic maintenance flux
JV = yVA * (kap*JA-JM);       % growth flux

if WV < WVp                   % below size at puberty ...
    JR = 0;                   % no reproduction flux
    JJ = sJJ * L^3;           % maturity maintenance flux
    % JH = (1-kap) * JA - JJ; % maturation flux (not used!)
else                          % above size at puberty ...
    JJ = sJJ * (WVp/dV);      % maturity maintenance flux
    JR = (1-kap) * JA - JJ;   % reproduction flux
end

% starvation rules may override these fluxes
if kap * JA < JM      % allocated flux to soma cannot pay maintenance
    if JA >= JM + JJ  % but still enough total assimilates to pay both maintenances
        JV = 0;       % stop growth
        if WV >= WVp  % for adults ...
            JR = JA - JM - JJ; % repro buffer gets what's left
        else
            % JH = JA - JM - JJ; % maturation gets what's left (not used!)
        end
    elseif JA >= JM   % only enough to pay somatic maintenance
        JV = 0;       % stop growth
        JR = 0;       % stop reproduction
        % JH = 0;     % stop maturation for juveniles (not used)
        JJ = JA - JM; % maturity maintenance flux gets what's left
    else              % we need to shrink
        JR = 0;       % stop reproduction
        JJ = 0;       % stop paying maturity maintenance
        % JH = 0;     % stop maturation for juveniles (not used)
        JV = (JA - JM) / yAV; % shrink; pay somatic maintenance from structure
    end
end

% calculate the derivatives
dWV = JV;             % change in body mass
dRc = yBA * JR / WB0; % continuous reproduction flux
dS  = -(h + hb)* S;   % change in survival probability (incl. background mort.)

% for the change in damage dynamics, there is a choice ...
switch glo.feedb
    case 1 % effects of: surface:volume, growth dilution, reproduction
        dDw = kd * (Lm/L) *(c-Dw) - Dw * dWV/WV - dRc * (WB0/WV) * KRV * Dw;
    case 2 % effects of: surface:volume, growth dilution
        dDw = kd * (Lm/L) *(c-Dw) - Dw * dWV/WV;
    case 3 % effects of: growth dilution
        dDw = kd * (c-Dw) - Dw * dWV/WV;
    case 4 % no feedback from traits on TK
        dDw = kd * (c-Dw);
end

% to avoid shrinking to lead to negative body size ...
if WV < WVb/4   % dont shrink below a quarter of the size at birth
    dWV = 0;    % simply stop shrinking ...
end

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