# Getting started with BYOM

Back to the overview of the manual

## Contents

## Installation

To install BYOM, unpack the ZIP file to a location of your choice. Do *not* rename any of the directories. The first files to study are in the `examples` directory. It is best not to modify these files, so you can always return to a working model. To create your own model, copy the files in `examples` to a new directory, somewhere as sub- or sub-sub-directory of the `BYOM` directory (do *not* rename the `BYOM` directory, put all new directories *below* this one, and don't start any of the directory names with `BYOM`.

Make sure that the new directory includes a copy of the functions `derivatives.m` (which holds your model equations), `call_deri.m` (calls the derivatives functions with initial state values and includes the events function), and `pathdefine.m` (makes sure the directory `../BYOM/engine` is added to the Matlab path). Modify `derivatives.m` to represent your own models. In several cases, (slight) modifications of `call_deri.m` will be needed. There is no need to change the functions in the engine folder, and in general I would advise to leave them alone unless you know what you're doing.

I suggest naming your script files starting with `byom_` to clearly distinguish them from the functions. Do *not* include spaces in script or directory names (Matlab does not like that). The engine directory contains the files for fitting the model to data; in general, there should be no need to modify these files.

## The example in this manual

To explain how BYOM works, I am first running through a very simple example. This example concerns biocontration with the standard one-compartment model. The internal concentration (*Ci*) is given by this ODE:

Here, *Cw* is the external concentration in water, *Kiw* is the partition coefficient between water and the organism (the bioconcentration factor), and *ke* is the elimination rate.

We are going to fit this model to a data set from Gomez et al. (2012), http://dx.doi.org/10.1007/s11356-012-0964-3 on accumulation of two exposure levels of tetrazepam in the mussel *Mytilus galloprovincialis*. Seven days of exposure exposure is followed by seven days of depuration in clean water. The same data set was used for the modelling chapter in the Marine Ecotoxicology textbook, but here we use the means rather than the individual replicates (see notes at the end of this page for comments on the statistics).

First, I will go through the the important elements of BYOM in isolation, and then show what it looks like in the BYOM script.

lets_get_started

## What you need to put in

**Entering data.** The data are entered in the main script in a cell array, like so:

% Internal concentration in mussels in ug/kg dwt, time in days, exposure % concentrations in ug/L. DATA{1} = [1 2.3 14.5 0 0 0 1 52.542 634.1 3 92.128 917.68 7 192.17 1438.8 8 71.975 783.81 10 33.108 356.28 14 6.4778 128.48];

The structure `DATA` holds all of your data in a specific format. In BYOM data sets are always linked to state variables. Our model has only one state (the internal concentration *Ci*), so we place our data set in location one in the cell array (hence `DATA{1}`). If the model has more state variables, you can add additional data sets as `DATA{2}` etc. As `DATA` is a cell array, the different data sets do not need to have the same size. If you have no data for a state, leave it out or just set it to zero.

The first element of the data matrix (top left) is for a flag how to treat the data. The 1 implies that it is continuous data, and that we don't want to use data transformation. The optimisation criterion is derived from the normal likelihood function, and thus assumes independent normal distributions, with constant variance, for the residuals. In this simple case, this approach is fully equivalent to least squares.

The other elements of the top row are the identifiers for the two exposure scenarios; in this case, the external concentrations in water during the exposure phase. These numbers are handed over to `call_deri` and `derivatives` as the variable `c`, and they are used for the legends in the plots. Below the top row are the time vector (in the first column) and the observations.

When we ask Matlab what the data is, we see that it is a cell array, containing a single cell, which in turn contains an 8x3 matrix with numbers:

DATA

DATA = cell [8×3 double]

Note that a data set may contain duplicate time points, and duplicate exposure treatments; no problem.

**Initial values.** Every state values needs to have initial values (the value at *t*=0) for each of the scenarios we like to run. These are the same identifiers `c` as we used in the data set. In BYOM, this is taken care of by the matrix `X0mat`:

X0mat = [2.3 14.5 % the scenarios (here exposure concentrations) 0 0 ]; % initial values state 1 (internal concentrations)

This tells BYOM that we want to run two scenarios (the two we also identified in `DATA`), and that the initial internal concentration is zero in both treatments. Note that `X0mat` may contain less scenarios than the data, or more. Only the scenarios in `X0mat` will be run, and the corresponding model curves plotted. If there is any data for these scenarios, they will be used for fitting and shown in the plots. If not, only a simulation is provided.

**Specifying parameters.** In BYOM, model parameters are part of a structure. This is very useful, as we can address them in `derivatives` and elsewhere by a name, which reduces the chance of errors. In your main script, you need to specify the parameters, provide them with a starting value, tell BYOM whether it needs to fit the parameter or fix it, and provide min-max bounds:

% syntax: par.name = [startvalue fit(0/1) minval maxval]; par.Kiw = [50 1 0 1e6]; % bioconcentration factor, L/kg par.ke = [0.2 1 1e-3 100]; % elimination rate constant, d-1

We use wide min-max bounds in this case. For rate constants, like *ke*, it is generally a good idea to make sure they do not run off to zero or infinity (that can lead to numerical issues such as extreme slowness of the calculations).

The nice thing about using a structure is that we can refer to parameters by their name, instead of its position in a vector. When we ask what `par` looks like after this:

par

par = struct with fields: Kiw: [50 1 0 1000000] ke: [0.2 1 0.001 100]

And we can also ask for a specific parameter by using its name:

par.ke

ans = 0.2 1 0.001 100

**Settings for modelling and plotting.** BYOM will choose its own time vector, based on the time points in your data sets. What I do need to specify is what text I would like to see with the axes of my plots:

% specify the y-axis labels for each state variable glo.ylab{1} = ['internal concentration (',char(181),'g/kg dwt)']; % specify the x-axis label (same for all states) glo.xlab = 'time (days)';

The specification of the y-label is a cell (hence the `{1}`). This is handy as each state can have its own y-axis label. There is only one x-label, which is usually time. The `char(181)` is a way to get the micro symbol for the units in there.

Next, specify what you would like to do with the legends. The scenarios (the first row in `X0mat`) will always be part of the legend, but there can be some text in front and after it:

glo.leglab1 = 'conc. '; % legend label before the 'scenario' number glo.leglab2 = [char(181),'g/L']; % legend label after the 'scenario' number

This will give us a nice legend when we're plotting. Now, let's look what the global variable `glo` now contains:

glo

glo = struct with fields: ylab: {'internal concentration (µg/kg dwt)'} xlab: 'time (days)' leglab1: 'conc. ' leglab2: 'µg/L'

**The model.** All of the information above had to be entered in the main script. The model, however, goes in the file derivatives.m. There, you first need to unpack the parameters from their structure (that leads to a more readable code):

Kiw = par.Kiw(1); % bioconcentration factor, L/kg ke = par.ke(1); % elimination rate constant, d-1

Note that I only want to know the first element of the parameter vector; its value. The other information is used elsewhere in the model, and we don't have to care about that.

The scenario (the first rows `X0mat`) enters `derivatives` through the variable `c`. In this case, our scenario identifiers are the concentration in the exposure phase, so that makes our life easy. However, the exposure is not constant: there is only 7 days exposure, after which the exposure is zero. To accommodate this, we use an `IF THEN` construction (note that time enters as `t`):

if t>7 % for the depuration phase ... c = 0; % make the external concentration zero end

The model is then entered in the form of the ODE (as specified above), in Matlab style:

dCi = ke * (Kiw * c - Ci); % first-order bioconcentration dX = [dCi]; % collect derivatives in vector dX to return to call_deri

It should be noted that in `derivatives`, the concentration `c` is always a single number, and time `t` as well (they enter one by one). That's all! Now we are ready to fit the model and make nice plots.

## What you will get out

So, what will you get as output from BYOM with this model and this data set? To show that, I first need to run the first part of a BYOM script (don't worry, this will be explained later in detail):

```
byom_my_first_pre % run the initial parts of the BYOM analysis
```

**Simulate first.** Now, all of the input steps above are done, and we are ready to fit. However, before fitting, it is *always* a good idea to start with a simulation. This is a good way to see if you entered the model correctly, that your data are plotted as you wanted, and that the starting values are not too far off. We can tell BYOM to simulate by the calling the function `calc_and_plot`, and feed it with the parameter structure `par` and optional options in `opt_plot`:

```
calc_and_plot(par,opt_plot); % calculate model lines and plot them
```

Plots are simulations using the user-provided parameter set. Parameter values used ===================================== Kiw 50 ke 0.2 ===================================== Minus log-likelihood 106.698. ===================================== Filename: byom_my_first_pre

The model runs, and the plot makes sense. The model lines are clearly not optimal fits, but they are close enough for the optimisation routine to do its job. Note the output on screen: BYOM tells you what it did, with which parameters, and what the minus log-likelhood (the optimisation criterion) is at this point.

**Then, optimise.** Now, let's run the optimisation routine by calling `calc_optim` (the `opt_optim` is used for setting optional options, don't worry about them now):

opt_optim.it = 0; % show iterations of the optimisation (1, default) or not (0) % Don't show iterations in this demo, as it will fill the page with % optimisation steps. In general, it is good to leave them on (which is the % defaul setting), especially when the model runs are slow. % optimise and plot (fitted parameters in par_out) par_out = calc_optim(par,opt_optim); % start the optimisation

Provisional goodness-of-fit measures for each data set 0.97601 ================================================================================= Results of the parameter estimation with BYOM version 4.0 ================================================================================= Filename : byom_my_first_pre Analysis date : 02-May-2017 (17:37) Data entered : data state 1: 7x2, continuous data, no transf. Search method: Nelder-Mead simplex direct search, 2 rounds. The optimisation routine has converged to a solution Total 69 simplex iterations used to optimise. Minus log-likelihood has reached the value 84.1322 (AIC=172.264). ================================================================================= Kiw 95.99 (fit: 1, initial: 50) ke 0.4625 (fit: 1, initial: 0.2) ================================================================================= Time required: 1.3 secs

This is the output report of the optimisation routine; it gives you everything you wanted to know (and probably more). The same output is also saved to a log file (`results.out`), so you can always look back at older runs (good idea to throw this file away once in a while, as it continues to grow).

The goodness-of-fit measure shown is the r-square, after any transformation have been applied (in this case, we did not transform). The output provides a number of identifiers for your reference, such as BYOM version, file name and date/time. The output shows what type of data you entered, which optimisation was used, and the final minus-log-likelihood achieved (the lower, the better the fit).

Now, the variable `par_out` contains our optimised parameter set. We can check to make sure:

par_out

par_out = struct with fields: Kiw: [95.988 1 0 1000000 1] ke: [0.4625 1 0.001 100 1] tag_fitted: 'y'

Note the extra tag that is added to let subsequent routines know that these parameters resulted from an optimisation. You could use these values as the new starting values for another optimisation (in general, that is a good idea, especially for parameter-rich models). You could manually copy the output from screen to your script (in the `par` structure with initial values), but there is a smarter way, by invoking a little function:

```
print_par(par_out) % display the contents of par_out in a formatted way
```

par.Kiw = [ 95.988 1 0 1e+06 1]; par.ke = [ 0.4625 1 0.001 100 1];

Now you can directly use copy-paste to make a new `par` matrix in your script.

**Plot the best fit.** Now that we have the optimised parameter set in `par_out`, we can plot it with the data, by calling the plotting routine `calc_and_plot` again, but now with the optimised parameter structure `par_out`:

```
calc_and_plot(par_out,opt_plot); % calculate model lines and plot them
```

Plots result from the optimised parameter values.

The plot looks nice. By default, BYOM plots in color, using different colors for different scenarios (first row in `X0mat`). However, there is a plotting option to plot in black and white, using different markers for different scenarios.

**Confidence intervals.** Every parameter estimate requires a confidence interval. I advise you to use the profile likelihood to construct these intervals. This is a robust method that will not only give you intervals but can also often tell if there is a better optimum to be had elsewhere in parameter space.

To make a profile of parameter *Kiw*:

calc_proflik(par_out,'Kiw',opt_prof); % calculate and plot the profile likelihood

95% confidence interval from the profile ================================================================================= Kiw interval: 88.86 - 103.8 Time required: 7.0 secs

The routine returns the 95% confidence interval for *Kiw*, which it has calculated from this graph. An explanation will follow in a later chapter, but note that the value of *Kiw* is plotted on the x-axis. The open circle denotes the best value as identified from the optimisation. The black points are new optimisations, fixing *Kiw* each time to values away from the best fit value. The horizontal line is a cut-off criterion: all values that are below this line are part of the confidence interval.

## Statistical and other notes

The analysis above makes sense, but there are a few notes and statistical warnings that I would like to add for completeness:

- It is a bit silly to use an ODE solver to tackle the equation for the one-compartment model given above, as an analytical solution can easily be derived. However, as soon as models become more meaningful, they are rapidly beyon analytical solutions. And, this simple model is a good example to demonstrate how BYOM operates.
- The data points were means of 1-3 replicate measurements. Therefore, I should have added weights in the analysis (the number of individual measurements underlying each data point). However, it would be even better to use the individual measurements themselves. Both are possible in BYOM, but will be explained later.
- The assumption of indepence of the residuals is allright: each measurement results from a new individual as the analysis of body residues is destructive. Normality may also be allright, but constant variance (homoscedasticity) is doubtful. It is likely that small values have less variance than large ones. Therefore, I would suggest a form of transformation on this data set (e.g., log-transformation). This is also possible with BYOM, and explained later.
- The profiling is way more robust than asymptotic standard errors, and relies less stringently on 'asymptotic' theory (i.e., it works well with small data sets). However, the cut-off criterion
*is*based on asymptotic theory, so the intervals will be less reliable with smaller data sets. This is nit-picking, as people generally rely without second thought on the much less reliable 'standard error' and the intervals they yield. - The likelihood function for the normally-distributed residuals requires an estimate of the residual variance (which we don't know). This variance is treated as a nuisance parameter, and profiled out. The different treatments within a data set are assumed to share the same residual variance, but different data sets have their own variance. This generally works well, but given sufficient interest I will be tempted to include more refined options. Note that BYOM also allows you to supply a residual variance per data set.
- The optimisation routine is a Nelder-Mead simplex search (
`fminsearch`in Matlab). Two runs are performed in sequence, the first with a low number of maximum iterations, which passes its optimised parameters to a second (default setting) one. I found this to be a useful enough strategy to implement as default. - BYOM offers a lot of options to modify its functions. You saw them in the Matlab calls above as
`opt_optim`,`opt_plot`and`opt_prof`. These options are structures that are filled automatically by the BYOM script. So, you only overwrite options if you like to depart from default behaviour. All options are described in a text file`all_options.txt`in the BYOM folder. - BYOM offers several other approaches for making confidence intervals and plotting results. E.g., there is a possibility to obtain a Bayesian sample from the posterior distribution, and create intervals on model curves. These will be treated later.

## Continu to next step

Next: byom_my_first.m to walk through this analysis in a formal way (Matlab-published version of the full script in a BYOM setting).