BYOM, byom_debkisstox_daphnia.m

BYOM is a General framework for simulating model systems in terms of ordinary differential equations (ODEs). The model itself needs to be specified in derivatives.m, and call_deri.m may need to be modified to the particular problem as well. The files in the engine directory are needed for fitting and plotting. Results are shown on screen but also saved to a log file (results.out).

The model: DEBkiss with modules to include toxicokinetics and toxic effects, as presented in detail in the DEBkiss book (see http://www.debtox.info/book_debkiss.html). Note that the TK module includes losses due to reproduction (which can be changed by changing parameter PRV).

This script: The water flea Daphnia magna exposed to fluoranthene; data set from Jager et al (2010), http://dx.doi.org/10.1007/s10646-009-0417-z. Published as case study in Jager & Zimmer (2012). Simultaneous fit on growth, reproduction and survival. Parameter estimates differ slightly from the ones in Jager & Zimmer as we here take the two controls as seperate treatments (with the same parameters).

Contents

Initial things

Make sure that this script is in a directory somewhere below the BYOM folder.

clear, clear global % clear the workspace and globals
global DATA W X0mat % make the data set and initial states global variables
global glo          % allow for global parameters in structure glo
global pri zvd      % global structures for optional priors and zero-variate data
diary off           % turn of the diary function (if it is accidentaly on)
% set(0,'DefaultFigureWindowStyle','docked'); % collect all figure into one window with tab controls
set(0,'DefaultFigureWindowStyle','normal'); % separate figure windows

pathdefine % set path to the BYOM/engine directory
glo.basenm  = mfilename; % remember the filename for THIS file for the plots
glo.saveplt = 0; % save all plots as (1) Matlab figures, (2) JPEG file or (3) PDF (see all_options.txt)

The data set

Data are entered in matrix form, time in rows, scenarios (exposure concentrations) in columns. First column are the exposure times, first row are the concentrations or scenario numbers. The number in the top left of the matrix indicates how to calculate the likelihood:

% scaled internal concentrations
DATA{1} = [0]; % there are never data for this state

% body length (mm) on each observation time (d), concentrations in uM
DATA{2} = [1 0 0 0.213	0.426	0.853
    0	0.88	0.88	0.88	0.88	0.88
    2	1.38	1.34	1.4     1.44	1.32
    4	1.96	1.72	1.9     1.88	1.74
    6	2.28	2.38	2.14	2.14	2.08
    8	2.34	2.52	2.56	2.36	2.16
    10	2.64	2.56	2.46	2.46	2.48
    12	2.66	2.56	2.6     2.58	2.7
    14	2.68	2.7     2.76	2.78	2.8
    16	2.88	2.78	2.82	NaN     NaN
    18	2.94	2.84	2.92	3.08	2.9
    20	3.06	3.02	3.02	3.22	2.76
    21	3.26	3.04	3.06	3       2.84];

W{2} = 5 * ones(size(DATA{2})-1); % weights: 5 animals in each treatment

% cumulative reproduction (nr. offspring) (mm) on each observation time (d), concentrations in uM
DATA{3} = [1 0 0 	0.213	0.426	0.853
    0	0.0	0.0	0.0	0.0	0.0
    2	0.0	0.0	0.0	0.0	0.0
    4	0.0	0.0	0.0	0.0	0.0
    6	0.0	0.0	0.0	0.0	0.0
    8	1.2     2.1     0.0     0.0     0.0
    10	6.5     7.7     10.0	0.8     1.7
    12	17.9	25.4	16.3	2.0     1.7
    14	31.8	29.5	33.3	6.4     1.7
    16	50.6	48.8	56.2	9.4     1.7
    18	61.3	62.5	62.5	13.8	1.7
    20	81.0	76.4	64.1	16.4	1.7
    21	81.0	76.4	64.1	16.4	1.7];

W{3} = [10	10	10	10	10
	10	10	10	10	10
	10	10	10	10	10
	10	10	10	10	10
	10	10	10	10	10
	10	10	10	10	9
	10	10	10	10	7
	10	10	10	10	5
	10	10	10	10	5
	10	9	10	10	3
	10	9	10	9	3
	10	8	10	9	2]; % weights: nr. of mothers alive

% in the next part, I subtract 2.5 days from the time vector for the repro
% data; reason is that the eggs are produced some 2.5 days before the
% offspring are counted.
DATA{3}(2:end,1) = DATA{3}(2:end,1) - 2.5;
DATA{3}([2 3],:) = [];
W{3}([1 2],:) = [];

% survivors on each observation time (d), concentrations in uM
DATA{4} = [-1	0	0	0.213	0.426	0.853
    0	10	10	10	10	10
    2	10	10	10	10	10
    4	10	10	10	10	10
    6	10	10	10	10	10
    8	10	10	10	10	10
    10	10	10	10	10	9
    12	10	10	10	10	7
    14	10	10	10	10	5
    16	10	10	10	10	5
    18	10	9	10	10	3
    20	10	9	10	9	3
    21	10	8	10	9	2];

% if weight factors are not specified, ones are assumed in start_calc.m

Initial values for the state variables

Initial states, scenarios in columns, states in rows. First row are the 'names' of all scenarios.

X0mat = [0 0.213 0.426 0.853   % the scenarios (here nominal concentrations)
         0     0     0     0   % initial values state 1
         0     0     0     0   % initial values state 2 (overwritten by L0)
         0     0     0     0   % initial values state 3
         1     1     1     1]; % initial values state 4

Initial values for the model parameters

Model parameters are part of a 'structure' for easy reference.

% global parameters as part of the structure glo
glo.delM  = 0.27; % shape corrector (used in call_deri.m)
glo.dV    = 0.2;  % dry weight density (used in call_deri.m)
glo.len   = 1;    % switch to fit length (0=dwt, 1=lenght, 2=length and no shrinking) (used in call_deri.m)
glo.mat   = 1;    % include maturity maint. (0=off, 1=include)
glo.moa   = 4;    % mode of action of toxicant

% syntax: par.name = [startvalue fit(0/1) minval maxval optional:log/normal scale (0/1)];
par.sJAm = [0.21  1 0 1e6 1]; % specific assimilation rate
par.sJM  = [0.10  1 0 1e6 1]; % specific maintenance costs
par.kap  = [0.41  1 0   1 1]; % allocation fraction to soma
par.WB0  = [8e-3  0 0 1e6 1]; % initial dry weight of egg
par.Lwp  = [1.6   1 0 1e6 1]; % body length at puberty
par.Lwf  = [0     0 0 1e6 1]; % body length at half-saturation of feeding (zero to ignore)
par.Lw0  = [0.97  0 0 1e6 1]; % body length at start experiment (not used as initial length is calculated from egg weight)

par.yAV  = [0.8   0 0 1 1];   % yield of assimilates on structure (starvation)
par.yBA  = [0.95  0 0 1 1];   % yield of egg buffer on assimilates
par.yVA  = [0.8   0 0 1 1];   % yield of structure on assimilates (growth)

par.ke   = [0.090  1 0  10 0]; % elimination rate constant, d-1
par.c0   = [0.030  1 0 1e6 1]; % no-effect concentration
par.cT   = [0.025  1 0 1e6 1]; % tolerance concentration
par.c0s  = [0.27   1 0 1e6 1]; % no-effect concentration survival
par.b    = [1.3    1 0 1e6 0]; % tolerance concentration survival
par.h0   = [0.0035 1 0 1e6 1]; % background hazard rate (d-1)
par.PRV  = [1      0 0 10  1]; % part. coeff. repro buffer and structure (weight basis)
% set to zero for no losses with repro, and 1 for same affinity structure and buffer

par.f    = [1   0 0  2 1]; % scaled functional response
par.Tlag = [0   0 0 30 1]; % lag time for start development

Time vector and labels for plots

Specify what to plot. If time vector glo.t is not specified, a default is used, based on the data set

% specify the y-axis labels for each state variable
glo.ylab{1} = ['scaled internal concentration (',char(181),'M)'];
glo.ylab{2} = 'body length (mm)';
glo.ylab{3} = 'cumul. repro. (predicted egg formation)';
glo.ylab{4} = 'survival fraction (-)';

% specify the x-axis label (same for all states)
glo.xlab    = 'time (days)';
glo.leglab1 = 'conc. '; % legend label before the 'scenario' number
glo.leglab2 = [char(181),'M']; % legend label after the 'scenario' number

prelim_checks % script to perform some preliminary checks and set things up
% Note: prelim_checks also fills all the options (opt_...) with defauls, so
% modify options after this call, if needed.

glo.eventson     = 1; % events function on (1) or off (0)

opt_optim.fit    = 1; % fit the parameters (1), or don't (0)
opt_optim.it     = 0; % show iterations of the optimisation (1, default) or not (0)
opt_plot.bw      = 0; % if set to 1, plots in black and white with different plot symbols
opt_plot.annot   = 2; % annotations in multiplot for fits: 1) box with parameter estimates 2) single legend
% opt_plot.statsup = [1 2]; % vector with states to suppress in plotting fits

