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 in the directory 'compound_model'. 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.
Copyright (c) 2012-2019, Tjalling Jager, all rights reserved. This source code is licensed under the MIT-style license found in the LICENSE.txt file in the root directory of BYOM.
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.
Dw = X(1); % state 1 is the scaled damage (referenced to water) WV = X(2); % state 2 is dry structural body mass Rc = X(3); % state 3 is cumulative reproduction (not used) S = X(4); % state 4 is survival probability
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 (mg/mm3) delM = glo.delM; % shape correction coefficient (-) yAV = glo.yAV; % yield of assimilates on structure (starvation) (-) yBA = glo.yBA; % yield of egg buffer on assimilates (-) yVA = glo.yVA; % yield of structure on assimilates (growth) (-) LpM = par.LpM(1); % actual body length at puberty (mm) LmM = par.LmM(1); % actual maximum body length (mm) rB = par.rB(1); % von Bertalanffy growth rate constant (1/d) Rm = par.Rm(1); % maximum reproduction rate (#/d) WB0 = par.WB0(1); % initial dry weight of egg (mg) f = par.f(1); % scaled functional response (-) hb = par.hb(1); % background hazard rate (d-1) LfM = par.LfM(1); % actual body length at half-saturation feeding (mm) Tlag = par.Tlag(1); % lag time before everything starts (d) kd = par.kd(1); % dominant rate constant (d-1) zb = par.zb(1); % effect threshold energy budget ([C]) bb = par.bb(1); % effect strength energy-budget effects (1/[C]) zs = par.zs(1); % effect threshold survival ([C]) bs = par.bs(1); % effect strength survival (1/([C] d)) KRV = 1; % part. coeff. repro buffer and structure (kg/kg) WVp = dV * (LpM * delM)^3; % translate puberty length to dry weight WVf = dV * (LfM * delM)^3; % translate feeding-saturation length to dry weight Lm = delM * LmM; % volumetric maximum length Lp = delM * LpM; % volumetric length at puberty % recalculate compound to primary parameters (see DEBkiss e-book, Section 5.5.2) sJM = rB * (3*dV)/yVA; % specific maintenance rate (mg/mm^3/d) sJAm = Rm * (WB0/yBA) * Lm/(Lm^3-Lp^3) + sJM * Lm; % specific assimilation rate (mg/mm^2/d) kap = Lm * sJM/sJAm; % allocation fraction to soma (-)
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. This is the full DEBkiss model, expressed in primary parameters, extended with damage dynamics and effects. See DEBkiss e-book Chapter 5. Even though the input parameters are now compound parameters, they have been translated to primary ones. 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 current volumetric length from dry weight if WVf > 0 % to include feeding limitation for juveniles ... sf = 1 / (1+WVf/WV); % hyperbolic relationship for f with body weight sJAm = sJAm * sf; % reduce the max assimilation rate kd = kd*sf; % also reduce elimination rate by same factor (open for discussion!) 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 completely (simplest DEBkiss model) end WVb = WB0 * yVA * kap; % body mass at birth (to put a stop on shrinking) s = bb*max(0,Dw-zb); % stress level for metabolic effects h = bs*max(0,Dw-zs); % hazard rate for effects on survival % for the mode of action, there is a choice ... switch glo.moa % switch for mode of action (feel free to add more options, such as an effect on kappa) case 1 % assimilation/feeding sJAm = sJAm * max(0,1-s); % effect on assimilation, alternative: /(1+s) increase in handling time case 2 % maintenance costs (somatic and maturity) sJM = sJM * (1+s); % effect on somatic maintenance sJJ = sJJ * (1+s); % effect on maturity maintenance case 3 % growth and repro costs yVA = yVA / (1+s); % effect on growth costs, alternative: *(1-s) yBA = yBA / (1+s); % add effect on repro costs case 4 % repro costs yBA = yBA / (1+s); % effect on repro costs end % calculate the main fluxes 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 % calculate the derivatives dWV = JV; % change in body mass dRc = yBA * JR / WB0; % continuous reproduction flux dS = -(h + hb)* S; % change in survival probability (incl. background mort.) % for the change in damage dynamics, there is a choice ... switch glo.feedb case 1 % effects of: surface:volume, growth dilution, reproduction dDw = kd * (Lm/L) *(c-Dw) - Dw * dWV/WV - dRc * (WB0/WV) * KRV * Dw; case 2 % effects of: surface:volume, growth dilution dDw = kd * (Lm/L) *(c-Dw) - Dw * dWV/WV; case 3 % effects of: growth dilution dDw = kd * (c-Dw) - Dw * dWV/WV; case 4 % no feedback from traits on TK dDw = kd * (c-Dw); end % to avoid shrinking to lead to negative body size ... if WV < WVb/4 % dont shrink below a quarter of the size at birth dWV = 0; % simply stop shrinking ... end if t < Tlag % when we use a lag time ... dX = [0;0;0;0]; % set all derivatives in vector dX to zero else dX = [dDw;dWV;dRc;dS]; % collect all derivatives in one vector dX end