BYOM, byom_debtox_daphnia.m, DEBtox with Daphnia

Author: Tjalling Jager
Date: November 2021
Web support: http://www.debtox.info/byom.html
Back to index walkthrough_debtox.html

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: DEBtox following Jager & Zimmer (2012), http://dx.doi.org/10.1016/j.ecolmodel.2011.11.012.

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

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.

Initial things
The data set
Initial values for the state variables
Initial values for the model parameters
Time vector and labels for plots
Calculations and plotting
Profiling the likelihood
Population growth rate
Other files: derivatives
Other files: call_deri

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
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(0) % set path to the BYOM/engine directory (option 1 uses parallel toolbox)
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:

-1 for multinomial likelihood (for survival data)
0 for log-transform the data, then normal likelihood
0.5 for square-root transform the data, then normal likelihood
1 for no transformation of the data, then normal likelihood

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

% scaled reserve density
DATA{2} = [0]; % there are never data for this state

% body length (mm) on each observation time (d), concentrations in uM
DATA{3} = [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{3} = 5 * ones(size(DATA{3})-1); % weights: 5 animals in each treatment

% cumulative reproduction (nr. offspring) (mm) on each observation time (d), concentrations in uM
DATA{4} = [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{4} = [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

glo.Tbp = 2.5; % shift model output for brood-pouch delay (rather than change the data set)

% survivors on each observation time (d), concentrations in uM
DATA{5} = [-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];

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
         1 1   1   1  % initial values state 2
         0 0   0   0  % initial values state 3 (overwritten by L0)
         0 0   0   0  % initial values state 4
         1 1   1   1];% initial values state 5

% Put the position of the various states in globals, to make sure correct
% one is selected for extra things (e.g., using plot_tktd for nicer plots);
% note that this must match the positions in the state vector in call_deri
% and derivatives.
glo.locD = 1; % location of scaled conc./damage in the state variable list
glo.locL = 3; % location of body size in the state variable list
glo.locR = 4; % location of cumulative reproduction in the state variable list
glo.locS = 5; % location of survival probability in the state variable list

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.moa = 4; % mode of action (see derivatives)

% syntax: par.name = [startvalue fit(0/1) minval maxval optional:log/normal scale (0/1)];
par.ke  = [0.039  1 0.01 10  0]; % elimination rate constant (d-1)
par.c0  = [0.063  1 0    1e6 1]; % no-effect concentration sub-lethal (uM)
par.cT  = [7e-3   1 0    1e6 1]; % tolerance concentration (uM)
par.c0s = [0.15   1 0    1e6 1]; % no-effect concentration survival (uM)
par.b   = [1.1    1 1e-6 1e6 0]; % killing rate (1/(d uM))

par.L0  = [0.88   0 1e-3 100 1]; % initial body length (mm)
par.Lp  = [1.6    1 1e-3 100 1]; % length at puberty (mm)
par.Lm  = [3.1    1 1e-3 100 1]; % maximum length (mm)
par.rB  = [0.14   1 0 100 1];    % von Bertalanffy growth rate (d-1)
par.Rm  = [11     1 0 1e6 1];    % maximum reproduction rate (eggs/d)
par.g   = [0.422  0 0 100 1];    % energy-investment ratio (-)
par.f   = [1      0 0 2   1];    % scaled functional response (-)
par.h0  = [0.0024 1 0 1e6 1];    % background hazard rate (d-1)

Time vector and labels for plots

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

% specify the y-axis labels for each state variable
glo.ylab{1} = ['scaled int. conc. (',char(181),'M)'];
glo.ylab{2} = 'scaled reserve dens. (-)';
glo.ylab{3} = 'body length (mm)';
if isfield(glo,'Tbp') && glo.Tbp > 0
    glo.ylab{4} = ['cumul. repro. (shift ',num2str(glo.Tbp),'d)'];
else
    glo.ylab{4} = 'cumul. repro. (no shift)';
end
glo.ylab{5} = '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.

NOTE: for this package, the options useode and eventson in glo will not be functional: the ODE solver is always used, and the events function as well.

glo.stiff = [0 3]; % ODE solver 0) ode45 (standard), 1) ode113 (moderately stiff), 2) ode15s (stiff)
% Second argument is for normally tight (1), tighter (2), or very tight (3)
% tolerances. Use 1 for quick analyses, but check with 3 to see if there is
% a difference! Especially for time-varying exposure, there can be large
% differences between the settings!
glo.break_time = 0; % break time vector up for ODE solver (1) or don't (0)
% Note: breaking the time vector is a good idea when the exposure scenario
% contains discontinuities. Don't use for continuous splines (type 1) as it
% will be much slower. For FOCUS scenarios (high time resolution), breaking
% up is not efficient and does not appear to be necessary.

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