% optimise and plot (fitted parameters in par_out)
par_out = calc_optim(par,opt_optim); % start the optimisation
calc_and_plot(par_out,opt_plot); % calculate model lines and plot them
Provisional goodness-of-fit measures for each data set
       NaN    0.9768    0.9763    0.9922

Warning: Goodness of fit measure for survival data is meaningless when more data
sets are fitted simultaneously! 
 
=================================================================================
Results of the parameter estimation with BYOM version 4.0
=================================================================================
   Filename      : byom_debkisstox_daphnia 
   Analysis date : 28-Apr-2017 (16:44) 
   Data entered  :
     data state 1: 0x0, no data.
     data state 2: 12x5, continuous data, no transf.
     data state 3: 10x5, continuous data, no transf.
     data state 4: 12x5, survival data.
   Search method: Nelder-Mead simplex direct search, 2 rounds. 
     The optimisation routine has converged to a solution
     Total 246 simplex iterations used to optimise. 
     Minus log-likelihood has reached the value 452.291 (AIC=924.582). 
=================================================================================
sJAm       0.2087 (fit: 1, initial: 0.21) 
sJM        0.1021 (fit: 1, initial: 0.1) 
kap        0.4125 (fit: 1, initial: 0.41) 
WB0         0.008 (fit: 0, initial: 0.008) 
Lwp         1.648 (fit: 1, initial: 1.6) 
Lwf             0 (fit: 0, initial:  0) 
Lw0          0.97 (fit: 0, initial: 0.97) 
yAV           0.8 (fit: 0, initial: 0.8) 
yBA          0.95 (fit: 0, initial: 0.95) 
yVA           0.8 (fit: 0, initial: 0.8) 
ke        0.08878 (fit: 1, initial: 0.09) 
c0        0.03005 (fit: 1, initial: 0.03) 
cT        0.02487 (fit: 1, initial: 0.025) 
c0s        0.2664 (fit: 1, initial: 0.27) 
b           1.309 (fit: 1, initial: 1.3) 
h0       0.003516 (fit: 1, initial: 0.0035) 
PRV             1 (fit: 0, initial:  1) 
f               1 (fit: 0, initial:  1) 
Tlag            0 (fit: 0, initial:  0) 
=================================================================================
Parameters ke and b are fitted on log-scale.
=================================================================================
Time required: 19.2 secs
Plots result from the optimised parameter values. 

