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:
- 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.
function dX = derivatives(t,X,par,c)
global glo % allow for global parameters in structure glo (handy for switches)
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
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 = par.mw(1); % median of threshold distribution, ug/L bw = par.bw(1); % 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!') end dX = [dS;dDw]; % collect both derivatives in one vector