BYOM, byom_bioconc_extra.m, a detailed example

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. Fitting relies on the multinomial likelihood for survival and independent normal distributions (if needed after transformation) for continuous data. Results are shown on screen but also saved to a log file (results.out).

In this example file (byom_bioconc_extra), we will walk step-by-step through a byom script file, explaining what happens and showing the outputs.

The model: An organism is exposed to a chemical in its surrounding medium. The animal accumulates the chemical according to standard one-compartment first-order kinetics, specified by a bioconcentration factor (Piw) and an elimination rate (ke). The chemical degrades at a certain rate (kd). When the external concentration reaches a certain concentration (Ct), degradation stops. This is useful to demonstrate the events function in call_deri.m, which catches this discontuity graciously. In the form of ODE's:

$$ \frac{d}{dt}C_w=-k_d C_w \quad \textrm{as long as } C_w>C_t $$

$$ \frac{d}{dt}C_i=k_e(P_{iw}C_w-C_i) $$

This script: byom_bioconc_extra demonstrates fitting using replicated data, with different concentration vectors in both data sets. Furthermore, several options are explained, such as the inclusion of zero-variate data and priors for Bayesian analysis. More support in the BYOM manual downloadable from here. A stripped version (less text, less options) can be found as byom_bioconc_start.m.

Contents

Initial things

Before we start, the memory is cleared, globals are defined, and the diary is turned off (output will be collected in the file results.out). The function pathdefine.m makes sure that the engine directory is added to the path. 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:

For each state variable there should be a data set in the correct position. So, DATA{1} should be the data for state variable 1, etc. The curly braces indicate that this is a 'cell array'. If there are no data for a state, use DATA{i}=0 (if you forget, this will be automatically done in the engine script prelim_checks.m). Use NaN for missing data.

Multiple data sets for each state can now be used. The first data set for state 1 becomes DATA{1,1} and the second DATA{2,1}, etc. Note that each data set is treated as completely independent. For continuous data that means that each data set has its own error distribution (either treated as 'nuisance parameter' or provided by the user in glo.var). You can enter replicated data by adding columns with the same scenario number. The scenario number should occur only once in the initial values matrix X0mat.

Note: for using survival data: enter the survival data as numbers of survivors, and not as survival probability. The model should be set up to calculate probabilities though. Therefore, the state variable is a probability, and hence the initial value in X0mat below should be a probability too (generally 1). This deluxe script automatically translates the data into probabilities for plotting the results. For survival data, the weights matrix has a different meaning: it is used to specify the number of animals that went missing or were removed during the experiment (enter the number of missing/removed animals at the time they were last seen alive in the test).

Note: In this example, the time vectors are not the same for both data sets. This is no problem for the analysis. Also the concentration vectors differ, which is also accounted for, but take care to construct the matrix with initial states carefully to catch ALL scenarios.

% observed external concentrations (two series per concentration)
DATA{1} = [ 1  10 10 20 20 30 30
            0  11 12 20 19 28 29
            10  6  8 13 12 20 18
            20  4  5  9  7 10 12
            30  5  7  6  4  6  7
            40  4  5  6  5  4  3];

% weight factors (number of replicates per observation in DATA{1})
W{1} = [10 10 10 10 10 10
        10 10 10 10 10 10
         8  8  8  8  8  9
         8  7  7  8  8  9
         6  6  6  6  6  7];

% observed internal concentrations (two series per concentration)
DATA{2} = [0.5   10  10   20   20   50   50
           0    0   0    0    0    0    0
           10 570 600 1150 1100 2500 3000
           20 610 650 1300 1210 2700 2800
           25 590 620  960  900 2300 2200
           40 580 560  920  650 1200 1700];

% if weight factors are not specified, ones are used for continuous data,
% and zero's for survival data.

Initial values for the state variables

For each state variable, we need to specify initial values for each scenario that we want to fit or simulate in the matrix X0mat. The first row specifies the scenarios (here: exposure treatments) that we want to model. Here, there are 4 scenarios in Xmat, but each data set has only 3 of them. Watch the plot to see how BYOM deals with that situation. If you do not want to fit certain scenarios (exposure treatments) from the data, simply leave them out of X0mat.

