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 in the directory 'simple_compound'. 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) L = X(2); % state 2 is body length 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.
% unpack globals FBV = glo.FBV; % dry weight of egg as fraction of dry body weight (-) KRV = glo.KRV; % part. coeff. repro buffer and structure (kg/kg) kap = glo.kap; % approximation for kappa (-) yP = glo.yP; % product of yVA and yAV (-) % unpack model parameters for the basic life history L0 = par.L0(1); % body length at start (mm) Lp = par.Lp(1); % body length at puberty (mm) Lm = par.Lm(1); % maximum body length (mm) rB = par.rB(1); % von Bertalanffy growth rate constant (1/d) Rm = par.Rm(1); % maximum reproduction rate (#/d) f = par.f(1); % scaled functional response (-) hb = par.hb(1); % background hazard rate (d-1) % unpack extra parameters for specific cases Lf = par.Lf(1); % body length at half-saturation feeding (mm) Tlag = par.Tlag(1); % lag time for start development (d) % unpack model parameters for the response to toxicants 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))
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 DEBkiss model, with toxicant effects and starvation module, expressed in terms of compound parameters, as presented in the DEBkiss e-book Section 5.5.5. (with updates that will be presented in a publication soon!).
L = max(1e-3*L0,L); % make sure that body length is not negative or almost zero (extreme shrinking may do that) if Lf > 0 % to include feeding limitation for juveniles ... f = f / (1+(Lf^3)/(L^3)); % hyperbolic relationship for f with body volume % kd = kd*f; % also reduce dominant rate by same factor? (open for discussion!) end % calculate stress factor and hazard rate s = bb*max(0,Dw-zb); % stress level for metabolic effects h = bs*max(0,Dw-zs); % hazard rate for effects on survival % Define specific stress factors s*, depending on the mode of action as % specified in the vector with switches glo.moa. Si = glo.moa * s; % vector with specific stress factors from switches for mode of action sA = min(1,Si(1)); % assimilation/feeding (maximise to 1 to avoid negative values for 1-sA) sM = Si(2); % maintenance (somatic and maturity) sG = Si(3); % growth costs sR = Si(4); % reproduction costs % calcululate the actual derivatives, with stress implemented dL = rB * ((1+sM)/(1+sG)) * (f*Lm*((1-sA)/(1+sM)) - L); % ODE for body length fR = f; % if there is no starvation, f for reproduction is the standard f % starvation rules can modify the outputs here if dL < 0 % then we are looking at starvation and need to correct things fR = (f - kap * (L/Lm) * ((1+sM)/(1-sA)))/(1-kap); % new f for reproduction alone if fR >= 0 % then we are in the first stage of starvation: 1-kappa branch can help pay maintenance dL = 0; % stop growth, but don't shrink else % we are in stage 2 of starvation and need to shrink to pay maintenance fR = 0; % nothing left for reproduction dL = (rB*(1+sM)/yP) * ((f*Lm/kap)*((1-sA)/(1+sM)) - L); % shrinking rate end end R = 0; % reproduction rate is zero, unless ... if L >= Lp % if we are above the length at puberty, reproduce R = max(0,(Rm/(1+sR)) * (fR*Lm*(L^2)*(1-sA) - (Lp^3)*(1+sM))/(Lm^3 - Lp^3)); end dRc = R; % cumulative reproduction rate dS = -(h + hb) * S; % change in survival probability (incl. background mort.) % For the damage dynamics, there are four feedback factors x* that obtain a % value based on the settings in the configuration vector glo.feedb: a % vector with switches for various feedbacks: [surface:volume on uptake, % surface:volume on elimination, growth dilution, losses with % reproduction]. Xi = glo.feedb .* [Lm/L,Lm/L,(3/L)*dL,R*FBV*KRV]; % multiply switch factor with feedbacks xu = max(1,Xi(1)); % if switch for surf:vol scaling is zero, the factor must be 1 and not 0! xe = max(1,Xi(2)); % if switch for surf:vol scaling is zero, the factor must be 1 and not 0! xG = Xi(3); % factor for growth dilution xR = Xi(4); % factor for losses with repro dDw = kd * (xu * c - xe * Dw) - (xG + xR) * Dw; % ODE for scaled damage if L <= 1e-3*L0 % if an animal has size close to zero ... dL = 0; % don't let it grow or shrink any further (to avoid numerical issues) end if t < Tlag % when we use a lag time ... dX = [0;0;0;0;0]; % set all derivatives in vector dX to zero else dX = [dDw;dL;dRc;dS]; % collect all derivatives in one vector dX end