% print_par(par_out,-2) % this prints out the optimised parameter values in a
% % formatted way so they can be directly copied into the code above; the
% % -2 means that only fitted parameters are displayed

% these plots are more readable, especially when plotting CIs
opt_tktd.repls = 0; % plot individual replicates (1) or means (0)
plot_tktd(par_out,opt_tktd,[]);
% Leave the options opt_conf empty to suppress all CIs for these plots.
% Plotting CIs requires a sample as saved by the various methods available
% in BYOM (see examples directory).

Goodness-of-fit measures for each state and each data set (R-square)
Based on the individual replicates, accounting for transformations, and including t=0
state 1.    NaN
state 2.    NaN
state 3.   0.98
state 4.   0.99
state 5.   0.96
Warning: R-square is not appropriate for survival data so interpret
qualitatively! 
Warning: (R-square not calculated when missing/removed animals; use plot_tktd) 
 
=================================================================================
Results of the parameter estimation with BYOM version 6.0 BETA 7
=================================================================================
   Filename      : byom_debtox_daphnia 
   Analysis date : 30-Nov-2021 (11:39) 
   Data entered  :
     data state 1: 0x0, no data.
     data state 2: 0x0, no data.
     data state 3: 12x5, continuous data, no transf.
     data state 4: 12x5, continuous data, no transf.
     data state 5: 12x5, survival data.
   Search method: Nelder-Mead simplex direct search, 2 rounds. 
     The optimisation routine has converged to a solution
     Total 1075 simplex iterations used to optimise. 
     Minus log-likelihood has reached the value 498.76 (AIC=1017.51). 
=================================================================================
ke        0.06311 (fit: 1, initial: 0.039) 
c0        0.08687 (fit: 1, initial: 0.063) 
cT        0.01170 (fit: 1, initial: 0.007) 
c0s        0.2044 (fit: 1, initial: 0.15) 
b          0.7688 (fit: 1, initial: 1.1) 
L0           0.88 (fit: 0, initial: 0.88) 
Lp          1.607 (fit: 1, initial: 1.6) 
Lm          3.101 (fit: 1, initial: 3.1) 
rB         0.1402 (fit: 1, initial: 0.14) 
Rm          10.60 (fit: 1, initial: 11) 
g           0.422 (fit: 0, initial: 0.422) 
f               1 (fit: 0, initial:  1) 
h0       0.002363 (fit: 1, initial: 0.0024) 
=================================================================================
Parameters ke and b are fitted on log-scale.
=================================================================================
Mode of action used: 4
=================================================================================
Time required: 1 min, 13.6 secs
Plots result from the optimised parameter values. 
 
Plot_tktd uses uses parameter set from user input for best-fit curves.
No confidence intervals are produced as opt_conf is empty
 
Goodness-of-fit metrics (treat with caution!)
 Values below are calculated on means, incl. t=0, from the residuals
 Residuals are used without accounting for any transformations
 However, means are calculated with transformation (if specified in DATA)
 
            Model efficiency (r2)   NRMSE
