BYOM, simple-compound model: byom_debtox_daphnia.m

Table of contents

Contents

About

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: simple DEBtox model for toxicants, based on DEBkiss and formulated in compound parameters. The model includes flexible modules for toxicokinetics/damage dynamics and toxic effects. The DEBkiss e-book (see http://www.debtox.info/book_debkiss.html) provides a partial description of the model; a publication is in preparation that contains the full details.

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 due to several model changes.

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.

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 damage
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 per mother) (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
    2	0	0	0	0	0
    4	0	0	0	0	0
    6	0	0	0	0	0
    8	1.2     2.1     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];

% weights: number of individuals alive
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];

% In the next part, I subtract 3 days from the time vector for the repro
% data; reason is that the eggs are produced some 3 days before the
% offspring are counted. After this shift, the time vector thus represents
% egg formation, which is the relevant trait for DEB analysis.
DATA{3}(2:end,1) = DATA{3}(2:end,1) - 3; % extract x days from time vector
DATA{3}([2 3],:) = []; % remove first two time points from DATA (t<0)
W{3}([1 2],:)    = []; % remove first two time points from W (t<0)

% 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)
X0mat(2,:) = 0; % initial values state 1 (scaled damage)
X0mat(3,:) = 0; % initial values state 2 (body length, initial value overwritten by L0)
X0mat(4,:) = 0; % initial values state 3 (cumulative reproduction)
X0mat(5,:) = 1; % initial values state 4 (survival probability)

% put the position of the various states in globals, to make sure correct
% one is selected for extra things (e.g., in call_deri for accommodating
% 'no shrinking' and for population growth rate)
glo.locD = 1; % location of scaled damage in the state variable list
glo.locL = 2; % location of body size in the state variable list
glo.locR = 3; % location of cumulative reproduction in the state variable list
glo.locS = 4; % location of survival probability in the state variable list (use -1 if no survival state available/used)

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.FBV = 0.02;    % dry weight of egg as fraction of structural body weight (-) (for losses with reproduction)
glo.KRV = 1;       % part. coeff. repro buffer and structure (kg/kg) (for losses with reproduction)
glo.kap = 0.8;     % approximation for kappa (for starvation response)
glo.yP  = 0.8*0.8; % product of yVA and yAV (for starvation response)
glo.len = 1;       % set switch to 2 to prevent shrinking in length (used in call_deri.m)

% select mode of action of toxicant as set of switches:
% [assimilation/feeding, maintenance costs (somatic and maturity), growth costs, repro costs]
% glo.moa = [1 0 0 0]; % assimilation/feeding
% glo.moa = [0 1 0 0]; % costs for maintenance
% glo.moa = [0 0 1 1]; % costs for growth and reproduction
glo.moa = [0 0 0 1]; % costs for reproduction

% select which feedbacks to use on damage dynamics as set of switches:
% [surface:volume on uptake, surface:volume on elimination, growth dilution, losses with reproduction]
% glo.feedb = [1 1 1 1]; % all feedbacks
glo.feedb = [1 1 1 0]; % classic DEBtox (no losses with repro)
% glo.feedb = [0 0 1 0]; % damage that is diluted by growth
% glo.feedb = [0 0 0 0]; % damage that is not diluted by growth

