# BYOM, byom_bioconc_extra.m, a detailed example

- Author: Tjalling Jager
- Date: April 2017
- Web support: http://www.debtox.info/byom.html
- Back to index walkthrough_byom.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. 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:

**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
- The data set
- Initial values for the state variables
- Initial values for the model parameters
- More options for the analysis
- Zero-variate data and priors for Bayesian analyses
- Time vector and labels for plots
- Calculations and plotting
- Local sensitivity analysis
- Asymptotic standard errors
- Profiling the likelihood
- Slice sampler
- Likelihood region
- Plot results with confidence intervals
- Other files: derivatives
- Other files: call_deri

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

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

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:

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