Profiling the likelihood

By profiling you make robust confidence intervals for one or more of your parameters. Use the name of the parameter as it occurs in your parameter structure par above. You do not need to run the entire script before you can make a profile.

Options for the profiling can be set using opt_prof (see prelim_checks). For this demo, the level of detail of the profiling is set to 'course' and no sub-optimisations are used. However, consider these options when working on your own data (e.g., set opt_prof.subopt=10).

opt_prof.detail   = 2; % detailed (1) or a course (2) calculation
opt_prof.subopt   = 30; % number of sub-optimisations to perform to increase robustness

% % UNCOMMENT FOLLOWING LINE(S) TO CALCULATE
% % run a profile for selected parameters ...
% calc_proflik(par_out,'ke',opt_prof);  % calculate a profile
% calc_proflik(par_out,'c0',opt_prof);  % calculate a profile

% % UNCOMMENT FOLLOWING LINE(S) TO CALCULATE
% % or run profiles for all fitted parameters.
% all_profiles(par_out,opt_prof);

% % UNCOMMENT FOLLOWING LINE(S) TO CALCULATE
% % Automatically calculate profiles for all parameters, and redo
% % optimisation when a better value is found.
% opt_prof.subopt   = -1; % number of sub-optimisations to perform to increase robustness
% % Note: the -1 makes a comparison between 0 and 10 sub-optimisations.
% % In this example, sub-optimisations are helpful!
% par_out = auto_profiles(par_out,opt_prof,opt_optim); % Experimental!

Population growth rate

The function calc_pop allows for a calculation of the intrinsic rate of population increase (Euler-Lotka equation).

