BYOM, byom_gutsred_timevar_diazinon.m, reduced GUTS model

BYOM is a General framework for simulating model systems.

The model: fitting survival data with the GUTS special cases based on the reduced model (TK and damage dynamics lumped): SD, IT and mixed (or GUTS proper).

This script: dataset from Ashauer et al. (2010), http://dx.doi.org/10.1021/es903478b. Gammarus pulex exposed to pulses of diazinon. One scenario from this data set (two pulses) was used in Jager (2014), http://dx.doi.org/10.1007/s10646-013-1149-7. This script demonstrates two ways to define time-varying exposure scenarios: as a continuous spline or as a series of episodes.

Contents

Initial things

Make sure that this script is in a directory somewhere below the BYOM folder.

clear, clear global % clear the workspace and globals
global DATA W X0mat % make the data set and initial states global variables
global glo          % allow for global parameters in structure glo
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:

% observed number of survivors, time in days, conc in nM
DATA{1} = [-1     0     1     2     3
     0    60    70    70    70
     1    55    66    65    65
     2    53    61    59    59
     3    51    55    56    55
     4    51    31    54    53
     5    51    31    50    51
     6    51    29    47    48
     7    48    26    46    46
     8    48    24    46    46
     9    46    22    40    46
    10    45    21    23    44
    11    44    19    22    41
    12    42    17    22    40
    13    41    14    21    40
    14    41    14    18    40
    15    40    13    17    39
    16    38    11    17    38
    17    38    11    13    36
    18    37    10    13    33
    19    37     9    13    28
    20    37     8    11    24
    21    37     8    11    23
    22    36     8    11    19];

% In this data set, exposure was time-varying and reported as a series of
% concentrations over time. BYOM offers several ways to deal with it, as
% show below.
%
% Define the exposure scenarios as a series of pulses with a new constant
% exposure in each time period. This allow use of an analytical solution
% for the TK/damage dynamics in simplefun. This is a simplification of the
% reported exposure scenario, which has almost, but not quite, constant
% pulses in each episode.

% Cw1 = [0	1 % scenario 1
%     0	100.12
%     1	0
%     3	101.04
%     4	0];
% Cw2 = [0	2 % scenario 2
%     0	103.55
%     1	0
%     8	99.69
%     9	0];
% Cw3 = [0	3 % scenario 3
%     0	97.61
%     1	0
%     16	98.55
%     17	0];
%
% make_scen(2,Cw1,Cw2,Cw3); % create the globals to define the forcing function

% We can also use the exposure pattern exactly as reported. There are two
% ways to use this in BYOM: create a smooth spline (and use the ODE solver)
% or as a linear forcing series (which has an anlytical solution, and is
% thus much faster than the previous one).
Cw1 = [0    1
    0     102.65
    1.02	97.59
    1.03	0
    2.99	0
    3.01	103.88
    4.01	98.19
    4.02	0
    22.01	0];
Cw2 = [0    2
    0     100.78
    1.02	106.32
    1.03	0
    8     0
    8.01	103.56
    9     95.82
    9.01	0
    22.01	0];
Cw3 = [0    3
    0     100.6
    1.02	94.61
    1.03	0
    16    0
    16.01	100.58
    17    96.51
    17.01	0
    22.01	9.85];

make_scen(4,Cw1,Cw2,Cw3); % prepare as linear-forcing function interpolation (set glo.use_ode = 0)
% make_scen(1,Cw1,Cw2,Cw3); % prepare as cubic spline for interpolation (set glo.use_ode = 1)

% internal concentrations
DATA{2} = 0; % no data for scaled body residues

Initial values for the state variables

Initial states, scenarios in columns, states in rows. First row are the 'names' of all scenarios.

X0mat(1,:) = [0 1 2 3]; % scenarios (pulses)
X0mat(2,:) = 1;         % initial survival probability
X0mat(3,:) = 0;         % initial internal concentration

Initial values for the model parameters

Model parameters are part of a 'structure' for easy reference.

