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

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.

cV = X(1); % state 1 is the scaled internal concentration
WV = X(2); % state 2 is dry structural body mass
cR = X(3); % state 3 is cumulative reproduction
S  = X(4); % state 4 is survival probability
Not enough input arguments.

Error in derivatives (line 34)
cV = X(1); % state 1 is the scaled internal concentration

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
delM  = glo.delM;     % shape correction coefficient

sJAm  = par.sJAm(1);  % specific assimilation rate
sJM   = par.sJM(1);   % specific maintenance costs
WB0   = par.WB0(1);   % initial dry weight of egg
Lwf   = par.Lwf(1);   % body length at half-saturation feeding
Lwp   = par.Lwp(1);   % body length at puberty
f     = par.f(1);     % scaled functional response

yAV   = par.yAV(1);   % yield of assimilates on structure (starvation)
yBA   = par.yBA(1);   % yield of egg buffer on assimilates
yVA   = par.yVA(1);   % yield of structure on assimilates (growth)
kap   = par.kap(1);   % allocation fraction to soma

ke    = par.ke(1);    % dominant elimination rate
c0    = par.c0(1);    % no-effect concentration
cT    = par.cT(1);    % tolerance concentration
c0s   = par.c0s(1);   % no-effect conc. for survival
b     = par.b(1);     % killing rate for survival
h0    = par.h0(1);    % background hazard rate (d-1)
PRV   = par.PRV(1);   % partition coeff between repro buffer and structure (weight basis)

WVp = dV * (Lwp * delM)^3; % translate puberty length todry  weight
WVf = dV * (Lwf * delM)^3; % translate feeding length to dry weight

Tlag = par.Tlag(1); % lag time before everything starts

Calculate the derivatives

This is the actual model, specified as a system of ODEs. 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 volumetric length from dry weight
Lm = kap*sJAm/sJM; % maximum structural body length in control (for scaling ke)

if WVf > 0
    sf = 1 / (1+WVf/WV); % hyperbolic relationship for f with body weight
    sJAm = sJAm * sf; % reduce the max assimilation rate
    ke = ke*sf; % also reduce elimination rate by same factor
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
end
WVb = WB0 *yVA * kap; % body mass at birth (to put a stop on shrinking)

s = max(cV-c0,0)/cT; % stress level for metabolic effects
h = b*max(0,cV-c0s); % hazard rate for effects on survival

switch glo.moa % switch for mode of action
    case 1 % assimilation/feeding
        sJAm = sJAm * max(0,1-s); % effect on assimilation, alternative: /(1+s) increase in handling time
    case 2 % maintenance
        sJM = sJM * (1+s);   % effect on somatic maintenance
        sJJ = sJJ * (1+s);   % effect on maturity maintenance
    case 3 % growth costs
        yVA = yVA / (1+s);   % effect on growth costs, alternative: *(1-s)
    case 4 % repro costs
        yBA = yBA / (1+s);   % effect on repro costs
    case 5 % hazard to embryo
        yBA = yBA * exp(-s); % hazard on embryo
    case 6 % growth and repro costs
        yVA = yVA / (1+s);   % effect on growth costs, alternative: *(1-s)
        yBA = yBA / (1+s);   % add effect on repro costs

end

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

% calcululate the derivatives
dWV = JV;             % change in body mass
dcR = yBA * JR / WB0; % continuous reproduction flux
dcV = ke * (Lm/L) * (c-cV) - cV * dWV/WV; % change in scaled internal conc.
dcV = dcV - yBA * JR * PRV * cV / WV; % include continuous losses due to repro

dS = -(h + h0)* S; % change in survival probability (incl. background mort.)

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

if t < Tlag
    dX = [0;0;0;0]; % set all derivatives in vector dX to zero
else
    dX = [dcV;dWV;dcR;dS]; % collect all derivatives in one vector dX
end