Note: if you do not want to start at t=0, specify the exact time in the global variable glo.Tinit here (e.g., glo.Tinit = 100;). If it is not specified, zero is used. The X0mat thus defines the states at glo.Tinit, and not necessarily the value at the first data point. Plotting also always starts from glo.Tinit.

Initial states, scenarios in columns, states in rows. First row are the identifiers of all scenarios (here: nominal concentrations). Second row is external concentration in mg/L, and third row body residues. For replicated data, scenarios should occur only once in X0mat.

X0mat = [10 20 30 50    % the scenarios (here nominal concentrations)
          9 18 27 47    % initial values state 1 (actual external concentrations)
          0  0  0  0];  % initial values state 2 (internal concentrations)

% Note: The initial external concentrations are slightly different from
% the "nominal" ones.

% X0mat(:,3) = []; % for example, quickly remove the third scenario (conc 30)

Initial values for the model parameters

Model parameters are part of a 'structure' for easy reference. This means that you can address parameters by their name, instead of their position in a parameter vector. This makes the model definition in derivatives.m a lot easier. For each parameter, provide the initial value, whether you want to fit it or fix it to the initial value, the minimum bound, the maximum bound, and whether to fit the parameter on log10-scale or on normal scale. Fitting on log-scale is advisable for parameters that can span a very wide range; the optimisation routine can search this range more effectively if it is on log-scale. If you do not include this last value, a one will be filled in by prelim_checks.m (which means fitting on normal scale).

% syntax: par.name = [startvalue fit(0/1) minval maxval optional:log/normal scale (0/1)];
par.kd    = [0.1  1 1e-3 10  1];  % degradation rate constant, d-1
par.ke    = [0.2  1 1e-3 10  1];  % elimination rate constant, d-1
par.Piw   = [200  1 100  1e3 1];  % bioconcentration factor, L/kg
par.Ct    = [5    0 0    10  1];  % threshold external concentration where degradation stops, mg/L

% % Optionally, include an initial state as parameter (requires changes in call_deri.m too)
% par.Ci0   = [100  0 0 1e6];  % initial internal concentration (mg/L)

More options for the analysis

Optionally, global parameters can be used in the structure glo (not fitted). However, see the text file reserved_globals.txt for names to avoid.

% % specify globals for parameters that sre never fitted
% glo.Ct = 5; % the threshold could be made into a fixed global (also modify derivatives.m)
% % Special globals to modify the optimisation (used in transfer.m)
% glo.wts = [1 10];  % different weights for different data sets (same size as DATA)
% glo.var = [10 20]; % supply residual variance for each data set, after transformation (same size as DATA)

Zero-variate data and priors for Bayesian analyses

Optionally, zero-variate data can be included. The corresponding model value needs to be calculated in call_deri.m, so modify that file too! You can make up your own parameter names in the structure zvd. Optionally, prior distributions can be specified for parameters, see the file calc_prior.m for the definition of the distributions. You must use the exact same names for the prior parameters as used in the par structure. If you do not specify zvd and/or pri in your scripts, these options are simply bot used in the analysis.

zvd.ku = [10 0.5]; % example zero-variate data point for uptake rate, with normal s.d.

% First element in pri is the choice of distribution.
pri.Piw   = [2 116 121 118]; % triangular, min, max and center
pri.ke    = [3 0.2 0.1];     % normal with mean 0.2 and sd 0.1
% If no prior is defined, the min-max bounds will define a uniform one.
% Note that the prior is always defined on normal scale, also when the
% parameter will be fitted on log scale.

Time vector and labels for plots

Specify what to plot. This also involves the time vector for the model lines (which may be taken differently than the time vector in the data set). You can specify the text that you would like to use for the axes, and the text used for the legend (which always includes the identifier for the scenario: the values in the first row of X0mat).

glo.t = linspace(0,50,100); % time vector for the model curves in days

