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