Tpop = linspace(0,42,100); % time vector for the population calculations
Cpop = linspace(0,1,100); % concentration range for population calculations
% Cpop = -1; % use only concentrations as given in X0mat
Th = 2.5; % time from fresh egg to t=0 in time vector (e.g., hatching time)
% Here, use 2.5 days, as this value was earlier subtracted from the time
% vector in the data set. So, reproduction now represents egg formation,
% and we need to account for hatching time in the population calculation.
trep = [0]; % vector with reproduction events (first one must be a zero)
% here, zero to base calculation on continuous reproduction
calc_pop(par_out,X0mat,[4,3],Tpop,Cpop,Th,trep)

Other files: derivatives

To archive analyses, publishing them with Matlab is convenient. To keep track of what was done, the file derivatives.m can be included in the published result.

%% 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.html
% 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.html 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: April 2017
% * Web support: <http://www.debtox.info/byom.html>
% * Back to index <walkthrough_debkisstox.html>

%% 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
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

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

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

Other files: call_deri

To archive analyses, publishing them with Matlab is convenient. To keep track of what was done, the file call_deri.m can be included in the published result.

%% BYOM function call_deri.m (calculates the model output)
%
%  Syntax: [Xout TE] = call_deri(t,par,X0v)
%
% This function calls the ODE solver to solve the system of differential
% equations specified in <derivatives.html derivatives.m>, or the explicit
% function(s) in <simplefun.html simplefun.m>. As input, it gets:
%
% * _t_   the time vector
% * _par_ the parameter structure
% * _X0v_   a vector with initial states and one concentration (scenario number)
%
% The output _Xout_ provides a matrix with time in rows, and states in
% columns. This function calls <derivatives.html derivatives.m>. The
% optional output _TE_ is the time at which an event takes place (specified
% using the events function). The events function is set up to catch
% discontinuities. It should be specified according to the problem you are
% simulating. If you want to use parameters that are (or influence) initial
% states, they have to be included in this function.
%
% * Author: Tjalling Jager
% * Date: April 2017
% * Web support: <http://www.debtox.info/byom.html>
% * Back to index <walkthrough_debkisstox.html>

%% Start

function [Xout TE] = call_deri(t,par,X0v)

global glo   % allow for global parameters in structure glo
global zvd   % global structure for zero-variate data

%% Initial settings
% This part extracts optional settings for the ODE solver that can be set
% in the main script (defaults are set in prelim_checks). The useode option
% decides whether to calculate the model results using the ODEs in
% <derivatives.html derivatives.m>, or the analytical solution in
% <simplefun.html simplefun.m>. Using eventson=1 turns on the events
% handling. Also modify the sub-function at the bottom of this function!
% Further in this section, initial values can be determined by a parameter
% (overwrite parts of X0), and zero-variate data can be calculated. See the
% example BYOM files for more information.

useode   = glo.useode; % calculate model using ODE solver (1) or analytical solution (0)
eventson = glo.eventson; % events function on (1) or off (0)
stiff    = glo.stiff; % set to 1 to use a stiff solver instead of the standard one

% Unpack the vector X0v, which is X0mat for one scenario
X0 = X0v(2:end); % these are the intitial states for a scenario

% % if needed, extract parameters from par that influence initial states in X0
% % start from specified initial size
% Lw0 = par.Lw0(1);  % initial body length (mm)
% X0(2) = glo.dV*((Lw0*glo.delM)^3); % recalculate to dwt, and put this parameter in the correct location of the initial vector

% % start from the value provided in X0mat
% if glo.len ~= 0 % only if the switch is set to 1 or 2! 
%     % Assume that X0mat includes initial length; if it includes mass, no correction is needed here
%     % The model is set up in masses, so we need to translate the initial length
%     % (in mm body length) to body mass using delM. Assume that size is second state.
%     X0(2) = glo.dV*((X0(2) * glo.delM)^3); % initial length to initial dwt
% end

% start at hatching, so calculate initial length from egg weight
WB0   = par.WB0(1);     % initial dry weight of egg
yVA   = par.yVA(1);     % yield of structure on assimilates (growth)
kap   = par.kap(1);     % allocation fraction to soma
WVb   = WB0 *yVA * kap; % dry body mass at birth
X0(2) = WVb;            % put this estimate in the correct location of the initial vector

% if needed, calculate model values for zero-variate data from parameter set
if ~isempty(zvd)
    zvd.Lwm(3) = (par.sJAm(1) * par.kap(1)/par.sJM(1))/glo.delM; % add model prediction as third value
end

