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_debtox_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
e  = X(2); % state 2 is the scaled reserve density
L  = X(3); % state 3 is body length
Rc = X(4); % state 4 is cumulative reproduction
S  = X(5); % state 5 is the survival probability at previous time point
Error using derivatives (line 34)
Not enough input arguments.

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.

ke  = par.ke(1);  % elimination rate constant, d-1
c0  = par.c0(1);  % no-effect concentration sub-lethal
cT  = par.cT(1);  % tolerance concentration
c0s = par.c0s(1); % no-effect concentration survival
b   = par.b(1);   % killing rate

L0 = par.L0(1);  % initial body length (mm)
Lp = par.Lp(1);  % length at puberty (mm)
Lm = par.Lm(1);  % maximum length (mm)
rB = par.rB(1);  % von Bertalanffy growth rate
Rm = par.Rm(1);  % maximum reproduction rate
g  = par.g(1);   % energy-investment ratio
f  = par.f(1);   % scaled functional response
h0 = par.h0(1);  % background hazard rate

Calculate the derivatives

This is the actual model, specified as a system of ODEs:

kM = rB * 3 * (1+g)/g;    % maintenance rate coefficient
v  = Lm * kM * g;         % energy conductance
s  = (1/cT)*max(0,cV-c0); % stress factor

switch glo.moa % stress effect according to selected mechanism
    case 1 % effect on assimilation/feeding
        f = f*max(0,1-s);
    case 2 % effect on maintenance costs
        kM = kM * (1+s);
        Rm = Rm * (1+s);
    case 3 % effect on growth costs
        g  = g * (1+s);
        kM = kM /(1+s);
    case 4 % effect on repro costs
        Rm = Rm / (1+s);
    case 5 % hazard to early embryo
        Rm = Rm * exp(-s);
    case 6 % costs for growth and repro
        g  = g * (1+s);
        kM = kM /(1+s);
        Rm = Rm / (1+s);
end

de = (f-e)*v/L; % change in scaled energy density
dL = (kM*g/(3*(e+g))) * (e*v/(kM*g) - L); % change in length
if L<Lp % before length at puberty is reached ...
    R = 0; % no reproduction
else
    R = (Rm/(Lm^3-Lp^3))*((v*(L^2)/kM + L^3)*e/(e+g) - Lp^3); % repro rate
    R = max(0,R); % make sure it is always positive
end
dRc = R; % change in cumulative repro is the repro rate
dcV = ke*(Lm/L)*(c-cV) - cV*(3/L)*dL; % change in scaled internal conc.

h  = b * max(0,cV-c0s); % calculate the hazard rate
dS = -(h + h0)* S;      % change in survival probability (incl. background mort.)

dX = [dcV;de;dL;dRc;dS]; % collect all derivatives in one vector dX