================================================================================
  survival    : 0.9736              0.0355 
  body length : 0.9796              0.0406 
  reproduction: 0.9863              0.1872 
================================================================================

Profiling the likelihood

By profiling you make robust confidence intervals for one or more of your parameters. Use the names of the parameters as they occurs in your parameter structure par above. This can be a single string (e.g., 'kd'), a cell array of strings (e.g., {'ke','c0'}), or 'all' to profile all fitted parameters. If you have used the parameter-space explorer for fitting, you already have (even more robust) profiles!

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 'coarse' 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
opt_prof.brkprof  = 2; % when a better optimum is located, stop (1) or automatically refit (2)

calc_proflik(par_out,{'ke','c0'},opt_prof,opt_optim);  % calculate a profile
% Enter single parameter names, a cell array of names, or 'all' to profile
% all fitted parameters.

  
95% confidence interval(s) from the profile(s)
=================================================================================
ke         interval:    0.01047 -  0.1437 
c0         interval:    0.01995 -  0.1351 
=================================================================================
 
Time required: 28 mins, 29.8 secs

Population growth rate

The function calc_pop allows for a calculation of the intrinsic rate of population increase (Euler-Lotka equation). It calculates the growth rate at three food levels (set in the engine file calc_pop.m).

% Population calculations can be performed with CIs
opt_conf.type    = 0; % make intervals from 1) slice sampler, 2)likelihood region, 3) parspace explorer
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,42,100); % time vector for the population calculations
Cpop = linspace(0,1,50);   % concentration range for population calculations
% Cpop = -1; % use only concentrations as given in X0mat
Th = glo.Tbp; % time from fresh egg to t=0 in time vector (e.g., hatching time)
% Here, use glo.Tbp 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,glo)
%
% This function calculates the derivatives for the model system. It is
% linked to the script file <byom_debtox_daphnia.html
% 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)
% * _glo_ is the structure with information (normally global)
%
% 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 2020
% * Web support: <http://www.debtox.info/byom.html>
% * Back to index <walkthrough_debtox.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. 

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

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

%% 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) (only used in call_deri.m)
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

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

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

% OPTION 1: DEBtox default, growth dilution and surface:volume ratio
dcV = ke*(Lm/L)*(c-cV) - cV*(3/L)*dL; % change in scaled internal conc./damage
% % OPTION 2: growth dilution only
% dcV = ke*(c-cV) - cV*(3/L)*dL; % change in scaled internal conc./damage
% % OPTION 3: no growth dilution and no surface:volume ratio
% dcV = ke*(c-cV); % change in scaled internal conc./damage

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

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,Xout2,zvd] = call_deri(t,par,X0v,glo)
%
% This function calls the ODE solver to solve the system of differential
% equations specified in <derivatives.html derivatives.m>. It is specific
% for the model in the directory 'DEBtox_classic', and modified to deal
% with time-varying exposure by splitting up the ODE solving to avoid hard
% switches.
%
% As input, it gets:
% * _t_   the time vector
% * _par_ the parameter structure
% * _X0v_   a vector with initial states and one concentration (scenario number)
% * _glo_ the structure with various types of information (used to be global)
%
% 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: November 2021
% * Web support: <http://www.debtox.info/byom.html>
% * Back to index <walkthrough_debtox.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. 

%% Start

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

% These outputs need to be defined, even if they are not used
Xout2    = []; % additional uni-variate output, not used in this example
zvd      = []; % additional zero-variate output, not used in this example

%% 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). 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.
% 
% This package will always use the ODE solver; therefore the general
% settings in the global glo for the calculation (useode and eventson) are
% removed here. This packacge also always uses the events function, so the
% option eventson is not used. Furthermore, simplefun.m is removed.

stiff      = glo.stiff;      % ODE solver 0) ode45 (standard), 1) ode113 (moderately stiff), 2) ode15s (stiff)
break_time = glo.break_time; % break time vector up for ODE solver (1) or don't (0)
min_t      = 500; % minimum length of time vector (affects ODE stepsize only, when needed)

