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:

$$ \frac{dC_i}{dt} = k_e (C_w K_{iw} - C_i) \quad \textrm{with } C_i(0)=0 $$

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:

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