% syntax: par.name = [startvalue fit(0/1) minval maxval optional:log/normal scale (0/1)];
par.L0   = [0.857   1 0 1e6 1]; % body length at start experiment (mm)
par.Lp   = [1.78    1 0 1e6 1]; % body length at puberty (mm)
par.Lm   = [3.09    1 0 1e6 1]; % maximum body length (mm)
par.rB   = [0.150   1 0 1e6 1]; % von Bertalanffy growth rate constant (1/d)
par.Rm   = [9.40    1 0 1e6 1]; % maximum reproduction rate (#/d)
par.f    = [1       0 0   2 1]; % scaled functional response
par.hb   = [0.00481 1 0 1e6 1]; % background hazard rate (d-1)
% Note: it does not matter whether length measures are entered as actual
% length or as volumetric length (as long as the same measure is used consistently).

% extra parameters for special situations
par.Lf   = [0 0 0 1e6 1]; % body length at half-saturation feeding (mm)
par.Tlag = [0 0 0  30 1]; % lag time for start development

ind_tox = length(fieldnames(par))+1; % index where tox parameters start
% the parameters below this line are all treated as toxicity parameters!
par.kd   = [0.1  1 0  10 0]; % dominant rate constant (d-1)
par.zb   = [0.1  1 0 1e6 1]; % effect threshold energy budget ([C])
par.bb   = [100  1 0 1e6 0]; % effect strength energy-budget effect (1/[C])
par.zs   = [0.2  1 0 1e6 1]; % effect threshold survival ([C])
par.bs   = [0.6  1 0 1e6 0]; % effect strength survival (1/([C] d))

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 damage (',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.

Calculations and plotting

Here, the function is called that will do the calculation and the plotting. Options for the plotting can be set using opt_plot (see prelim_checks.m). Options for the optimsation routine can be set using opt_optim. Options for the ODE solver are part of the global glo.

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]; % vector with states to suppress in plotting fits

% Select whether to fit the control treatment (0), the toxicity treatments
% (1; the control is shown as well), or both together (2). A good strategy
% is to fit the basic parameters to the control data first (fit_tox=0),
% copy the best values into the parameter matrix above, and then fit the
% toxicity parameter to the complete data set (fit_tox=1). The code below
% automatically keeps the parameters fixed that need to be fixed.

fit_tox = 1; % set by user

if fit_tox == 0 % for fitting the control response only ...
    % This piece of code selects only the treatment with identifier 0
    % (zero), and turns fitting off for the toxicity parameters.
    X0mat = X0mat(:,X0mat(1,:)==0); % only keep controls
    pmat = packunpack(1,par,[]); % turn parameter structure into a matrix
    pmat(ind_tox:end,2) = 0; % don't fit the toxicity parameters now
    par = packunpack(2,[],pmat); % turn parameter matrix back into a structure
    glo2.ctot = 0; % make sure that the total concentration vector to analyse only has a zero
elseif fit_tox == 1 % for fitting tox parameters only ...
    % This piece of code turns fitting off for the basic parameters so only
    % the toxicity parameters are fitted.
    pmat = packunpack(1,par,[]); % turn output parameter structure into a matrix
    pmat(1:ind_tox-1,2) = 0; % turn fitting for the basic parameters off in this vector
    par = packunpack(2,[],pmat); % turn parameter matrix into a structure
end

% 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

% print_par(par_out) % this prints out the optimised parameter values in a
% % formatted way so they can be directly copied into the code above
Goodness-of-fit measures for each data set (R-square)
       NaN    0.9743    0.9870    0.9148

Warning: R-square is not appropriate for survival data and needs to be
interpreted more qualitatively! 
 
=================================================================================
Results of the parameter estimation with BYOM version 4.3
=================================================================================
   Filename      : byom_debtox_daphnia_tox 
   Analysis date : 02-Sep-2019 (17:15) 
   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 209 simplex iterations used to optimise. 
     Minus log-likelihood has reached the value 439.83563 (AIC=889.671). 
=================================================================================
L0          0.857 (fit: 0, initial: 0.857) 
Lp           1.78 (fit: 0, initial: 1.78) 
Lm           3.09 (fit: 0, initial: 3.09) 
rB           0.15 (fit: 0, initial: 0.15) 
Rm            9.4 (fit: 0, initial: 9.4) 
f               1 (fit: 0, initial:  1) 
hb        0.00481 (fit: 0, initial: 0.00481) 
Lf              0 (fit: 0, initial:  0) 
Tlag            0 (fit: 0, initial:  0) 
kd        0.07658 (fit: 1, initial: 0.1) 
zb        0.09632 (fit: 1, initial: 0.1) 
bb          74.25 (fit: 1, initial: 100) 
zs         0.2321 (fit: 1, initial: 0.2) 
bs         0.6606 (fit: 1, initial: 0.6) 
=================================================================================
Parameters kd, bb and bs are fitted on log-scale.
=================================================================================
Mode of action used: 0  0  0  1
Feedbacks used     : 1  1  1  0
=================================================================================
Time required: 30.4 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 coarse (2) calculation
opt_prof.subopt   = 0; % 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,'kd',opt_prof);  % calculate a profile
% calc_proflik(par_out,'zb',opt_prof);  % calculate a profile

% % UNCOMMENT FOLLOWING LINE(S) TO CALCULATE
% % Automatically calculate profiles for all parameters, and redo
% % optimisation when a better value is found.
% opt_prof.subopt   = 0; % number of sub-optimisations to perform to increase robustness
% % Note: a -1 makes a comparison between 0 and 10 sub-optimisations.
% 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). The function calc_pop needs checking ... so use with care! By default, the calculation is based on continuous reproduction (like the fit).

% Population calculations can be performed with CIs
opt_conf.type    = 0; % use the values from the slice sampler (1) or likelihood region (2) to make intervals
opt_conf.lim_set = 0; % for lik-region sample: use limited set of n_lim points (1) or outer hull (2) to create CIs

Tpop = linspace(0,50,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 = 3; % time from fresh egg to t=0 in time vector (e.g., hatching time)
% Here, use 3 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.

calc_pop(par_out,X0mat,Tpop,Cpop,Th,opt_conf,opt_pop)

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

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

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

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

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>. It is linked to the model in
% the directory 'simple_compound'. 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: September 2019
% * Web support: <http://www.debtox.info/byom.html>
% * Back to index <walkthrough_debtox2019.html>
% 
%  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. 

%% 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 in a model parameter
L0 = par.L0(1); % initial body length (mm) is a parameter
X0(glo.locL) = L0; % put this estimate in the correct location of the initial vector

%% 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',(t(end)-t(1))/1000,'MaxStep',(t(end)-t(1))/100); % specify smaller stepsize
% This small time step will make the simulations considerably slower.
% However, this extra precision is needed for time-varying concentrations,
% time-varying food levels, and possibly also for starvation. You can
% comment out this option setting for extra speed, if you know what you're
% doing.

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, especially 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. 

if glo.len == 2 % when animal cannot shrink in length (but does on weight!)
    L = Xout(:,glo.locL); % take correct state for body lenght
    maxL = 0; % initialise value for max size achieved
    for i = 1:length(L) % run through length data
        maxL = max([maxL;L(i)]); % new max is max of old max and previous size
        L(i) = maxL; % replace value in output vector
    end
    Xout(:,glo.locL) = L; % 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
Lp = par.Lp(1); % length at puberty (mm)
zb = par.zb(1); % effect threshold for the energy budget
zs = par.zs(1); % effect threshold for survival

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

value    = zeros(nevents,1); % initialise with zeros
value(1) = X(glo.locD) - zb; % follow when scaled damage exceeds the effect threshold for the energy budget
value(2) = X(glo.locD) - zs; % follow when scaled damage exceed the effect threshold for survival
value(3) = X(glo.locL) - Lp; % follow when body length 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