% specify the y-axis labels for each state variable
glo.ylab{1} = 'external concentration (mg/L)';
glo.ylab{2} = 'internal concentration (mg/kg)';
% specify the x-axis label (same for all states)
glo.xlab    = 'time (days)';
glo.leglab1 = 'conc. '; % legend label before the 'scenario' number
glo.leglab2 = 'mg/L'; % 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 functions are called that will do the calculation and the plotting. Note that calc_plot can provide all of the plotting information as output, so you can also make your own customised plots. This section, by default, makes a multiplot with all state variables (each in its own panel of the multiplot). When the model is fitted to data, output is provided on the screen (and added to the log-file results.out).

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.

You can turn on the events function there too, to smoothly catch the discontinuity in the model. For the demo, the iterations are turned off (opt_optim.it = 0).

glo.eventson   = 1; % events function on (1) or off (0)
opt_optim.it   = 0; % show iterations of the optimisation (1, default) or not (0)
opt_plot.annot = 1; % extra subplot in multiplot for fits: 1) box with parameter estimates, 2) overall legend

% 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

% save_plot(gcf,'fit_example',[],3) % save active figure as PDF (see save_plot for more options)
Provisional goodness-of-fit measures for each data set
    0.9613    0.9864

=================================================================================
Results of the parameter estimation with BYOM version 4.0
=================================================================================
   Filename      : byom_bioconc_extra 
   Analysis date : 27-Apr-2017 (11:52) 
   Data entered  :
     data state 1: 5x6, continuous data, no transf.
     data state 2: 5x6, continuous data, power 0.50 transf.
   Search method: Nelder-Mead simplex direct search, 2 rounds. 
     The optimisation routine has converged to a solution
     Total 150 simplex iterations used to optimise. 
     Minus log-likelihood has reached the value 228.377 (AIC=462.755). 
=================================================================================
kd        0.04356 (fit: 1, initial: 0.1) 	 prior: unif (0.001-10)
ke        0.09164 (fit: 1, initial: 0.2) 	 prior: norm (mn 0.2,sd 0.1)
Piw           118 (fit: 1, initial: 200) 	 prior: tri  (116-118-121)
Ct              5 (fit: 0, initial:  5) 	 prior: unif (0-10)
=================================================================================
Final values for zero-variate data points included:
ku          10.81 (data point: 10, sd: 0.5) 
=================================================================================
Time required: 4.5 secs
Plots result from the optimised parameter values. 

Local sensitivity analysis

Local sensitivity analysis of the model. All model parameters are increased one-by-one by a small percentage. The sensitivity score is by default scaled (dX/X p/dp) or alternatively absolute (dX p/dp).

Options for the sensitivity can be set using opt_sense (see prelim_checks.m).

% UNCOMMENT LINE(S) TO CALCULATE
% calc_localsens(par_out,opt_sens)

Asymptotic standard errors

Here, asymptotic standard errors are calculated from the curvature of the likelihood function in the best estimate. This involves a numerical method that is not very robust, and therefore, I do not recommend using this function for the final word on the parameter estimates (the two procedure below are far superior). However, it is a quick and dirty way to learn something about parameter uncertainty and correlations. The 'asymptotic' implies that these standard errors become more valid when the number of data points increases.

Options for the standard error can be set using opt_ase (see prelim_checks.m).

% UNCOMMENT LINE(S) TO CALCULATE
% calc_ase(par_out,opt_ase) % uncomment this line to calculate ase's

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. This example produces a profile for parameter kd and provides the 95% confidence interval (indicated by the horizontal broken line in the plot).

Notes: profiling for complex models is a slow process, so grab a coffee! If the profile finds a better solution, it immediately saves it to the log file profiles_newopt.out. You can break of the analysis by pressing ctrl-c anytime, and use the values from the log file to restart (copy the better values into your par structure).

Options for the profiling can be set using opt_prof (see prelim_checks.m).

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

% % UNCOMMENT LINE(S) TO CALCULATE
% % run a profile for selected parameters ...
% calc_proflik(par_out,'kd',opt_prof);  % uncomment to calculate a profile
% calc_proflik(par_out,'Piw',opt_prof); % uncomment to calculate a profile

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