% global parameters for GUTS purposes
glo.sel  = 1; % select death mechanism: 1) SD 2) IT 3) mixed
glo.locD = 2; % location of scaled damage in the state variable list
glo.locS = 1; % location of survival probability in the state variable list

% Best values SD
% syntax: par.name = [startvalue fit(0/1) minval maxval optional:log/normal scale (0/1)];
par.kd    = [0.083 1 1e-3  10 0]; % dominant rate constant (d-1)
par.mw    = [4.8   1    0 1e6 1]; % median threshold for survival (nM)
par.hb    = [0.027 1    0 1e6 1]; % background hazard rate (1/d)
par.bw    = [0.023 1    0 1e3 0]; % killing rate (nM-1 d-1) (SD and mixed)
par.Fs    = [2     1    1 100 1]; % fraction spread of NEC distribution (-) (IT and mixed)

switch glo.sel % make sure that right parameters are fitted
    case 1 % SD
        par.Fs(2) = 0; % never fit the NEC spread!
    case 2 % IT
        par.bw(2) = 0; % never fit the killing rate!
    case 3 % mixed
        % do nothing: fit all parameters
end

Time vector and labels for plots

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

% specify the y-axis labels for each state variable
glo.ylab{1} = 'survival probability';
glo.ylab{2} = 'scaled damage (nM)';
% specify the x-axis label (same for all states)
glo.xlab    = 'time (days)';
glo.leglab1 = 'scenario '; % legend label before the 'scenario' number
glo.leglab2 = ''; % legend label after the 'scenario' number

prelim_checks % script to perform some preliminary checks and set things up
% Note: prelim_checks also fills all the options (opt_...) with defauls, so
% modify options after this call, if needed.

Calculations and plotting

Here, the function is called that will do the calculation and the plotting.

glo.useode = 0; % use analytical solution in simplefun (0) or ODE solver (1)
glo.stiff  = 1; % use ode113: works best with small step sizes needed for splined scenarios
if glo.sel == 1 % for the pure SD model, turn on the events function to catch threshold exceedance
    glo.eventson = 1; % events function on (1) or off (0)
end

opt_plot.annot = 2; % annotations in multiplot for fits: 1) box with parameter estimates 2) single legend
opt_optim.fit  = 1; % fit the parameters (1), or don't (0)
opt_optim.it   = 0; % show iterations of the optimisation (1, default) or not (0)

% 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

% opt_guts.doseresp  = 0; % set to zero as a dose-response is meaningless for time-var. exposure
% plot_guts(par_out,[],opt_guts); % make GUTS multiplots, but without CIs
Goodness-of-fit measures for each data set (R-square)
    0.9763

Warning: R-square is not appropriate for survival data and needs to be
interpreted more qualitatively! 
 
=================================================================================
Results of the parameter estimation with BYOM version 4.2b
=================================================================================
   Filename      : byom_gutsred_timevar_diazinon 
   Analysis date : 06-Sep-2018 (13:47) 
   Data entered  :
     data state 1: 23x4, survival data.
     data state 2: 0x0, no data.
   Search method: Nelder-Mead simplex direct search, 2 rounds. 
     The optimisation routine has converged to a solution
     Total 114 simplex iterations used to optimise. 
     Minus log-likelihood has reached the value 692.624 (AIC=1393.25). 
=================================================================================
kd        0.08359 (fit: 1, initial: 0.083) 
mw          4.675 (fit: 1, initial: 4.8) 
hb        0.02601 (fit: 1, initial: 0.027) 
bw        0.02281 (fit: 1, initial: 0.023) 
Fs              2 (fit: 0, initial:  2) 
=================================================================================
Parameters kd and bw are fitted on log-scale.
=================================================================================
Time required: 0.6 secs
Plots result from the optimised parameter values. 

Profiling the likelihood

By profiling you make robust confidence intervals for one or more of your parameters.

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

% % UNCOMMENT FOLLOWING 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,'mw',opt_prof);  % uncomment to calculate a profile