%% Calculations
% This part calls the ODE solver (or the explicit model in <simplefun.html
% simplefun.m>) to calculate the output (the value of the state variables
% over time). There is generally no need to modify this part. The solver
% ode45 generally works well. For stiff problems, the solver might become
% very slow; you can try ode15s instead.

c  = X0v(1);     % the concentration (or scenario number)
t  = t(:);       % force t to be a row vector (needed when useode=0)

% specify options for the ODE solver
options = odeset; % start with default options
if eventson == 1
    options = odeset(options, 'Events',@eventsfun); % add an events function
end
options = odeset(options, 'RelTol',1e-4,'AbsTol',1e-7); % specify tightened tolerances
% options = odeset(options,'InitialStep',max(t)/1000,'MaxStep',max(t)/100); % specify smaller stepsize

TE = 0; % dummy for time of events

if useode == 1 % use the ODE solver to calculate the solution
    % call the ODE solver (use ode15s for stiff problems, escpecially for pulses)
    if isempty(options.Events) % if no events function is specified ...
        switch stiff
            case 0
                [~,Xout] = ode45(@derivatives,t,X0,options,par,c);
            case 1
                [~,Xout] = ode113(@derivatives,t,X0,options,par,c);
            case 2
                [~,Xout] = ode15s(@derivatives,t,X0,options,par,c);
        end
    else % with an events functions ... additional output arguments for events:
        % TE catches the time of an event, YE the states at the event, and IE the number of the event
        switch stiff
            case 0
                [~,Xout,TE,YE,IE] = ode45(@derivatives,t,X0,options,par,c);
            case 1
                [~,Xout,TE,YE,IE] = ode113(@derivatives,t,X0,options,par,c);
            case 2
                [~,Xout,TE,YE,IE] = ode15s(@derivatives,t,X0,options,par,c);
        end
    end
else % alternatively, use an explicit function provided in simplefun!
    Xout = simplefun(t,X0,par,c);
end

if isempty(TE) || all(TE == 0) % if there is no event caught
    TE = +inf; % return infinity
end

%% Output mapping
% _Xout_ contains a row for each state variable. It can be mapped to the
% data. If you need to transform the model values to match the data, do it
% here. 

% here, translate from the state body weight to output in length, if needed
if glo.len ~= 0 % only if the switch is set to 1 or 2!
    % Here, I translate the model output in body mass to estimated physical body length
    W = Xout(:,2); % Assume that dry body weight is the second state
    Lw = (W/glo.dV).^(1/3)/glo.delM; % estimated physical body length    
    if glo.len == 2 % when animal cannot shrink in length (but does on volume!)
        maxLw = 0;%zeros(size(Lw)); % remember max size achieved
        for i = 1:length(Lw) % run through length data
            maxLw = max([maxLw;Lw(i)]); % new max is max of old max and previous size
            Lw(i) = maxLw; % replace value in output vector
        end 
    end
    Xout(:,2) = Lw; % replace body weight by physical body length
end

% % To obtain the output of the derivatives at each time point. The values in
% % dXout might be used to replace values in Xout, if the data to be fitted
% % are the changes (rates) instead of the state variable itself.
% % dXout = zeros(size(Xout)); % initialise with zeros
% for i = 1:length(t) % run through all time points
%     dXout(i,:) = derivatives(t(i),Xout(i,:),par,c); 
%     % derivatives for each stage at each time
% end

%% Events function
% Modify this part of the code if _eventson_=1. This subfunction catches
% the 'events': in this case, this one catches when the scaled internal
% concentration exceeds one of the NECs, and catches the switch at puberty.
%
% Note that the eventsfun has the same inputs, in the same sequence, as
% <derivatives.html derivatives.m>.

function [value,isterminal,direction] = eventsfun(t,X,par,c)

global glo
Lwp = par.Lwp(1);  % length at puberty (mm)
c0  = par.c0(1);   % no-effect concentration
c0s = par.c0s(1);  % no-effect concentration survival

WVp = glo.dV * ((Lwp * glo.delM)^3); % translate physical length to dry weight

nevents = 3;         % number of events that we try to catch

value       = zeros(nevents,1); % initialise with zeros
value(1)    = X(1) - c0;        % follow when scaled internal conc exceed the NEC
value(2)    = X(1) - c0s;       % follow when scaled internal conc exceed the NEC survival
value(3)    = X(2) - WVp;       % follow when body weight exceeds length at puberty

isterminal  = zeros(nevents,1); % do NOT stop the solver at an event
direction   = zeros(nevents,1); % catch ALL zero crossing when function is increasing or decreasing