% 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 not needed. Note that the profile
% for Piw is cut off due to the triangular prior distribution.
par_out = auto_profiles(par_out,opt_prof,opt_optim); % Experimental!
 
Starting automatic calculations 
 
Starting automatic calculations without sub-optimisations.
No results will be printed on screen or plotted until the analysis is finished.
  Starting a profile for parameter kd (1 of 3)
  Starting a profile for parameter ke (2 of 3)
  Starting a profile for parameter Piw (3 of 3)
Starting second run with sub-optimisations.
  Starting a profile for parameter kd (1 of 3)
  Starting a profile for parameter ke (2 of 3)
  Starting a profile for parameter Piw (3 of 3)
 
Parameters and bounds are for the last runs with 10 sub-optimisations.
Parameter, best fit value, interval
==========================================================
kd        0.04356 interval:    0.04093 - 0.04638 
ke        0.09164 interval:    0.08458 - 0.09883 
Piw           118 interval:      116.8 - 120.4 
Ct              5 parameter not fitted 
==========================================================
Time required: 11 mins, 44.6 secs
 
very similar profile for parameter kd   (max. diff. 0.0066) 
very similar profile for parameter ke   (max. diff. 0.065) 
very similar profile for parameter Piw  (max. diff. 0.009) 
 

Slice sampler

The slice samples can be used for a Bayesian analysis as it provides a sample from the posterior distribution. A .mat file is saved which contains the sample, to use later for e.g., intervals on model predictions. The output includes the Markov chain and marginal distributions for each fitted parameter. The function calc_conf can be used to put confidence intervals on the model lines. The function calc_ellipse makes an error ellipse for any two parameters that you specify.

Notes: MCMC is also slow process for complex models, so consider another coffee! If a prior is specified, it is plotted in the marginal posterior plot as a red distribution. Two types of credible intervals are generated:

  1. Highest probability regions. This captures the 95% of the parameter values with the highest posterior probability. it is calculated by lowering the horizotal dotted line in the figures for each parameters, until the area between the cut-off points captures 95% of the density. This can lead to unequal probabilities in the tails.
  2. Quantiles. The region between the 2.5 and 97.5% quantiles is taken. This leads to 2.5% probability in each tail. This is shown in the plots by the broken vertical lines.

Use opt_conf.type to tell calc_conf which sample to use, and how to use it: 1) Bayesian MCMC sample (default); CI on predictions are 95% ranges on the model curves from the sample 2) parameter sets from the 95% joint likelihood region; CIs on predictions are min-max of the curves from sample of the conf. region 3) as 1 but only 95% of the sample with the highest posterior density used; CIs on predictions are min-max of the curves from sample selection

Options for the slice sampling can be set using opt_slice (see prelim_checks.m). Options for confidence bounds on model curves and calculation of error ellipse can be set in opt_conf (see prelim_checks).

Taking on parameters on log scale helps, as the slice sampler seems quite bad at obtaining a representative sample when the parameters have very different ranges. Thinning is needed to reduce autocorrelation. Note the plot that is produced with details of the sample!

% % UNCOMMENT LINE(S) TO CALCULATE
% opt_slice.thin     = 5; % thinning of the sample (keep one in every 'thin' samples)
% opt_slice.burn     = 100; % number of burn-in samples (0 is no burn in)
% opt_slice.alllog   = 1; % set to 1 to put all parameters on log-scale before taking the sample
% calc_slice(par_out,200,opt_slice); % second argument number of samples (-1 to re-use saved sample from previous runs)

Likelihood region

Another way to make intervals on model predictions is to use a sample of parameter sets taken from the joint likelihood-based conf. region. This is done by the function calc_likregion.m. It first does profiling of all fitted parameters to find the edges of the region. Then, Latine-Hypercube sampling, keeping only those parameter combinations that are not rejected at the 95% level in a lik.-rat. test.

Options for the likelihood region can be set using opt_likreg (see prelim_checks.m).

opt_likreg.detail   = 2; % detailed (1) or a course (2) calculation
opt_likreg.subopt   = 0; % number of sub-optimisations to perform to increase robustness
calc_likregion(par_out,200,opt_likreg); % second argument number of samples (-1 to re-use saved sample from previous runs)
 
