BYOM function derivatives.m (the model in ODEs)
Syntax: dX = derivatives(t,X,par,c)
This function calculates the derivatives for the reduced GUTS model system. Note that the survival probability due to chemical stress is all calculated in call_deri.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 % mw = par.mw(1); % median of threshold distribution (used in call_deri) % bw = par.bw(1); % killing rate (used in call_deri) % Fs = par.Fs(1); % fraction spread of threshold distribution, (-) (used in call_deri) hb = par.hb(1); % background hazard rate
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.
dDw = kd * (c - Dw); % first order damage build-up from c (scaled) dS = -hb* S; % only background hazard rate % mortality due to the chemical is included in call_deri! dX = [dS;dDw]; % collect derivatives in one vector