Likelihood region

Make intervals on model predictions by using a sample of parameter sets taken from the joint likelihood-based conf. region.

% UNCOMMENT FOLLOWING LINE(S) TO CALCULATE
opt_likreg.subopt = 10; % number of sub-optimisations to perform to increase robustness
opt_likreg.detail = 2; % detailed (1) or a coarse (2) calculation
opt_likreg.burst  = 2000; % number of random samples from parameter space taken every iteration
calc_likregion(par_out,500,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.05594 - 0.1411 
mw         interval:       3.01 - 6.203 
hb         interval:    0.01845 - 0.03496 
bw         interval:    0.01232 - 0.03564 
  
Edges of the joint 95% confidence region (using df=p)
=================================================================================
kd         interval:    0.04409 - 0.207 
mw         interval:     0.9108 - 8.043 
hb         interval:    0.01361 - 0.04091 
bw         interval:   0.006432 - 0.04467 
 
Starting with obtaining a sample from the joint confidence region.
Bursts of 2000 samples, until at least 500 samples are accepted in inner rim.
Currently accepted (inner rim): 0, 11, 22, 30, 37, 47, 51, 57, 66, 78, 
91, 97, 108, 118, 127, 136, 146, 157, 166, 177, 
187, 194, 210, 217, 229, 238, 246, 255, 264, 275, 
284, 293, 302, 311, 318, 331, 341, 346, 351, 362, 
370, 378, 387, 398, 406, 414, 423, 435, 445, 453, 
462, 473, 483, 494, 
 
Size of total sample (profile points have been added): 108085
of which accepted in the conf. region: 4458, and in inner rim: 551
 
Time required: 9 mins, 6.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; Bayes only).

The plot_guts function makes multiplots for the survival data, which are more readable when plotting with various intervals.

opt_conf.type    = 2; % use the values from the slice sampler (1) or likelihood region (2) to make intervals
opt_conf.lim_set = 2; % for lik-region sample: use limited set of n_lim points (1) or outer hull (2) to create CIs

% UNCOMMENT FOLLOWING LINE(S) TO CALCULATE AND PLOT
% Calculate and plot extra plots, specifically for GUTS.
opt_guts.doseresp = 0; % set to zero as a dose-response is meaningless for time-var. exposure
plot_guts(par_out,opt_conf,opt_guts);
 
Amount of 0.4 added to the chi-square criterion for inner rim
  Outer hull of 240 sets used.
Calculating CIs on model predictions, 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: 9 mins, 7.6 secs
 
Amount of 0.4 added to the chi-square criterion for inner rim
  Outer hull of 240 sets used.
Calculating CIs on model predictions, 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: 9 mins, 8.6 secs
============================================================================
Test for goodness-of-fit (Pearson chi-square with MC simulated distribution).
  Observed fit is better than this fraction of samples (p-value): 0.001
  Fit rejected at a=0.05; the data are not consistent with the calibrated model.
============================================================================

Calculate LPx

We can calculate an exposure multiplication factor: with which factor do we need to multiply an exposure scenario to obtain x% effect at the end of the scenario? Here, done for 10% effect after 2-22 days for scenario number 2 (as defined in the data set).

opt_conf.type    = 2; % use the values from the slice sampler (1) or likelihood region (2) to make intervals
opt_conf.lim_set = 2; % for lik-region sample: use limited set of n_lim points (1) or outer hull (2) to create CIs

% UNCOMMENT FOLLOWING LINE(S) TO CALCULATE
% this is the slow and general method that can work with a time vector
opt_lc50.Feff = 0.10; % effect level (>0 en <1), x/100 in LCx
opt_lc50.Fexp = [1 2]; % set first value to 1 to calculate exposure factor instead of LCx, second value is scenario identifier in X0mat
Tend          = [2 5 10 15 22]; % times at which to calculate exposure factor, relative to control
calc_conf_lc50(par_out,Tend,opt_conf,opt_lc50); % calculates LPx values, CI requires that there is a mat file with sample

