BYOM function derivatives.m (the model in ODEs)
Syntax: dX = derivatives(t,X,par,c,glo)
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)
- glo is the structure with information (normally global)
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.
- Author: Tjalling Jager
- Date: November 2021
- Web support: http://www.debtox.info/byom.html
- Back to index walkthrough_debtox2019.html
Copyright (c) 2012-2021, 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.
Contents
Start
function dX = derivatives(t,X,par,c,glo)
% NOTE: glo is now no longer a global here, but passed on in the function % call from call_deri. That saves 20% calculation time!
Unpack states
The state variables enter this function in the vector X. Here, we give them a more handy name.
X = max(X,0); % ensure the states cannot become negative due to numerical issues Dw = X(1); % state is the scaled damage (referenced to water) L = X(2); % state is body length % Rc = X(3); % state is cumulative reproduction (not used) S = X(4); % state is survival probability
Unpack parameters
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) Lj = par.Lj(1); % body length at which acceleration stops (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). Note that glo.timevar is used to signal that there is a time-varying concentration. This option is set in call_deri.
if glo.timevar(1) == 1 % if we are warned that we have a time-varying concentration ... c = read_scen(-1,c,t,glo); % use read_scen to derive actual exposure concentration % Note: this has glo as input to the function to save time! % the -1 lets read_scen know we are calling from derivatives (so need one c) 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 publication (Jager, 2020).
L = max(1e-3*L0,L); % make sure that body length is not negative or almost zero (extreme shrinking may do that) % This should not be needed as shrinking is limited at the bottom of this % function. 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 if Lj > 0 % to include acceleration until metamorphosis ... f = f * min(1,L/Lj); % this implies lower f for L<Lj 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 sH = Si(5); % also include hazard to reproduction % 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,(exp(-sH)*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 = Xi(1); % factor for surf:vol scaling uptake rate xe = Xi(2); % factor for surf:vol scaling elimination rate xG = Xi(3); % factor for growth dilution xR = Xi(4); % factor for losses with repro % if switch for surf:vol scaling is zero, the factor must be 1 and not 0! % Note: this was previously done with a max-to-1 command. However, that is % a bit dangerous when modifying the model. if Xi(1) == 0 xu = 1; end if Xi(2) == 0 xe = 1; end xG = max(0,xG); % NOTE NOTE: reverse growth dilution (concentration by shrinking) is now % turned OFF as it leads to runaway situations that lead to failure of the % ODE solvers. However, this needs some further thought! dDw = kd * (xu * c - xe * Dw) - (xG + xR) * Dw; % ODE for scaled damage if L <= 0.5 * L0 % if an animal has size less than half the start size ... dL = 0; % don't let it grow or shrink any further (to avoid numerical issues) end dX = zeros(size(X)); % initialise with zeros in right format if t >= Tlag % when we are past the lag time ... dX = [dDw;dL;dRc;dS]; % collect all derivatives in one vector dX end