if length(stiff) == 1
    stiff(2) = 1; % by default: normally tightened tolerances
end

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

% Start with a check on time-varying exposure scenarios. This is a
% time-saver! The isfield and ismember calls take some time that rapidly
% multiplies as derivatives is called many times.
glo.timevar = [0 0]; % flag for time-varying exposure (second element is to tell read_scen which interval we need)
if isfield(glo,'int_scen') && ismember(c,glo.int_scen) % is c in the scenario range of the global?
    glo.timevar = [1 0]; % tell derivatives that we have a time-varying treatment
end

% if needed, extract parameters from par that influence initial states in X0
L0    = par.L0(1);  % initial body length (mm)
X0(3) = L0;         % put this parameter in the correct location of the initial vector

% % if needed, calculate model values for zero-variate data from parameter set
% if ~isempty(zvd)
%     zvd.ku(3) = par.Piw(1) * par.ke(1); % add model prediction as third value
% end

%% Calculations
% This part calls the ODE solver 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.

t  = t(:);      % force t to be a row vector (needed when useode=0)
t_rem = t;      % remember the original time vector (as we will add to it)

% The code below is meant for discontinous time-varying exposure. It will
% also work for constant exposure. Breaking up will not be as effective for
% piecewise polynomials (type 1), and slow, so better turn that off
% manually when you try splining. The code will run each interval between
% two exposure 'events' (in Tev) separately, stop the solver, and restart.
% This ensures that the discontinuities are no problem anymore.

TE  = 0; % dummy for time of events in the events function
Tev = [0 c]; % exposure profile events setting: without anything else, assume it is constant
glo.timevar = [0 0]; % flag for time-varying exposure (second element is to tell read_scen which interval we need)

InitialStep = max(t)/100; % specify initial stepsize
MaxStep     = max(t)/10;  % specify maximum stepsize
% For constant concentrations, and when breaking the time vector, we can
% use a default; Matlab uses as default the length of the time vector
% divided by 10.

if glo.timevar(1) == 1 % time-varying concentrations?
    % if it exists: use it to derive current external conc.
    Tev = read_scen(-2,c,t,glo); % use read_scen to derive actual exposure concentration
    % the -2 lets read_scen know we need events
    
    if break_time == 0
        min_t = max(min_t,length(Tev(Tev(:,1)<t(end),1))*2);
        % The min_t is only used to limit the step size of the ODE
        % solver, since the all calculations are done in derivatives.m.
        InitialStep = t(end)/(10*min_t); % initial step size
        MaxStep     = t(end)/min_t;      % maximum step size
        % For the ODE solver, when we have a time-varying exposure set
        % here, we base minimum step size on min_t. Small stepsize is a
        % good idea for pulsed exposures; otherwise, stepsize may become so
        % large that a concentration change is missed completely. When we
        % break the time vector, limiting step size is not needed.
    end
end

% This is a means to include a delay caused by the brood pounch in species
% like Daphnia. The repro data are for the appearance of neonates, but egg
% production occurs earlier. This global shifts the model output in this
% function below. This way, the time vector in the data does not need to be
% manipulated, and the model plots show the neonate production as expected.
% (NEW)
bp = 0; % by default, no brood-pouch delay
if isfield(glo,'Tbp') && glo.Tbp > 0 % if there is a global specifying a brood-pouch delay ...
    tbp = t(t>glo.Tbp)-glo.Tbp; % extra times needed to calculate brood-pounch delay
    t   = unique([t;tbp]);      % add the shifted time points to account for brood-pouch delay
    bp  = 1;                    % signal rest of code that we need brood-pouch delay
end

