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



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.

S  = X(1); % state 1 is the survival probability at previous time point
Dw = X(2); % state 2 is the scaled damage at previous time point

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.

kd   = par.kd(1);   % dominant rate constant, d-1
mw   =;   % median of threshold distribution, ug/L
bw   =;   % killing rate, L/ug/d
Fs   = par.Fs(1);   % fraction spread of threshold distribution (used in call_deri)
hb   = par.hb(1);   % background hazard rate, d-1

Calculate the derivatives

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

Note: scaled damage and background hazard are calculated here for every special case. However, survival is done here only for the pure SD model. Pure IT is calculated in call_deri.m.

dDw = kd * (c - Dw); % first order damage build-up from c (scaled)

switch glo.sel
    case 1 % stochastic death
        h  = bw * max(0,Dw-mw); % calculate the hazard rate
        dS = -(h + hb)* S;      % change in survival probability (incl. background mort.)
    case 2 % individual tolerance
        dS = -hb* S;            % only background hazard rate
        % mortality due to the chemical is included in call_deri!
    case 3
        error('Use the GUTS models in the "full" or "reduced" folder for the IT+SD model!')

dX = [dS;dDw]; % collect both derivatives in one vector