# 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*.

- Author: Tjalling Jager
- Date: August 2015
- Web support: http://www.debtox.info/byom.html
- Back to index walkthrough_debtox.html

## Contents

## Start

```
function dX = derivatives(t,X,par,c)
```

global glo % allow for global parameters in structure glo (handy for switches)

## Unpack states

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.

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

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

## Calculate the derivatives

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