% this is the new fast method that only does on t at a time
opt_lpx.Feff      = 0.10; % effect level (>0 en <1), x/100 in LCx
opt_lpx.scen_type = 4; % type of definition for the exposure scenario (as used in make_scen)
Tend              = 22; % time at which to calculate exposure factor, relative to control (this function cannot work with a T vector!)
calc_conf_lpx_lim(par_out,Tend,2,opt_conf,opt_lpx);
 
Calculate LPx for scenario 2
  Note that first column in X0mat is used for starting values
=================================================================================
LP10 (   2 days): 1.029 
LP10 (   5 days): 0.7948 
LP10 (  10 days): 0.6357 
LP10 (  15 days): 0.5421 
LP10 (  22 days): 0.5421 
=================================================================================
 
Amount of 0.4 added to the chi-square criterion for inner rim
  Outer hull of 240 sets used.
Calculating CIs on model predictions, using sample from joint likelihood region ... please be patient.
 
Percentage of sets finished for CIs:    0  10  20  30  40  50  60  70  80  90
=================================================================================
LP10 (   2 days): 1.029 (0.8533-1.203) 
LP10 (   5 days): 0.7948 (0.5868-0.9166) 
LP10 (  10 days): 0.6357 (0.4649-0.7181) 
LP10 (  15 days): 0.5421 (0.3874-0.6141) 
LP10 (  22 days): 0.5421 (0.3874-0.6141) 
=================================================================================
Intervals: 95% confidence (likelihood region)
 
Time required: 9 mins, 30.6 secs
 
Calculate exposure multiplication factors
==================================================================
LP10 (  22 days): 0.5421 
==================================================================
 
Amount of 0.4 added to the chi-square criterion for inner rim
  Outer hull of 240 sets used.
Calculating CIs on model predictions, using sample from joint likelihood region ... please be patient.
 
Percentage of sets finished for CIs:   0  10  20  30  40  50  60  70  80  90
 
Results table for LPx (-) with 95% CI, scenario 2
==================================================================
LP10 (  22 days): 0.5421 (0.3873-0.6141) 
==================================================================
Intervals: 95% confidence (likelihood region)
 
 
Amount of 0.4 added to the chi-square criterion for inner rim
  Outer hull of 240 sets used.
Calculating CIs on model predictions, 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: 9 mins, 31.5 secs
Time required: 9 mins, 31.6 secs

Calculate LCx versus time

Here, the LCx (by default the LC50) is calculated at several time points.

opt_conf.type    = 2; % use the values from the slice sampler (1) or likelihood region (2) to make intervals
opt_conf.lim_set = 2; % for lik-region sample: use limited set of n_lim points (1) or outer hull (2) to create CIs

% UNCOMMENT FOLLOWING LINE(S) TO CALCULATE
% We can directly use the fast method specific for this reduced GUTS folder
opt_lcx.Feff     = 0.50; % effect level (>0 en <1), x/100 in LCx
Tend = [2:0.2:3 3.5 4:8 10 15 22]; % times at which to calculate LCx, relative to control
calc_conf_lcx_lim(par_out,Tend,opt_conf,opt_lcx); % calculates LCx values, CI requires that there is a mat file with sample

% We can also use the older slower version, but we need to remove all
% scenario information first. We used pulsed exposure scenarios to fit the
% parameters, but the LC50 calculations uses X0mat to transfer real
% concentrations to the model, not scenario numbers (this might interfere
% with previously defined scenarios).

% % UNCOMMENT FOLLOWING LINE(S) TO CALCULATE
% make_scen(-5,-1);  % remove all spline info from glo
% opt_lc50.Feff = 0.50; % effect level (>0 en <1), x/100 in LCx
% opt_lc50.Fexp = [0 2]; % turn off the calculation of the exposure factor, and move back to LCx
% Tend = [2:0.2:3 3.5 4:8 10 15 22]; % times at which to calculate LCx, relative to control
% calc_conf_lc50(par_out,Tend,opt_conf,opt_lc50); % calculates LCx values, CI requires that there is a mat file with sample
 