Calculating profiles and sample from confidence region ... please be patient.
  
Confidence interval from the profile (single parameter, 95% confidence)
=================================================================================
kd         interval:    0.04093 - 0.04638 
ke         interval:    0.08458 - 0.0988 
Piw        interval:      116.8 - 120.4 
  
Edges of the joint 95% confidence region (using df=p)
=================================================================================
kd         interval:     0.0398 - 0.04781 
ke         interval:    0.08132 - 0.1018 
Piw        interval:      116.8 - 120.4 
 
Starting with obtaining a sample from the joint confidence region.
Bursts of 100 samples, until at least 200 samples are accepted.
Currently accepted: 0, 52, 106, 158, 
 
Number of latin-hypercube samples taken: 400, of which accepted in the conf. region: 208
 
Time required: 12 mins, 55.2 secs

Plot results with confidence intervals

The following code can be used to make a standard plot (the same as for the fits), but with confidence intervals (and sampling error for survival data, if needed). Options for confidence bounds on model curves and calculation of error ellipse can be set using opt_conf (see prelim_checks).

Use opt_conf.type to tell calc_conf which sample to use, and how to use it: 1) Bayesian MCMC sample (default); CI on predictions are 95% ranges on the model curves from the sample 2) parameter sets from the 95% joint likelihood region; CIs on predictions are min-max of the curves from sample of the conf. region 3) as 1 but only 95% of the sample with the highest posterior density used; CIs on predictions are min-max of the curves from sample selection

opt_conf.type = 2; % use the values from the slice sampler (1) or likelihood region (2) to make intervals
out_conf = calc_conf(par_out,opt_conf); % calculate confidence intervals on model curves
calc_and_plot(par_out,opt_plot,out_conf); % call the plotting routine again to plot fits with CIs
 
Calculating CIs on model curves, using sample from joint likelihood region ... please be patient.
Percentage of samples finished:    0  10  20  30  40  50  60  70  80  90
Time required: 12 mins, 59.2 secs
Plots result from the optimised parameter values. 

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_bioconc_extra.html byom_bioconc_extra.m>
% and <byom_bioconc_start.html byom_bioconc_start.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_byom.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.

Cw = X(1); % state 1 is the external concentration at previous time point
Ci = X(2); % state 2 is the internal concentration 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.

kd   = par.kd(1);   % degradation rate constant, d-1
ke   = par.ke(1);   % elimination rate constant, d-1
Piw  = par.Piw(1);  % bioconcentration factor, L/kg
Ct   = par.Ct(1);   % threshold external concentration that stops degradation, mg/L

% % Optionally, include the threshold concentration as a global (see script)
% Ct   = glo.Ct;      % threshold external concentration that stops degradation, mg/L

%% Calculate the derivatives
% This is the actual model, specified as a system of two ODEs:
%
% $$ \frac{d}{dt}C_w=-k_d C_w \quad \textrm{as long as } C_w>C_t $$
%
% $$ \frac{d}{dt}C_i=k_e(P_{iw}C_w-C_i) $$

if Cw > Ct % if we are above the critical internal concentration ...
    dCw = -kd * Cw; % let the external concentration degrade
else           % otherwise ...
    dCw = 0; % make the change in external concentration zero
end

dCi = ke * (Piw*Cw-Ci); % first order bioconcentration

dX = [dCw;dCi]; % collect both 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] = 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_byom.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
% Ci0   = par.Ci0(1); % example: parameter for the internal concentration
% X0(2) = Ci0; % 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 (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. 

% 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); 
%     % 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, it looks for the
% external concentration where degradation stops. This function should be
% adapted to the problem you are modelling (this one matches the 
% byom_bioconc_... files). You can catch more events by making a vector out
% of _values_.
%
% 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)

Ct      = par.Ct(1); % threshold external concentration where degradation stops
nevents = 1;         % number of events that we try to catch

value       = zeros(nevents,1); % initialise with zeros
value(1)    = X(1) - Ct;        % thing to follow is external concentration (state 1) minus threshold
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