% specify options for the ODE solver
options = odeset; % start with default options for the ODE solver
% This needs further study ... events function removed. Events function
% needs to be considered very carefully for this model (and would only be
% useful for SD).
switch stiff(2)
    case 1 % for ODE15s, slightly tighter tolerances seem to suffice (for ODE113: not tested yet!)
        RelTol  = 1e-4; % relative tolerance (tightened)
        AbsTol  = 1e-7; % absolute tolerance (tightened)
    case 2 % somewhat tighter tolerances ...
        RelTol  = 1e-5; % relative tolerance (tightened)
        AbsTol  = 1e-8; % absolute tolerance (tightened)
    case 3 % for ODE45, very tight tolerances seem to be necessary in some cases
        RelTol  = 1e-9; % relative tolerance (tightened)
        AbsTol  = 1e-9; % absolute tolerance (tightened)
end
options = odeset(options,'RelTol',RelTol,'AbsTol',AbsTol,'Events',@eventsfun,'InitialStep',InitialStep,'MaxStep',MaxStep); % set options
% Note: setting tolerances is pretty tricky. For some cases, tighter
% tolerances are needed but not for others. For ODE45, tighter tolerances
% seem to work well, but not for ODE15s.

T = Tev(:,1); % time vector with events
if T(end) > t(end) % scenario may be longer than t(end)
    T(T>t(end)) = []; % remove all entries that are beyond the last time point
    % this may remove one point too many, but that will be added next
end
if T(end) < t(end) % scenario may (now) be shorter than we need
    T = cat(1,T,t(end)); % then add last point from t
end
% If a lag time is used, it is a good idea to add it as an event as ODE45
% does not like such a switch either.
if isfield(par,'Tlag') && par.Tlag(1) > 0
    T = unique([par.Tlag(1);T]);
end

% Breaking up the time vector is faster and more robust when calibrating on
% data from time-varying exposure that contain discontinuities. However,
% when running through very detailed exposure profiles (e.g., those from
% FOCUS), it is probably best to use break_time=0. Testing indicates that
% breaking up will be much slower than simply running the ODE solver across
% the entire profile (and produces very similar output).
% 
% Note: it is possible to initialise tout and Xout with NaNs, and avoid
% these matrices to grow. This may be a bit slower for exposure scenarios
% with few events, but faster for scenarios with many.

t = unique([T;t;(T(1:end-1)+T(2:end))/2]); % combine T, t, and halfway-T into new time vector
% this hopefully prevents the ODE solver from missing exposure pulses, and
% it makes sure that, for break_time=1, we never have a temporary time
% vector of length 2.

if break_time == 0 
    
    % simply use the ODE solver for the entire time vector
    switch stiff(1)
        case 0
            [tout,Xout,TE,~,~] = ode45(@derivatives,t,X0,options,par,c,glo);
        case 1
            [tout,Xout,TE,~,~] = ode113(@derivatives,t,X0,options,par,c,glo);
        case 2
            [tout,Xout,TE,~,~] = ode15s(@derivatives,t,X0,options,par,c,glo);
    end
    