Calculate LCx values
 
Amount of 0.4 added to the chi-square criterion for inner rim
  Outer hull of 240 sets used.
Calculating CIs on model predictions, using sample from joint likelihood region ... please be patient.
 
Percentage of sets finished for CIs:   0  10  20  30  40  50  60  70  80  90
 
Results table for LCx,t (), with 95% CI
==================================================================
LC50 (   2 days): 247.8 (200.2-333.3) 
LC50 ( 2.2 days): 210.3 (170.9-280.3) 
LC50 ( 2.4 days): 181.4 (148.2-239.6) 
LC50 ( 2.6 days): 158.5 (130.2-207.7) 
LC50 ( 2.8 days): 140.1 (115.7-182.1) 
LC50 (   3 days): 125.1 (103.8-161.3) 
LC50 ( 3.5 days): 97.51 (81.81-123.4) 
LC50 (   4 days): 79.01 (66.98-98.30) 
LC50 (   5 days): 56.25 (48.59-67.99) 
LC50 (   6 days): 43.14 (37.84-51.49) 
LC50 (   7 days): 34.80 (30.77-41.06) 
LC50 (   8 days): 29.10 (25.94-33.98) 
LC50 (  10 days): 21.94 (19.81-25.16) 
LC50 (  15 days): 13.93 (12.73-15.52) 
LC50 (  22 days): 9.806 (8.452-11.14) 
==================================================================
Intervals: 95% confidence (likelihood region)
 
Time required: 9 mins, 33.3 secs

Calculate LPx for a monitoring profile

We can use the function calc_conf_lpx_lim to calculate an exposure multiplication factor: with which factor do we need to multiply an exposure scenario to obtain x% effect at the end of the scenario? Here, done for 10% effect at the end of the exposure scenario.

This is the fast method, which is closely linked to this GUTS model (reduced model, standard directory). Don't use if you modify the GUTS model in simplefun.m or derivatives.m, unless you know what you are doing. The function calc_conf_lc50 can also calculate LPx values. This works more generally, but much slower, and it requires some more tricked-out coding. Contact me if you feel the need to use it.

% UNCOMMENT FOLLOWING LINE(S) TO CALCULATE
opt_conf.type     = 2; % use values from slice sampler (1), likelihood region(2) to make intervals (-1 to skip intervals)
opt_conf.lim_set  = 1; % for lik-region sample: use limited set of n_lim points (1) or outer hull (2) to create CIs
opt_lpx.Feff      = 0.10; % effect level (>0 en <1), x/100 in LCx
opt_lpx.scen_type = 4; % type of definition for the exposure scenario (as used in make_scen)
glo.useode        = 0; % no need for ODE solver, as long as scen_type = 4 (linear forcing)

Tend  = []; % time at which to calculate LPx (empty is end of profile)
LPx_mon = calc_conf_lpx_lim(par_out,Tend,'profile_monit.txt',opt_conf,opt_lpx);
 
Calculate exposure multiplication factors
==================================================================
LP10 ( 107 days): 9.690e+04 
==================================================================
 
Amount of 0.4 added to the chi-square criterion for inner rim
  Limited set of 200 sets used.
Calculating CIs on model predictions, using sample from joint likelihood region ... please be patient.
 
Percentage of sets finished for CIs:   0  10  20  30  40  50  60  70  80  90
 
Results table for LPx (-) with 95% CI, scenario 1
==================================================================
LP10 ( 107 days): 9.690e+04 (6.525e+04-1.188e+05) 
==================================================================
Intervals: 95% confidence (likelihood region)
 
 
Amount of 0.4 added to the chi-square criterion for inner rim
  Limited set of 200 sets used.
Calculating CIs on model predictions, 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: 9 mins, 55.5 secs
Time required: 9 mins, 55.5 secs