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_debkisstox_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 WV = X(2); % state 2 is dry structural body mass cR = X(3); % state 3 is cumulative reproduction S = X(4); % state 4 is survival probability
Not enough input arguments. Error in derivatives (line 34) cV = X(1); % state 1 is the scaled internal concentration
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.
dV = glo.dV; % dry weight density of structure delM = glo.delM; % shape correction coefficient sJAm = par.sJAm(1); % specific assimilation rate sJM = par.sJM(1); % specific maintenance costs WB0 = par.WB0(1); % initial dry weight of egg Lwf = par.Lwf(1); % body length at half-saturation feeding Lwp = par.Lwp(1); % body length at puberty f = par.f(1); % scaled functional response yAV = par.yAV(1); % yield of assimilates on structure (starvation) yBA = par.yBA(1); % yield of egg buffer on assimilates yVA = par.yVA(1); % yield of structure on assimilates (growth) kap = par.kap(1); % allocation fraction to soma ke = par.ke(1); % dominant elimination rate c0 = par.c0(1); % no-effect concentration cT = par.cT(1); % tolerance concentration c0s = par.c0s(1); % no-effect conc. for survival b = par.b(1); % killing rate for survival h0 = par.h0(1); % background hazard rate (d-1) PRV = par.PRV(1); % partition coeff between repro buffer and structure (weight basis) WVp = dV * (Lwp * delM)^3; % translate puberty length todry weight WVf = dV * (Lwf * delM)^3; % translate feeding length to dry weight Tlag = par.Tlag(1); % lag time before everything starts
Calculate the derivatives
This is the actual model, specified as a system of ODEs. Note: some unused fluxes are calculated because they may be needed for future modules (e.g., respiration and ageing).
L = (WV/dV)^(1/3); % calculate volumetric length from dry weight Lm = kap*sJAm/sJM; % maximum structural body length in control (for scaling ke) if WVf > 0 sf = 1 / (1+WVf/WV); % hyperbolic relationship for f with body weight sJAm = sJAm * sf; % reduce the max assimilation rate ke = ke*sf; % also reduce elimination rate by same factor end % select what to do with maturity maintenance if glo.mat == 1 sJJ = sJM * (1-kap)/kap; % add specific maturity maintenance with the suggested value else sJJ = 0; % or ignore it end WVb = WB0 *yVA * kap; % body mass at birth (to put a stop on shrinking) s = max(cV-c0,0)/cT; % stress level for metabolic effects h = b*max(0,cV-c0s); % hazard rate for effects on survival switch glo.moa % switch for mode of action case 1 % assimilation/feeding sJAm = sJAm * max(0,1-s); % effect on assimilation, alternative: /(1+s) increase in handling time case 2 % maintenance sJM = sJM * (1+s); % effect on somatic maintenance sJJ = sJJ * (1+s); % effect on maturity maintenance case 3 % growth costs yVA = yVA / (1+s); % effect on growth costs, alternative: *(1-s) case 4 % repro costs yBA = yBA / (1+s); % effect on repro costs case 5 % hazard to embryo yBA = yBA * exp(-s); % hazard on embryo case 6 % growth and repro costs yVA = yVA / (1+s); % effect on growth costs, alternative: *(1-s) yBA = yBA / (1+s); % add effect on repro costs end JA = f * sJAm * L^2; % assimilation flux JM = sJM * L^3; % somatic maintenance flux JV = yVA * (kap*JA-JM); % growth flux if WV < WVp % below size at puberty ... JR = 0; % no reproduction flux JJ = sJJ * L^3; % maturity maintenance flux % JH = (1-kap) * JA - JJ; % maturation flux (not used!) else % above size at puberty ... JJ = sJJ * (WVp/dV); % maturity maintenance flux JR = (1-kap) * JA - JJ; % reproduction flux end % starvation rules may override these fluxes if kap * JA < JM % allocated flux to soma cannot pay maintenance if JA >= JM + JJ % but still enough total assimilates to pay both maintenances JV = 0; % stop growth if WV >= WVp % for adults ... JR = JA - JM - JJ; % repro buffer gets what's left else % JH = JA - JM - JJ; % maturation gets what's left (not used!) end elseif JA >= JM % only enough to pay somatic maintenance JV = 0; % stop growth JR = 0; % stop reproduction % JH = 0; % stop maturation for juveniles (not used) JJ = JA - JM; % maturity maintenance flux gets what's left else % we need to shrink JR = 0; % stop reproduction JJ = 0; % stop paying maturity maintenance % JH = 0; % stop maturation for juveniles (not used) JV = (JA - JM) / yAV; % shrink; pay somatic maintenance from structure end end % calcululate the derivatives dWV = JV; % change in body mass dcR = yBA * JR / WB0; % continuous reproduction flux dcV = ke * (Lm/L) * (c-cV) - cV * dWV/WV; % change in scaled internal conc. dcV = dcV - yBA * JR * PRV * cV / WV; % include continuous losses due to repro dS = -(h + h0)* S; % change in survival probability (incl. background mort.) if WV < WVb/4 % dont shrink below a quarter of the size at birth dWV = 0; % simply stop shrinking ... end if t < Tlag dX = [0;0;0;0]; % set all derivatives in vector dX to zero else dX = [dcV;dWV;dcR;dS]; % collect all derivatives in one vector dX end