else
    
    % Transferring the interval to derivatives is a huge time saver! It is
    % not safe to use i as index to Tev since we may add Tlag to T (which
    % is not in Tev); so better calculate it explicitly. Using <find> may
    % be a bit slow, that is why this is done outside of the loop calling
    % the ODE solver.
    ind_Tev = zeros(length(T)-1,1);
    for i = 1:length(T)-1 % run through all intervals between events
        ind_Tev(i) = find(Tev(:,1) <= T(i),1,'last');
    end
    
    % NOTE: predefining *should* increase speed, but the speed gain is not
    % so clear in this case. I expect some gain when using many exposure
    % intervals.
    tout    = nan(length(t)+5,1); % this vector will collect the time output from the solver (5 more than needed)
    Xout    = nan(length(tout),length(X0)); % this matrix will collect the states output from the solver (5 more than needed)
    ind_i   = 1; % index for where we are in tout and Xout
    % Note: I initialise tout and Xout with more than the elements of t
    % because the events may be added as well (only stopping events, which
    % are not used yet). The number is rather arbitrary but should be equal
    % to, or more than, the number of events.

    % run the ODE solver piece-wise across all exposure events specified
    for i = 1:length(T)-1 % run through all intervals between events
        t_tmp = [T(i);T(i+1)]; % start with start and end time for this period
        t_tmp = unique([t_tmp;t(t<t_tmp(2) & t>t_tmp(1))]); % and add the time points from t that fit in there
        glo.timevar(2) = ind_Tev(i); % tell derivatives which part of Tev we're in
        
        % NOTE: adding a time point in between when t_tmp is length 2 is
        % not needed. We've added T and half-way-into T into time vector t.
        % Since we do NOT stop at events, there is no way to have a time
        % vector of two elements!

        % use ODE solver to find solution in this interval
        % Note: the switch-case may be compromising calculation speed; however,
        % it needs to be investigated which solver performs best
        switch stiff(1)
            case 0
                [tout_tmp,Xout_tmp,TE,~,~] = ode45(@derivatives,t_tmp,X0,options,par,c,glo);
            case 1
                [tout_tmp,Xout_tmp,TE,~,~] = ode113(@derivatives,t_tmp,X0,options,par,c,glo);
            case 2
                [tout_tmp,Xout_tmp,TE,~,~] = ode15s(@derivatives,t_tmp,X0,options,par,c,glo);
        end
        
        % collect output in correct location
        nt = length(tout_tmp); % length of the output time vector
        tout(ind_i:ind_i-1+nt)   = tout_tmp; % add tout_tmp to correct position in tout
        Xout(ind_i:ind_i-1+nt,:) = Xout_tmp; % add Xout_tmp to correct position in Xout

        ind_i = ind_i-1+nt;       % update ind_i for next round
        X0    = Xout_tmp(end,:)'; % update X0 for next round
        % Note: the way ind_i is used, there will be no double time points
        % in tout and Xout.
    end
    
end

if isempty(TE) || all(TE == 0) % if there is no event caught
    TE = +inf; % return infinity
end
% This is not so useful as only the TE found in the last time period is
% returned.

%% 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 bp == 1 % if we need a brood-pouch delay ... (NEW)
    [~,loct] = ismember(tbp,tout); % find where the extra brood-pouch time points are in the long Xout
    Xbp      = Xout(loct,glo.locR); % only keep the brood-pouch ones we asked for
end
t = t_rem; % return t to the original input vector

% Select the correct time points to return to the calling function
[~,loct] = ismember(t,tout); % find where the requested time points are in the long Xout
Xout     = Xout(loct,:);     % only keep the ones we asked for

if bp == 1 % if we need a brood-pouch delay ... (NEW)
    [~,loct] = ismember(tbp+glo.Tbp,t); % find where the extra brood-pouch time points SHOULD BE in the long Xout
    Xout(:,glo.locR)    = 0;   % clear the reproduction state variable
    Xout(loct,glo.locR) = Xbp; % put in the brood-pouch ones we asked for
end

% Xout(:,1) = Xout(:,1).^3; % e.g., do something on first column, like cube it ...

% % 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,glo); 
%     % derivatives for each stage at each time
% end

%% Events function
% This subfunction catches the 'events': in this case, this one catches
% when the scaled internal concentration/damage exceeds one of the
% thresholds, 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,glo)

% Note: glo needs to be an input to the events function as well, rather
% than a global since the inputs for the events function must match those
% for derivatives.

Lp  = par.Lp(1);  % length at puberty (mm)
c0  = par.c0(1);  % no-effect concentration sub-lethal
c0s = par.c0s(1); % no-effect concentration survival

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

value       = zeros(nevents,1); % initialise with zeros
value(1)    = X(3) - Lp;        % follow when body length exceeds length at puberty
value(2)    = X(1) - c0;        % follow when scaled conc./damage exceed the NEC
value(3)    = X(1) - c0s;       % follow when scaled conc./damage exceed the NEC
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