# BYOM, byom_guts_start.m

- Author: Tjalling Jager
- Date: April 2017
- Web support: http://www.debtox.info/byom.html
- Back to index walkthrough_guts.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. Results are shown on screen but also saved to a log file (`results.out`).

**The model:** fitting survival data with the two simplest GUTS special cases: the reduced model with pure SD and IT (log-logistic). For faster calculations (when exposure is constant), the analytical solution in simplefun.m is used. To use the ODE version, set the option `glo.useode= 1` (which uses the model in derivatives.m). For more elaborate GUTS scripts, take a look at the examples in the "reduced" and "full" directories.

The equation for scaled damage (referenced to water) is used both by SD and IT:

The SD model uses the hazard rate, which is calulated as:

The hazard rate is integrated over time to yield the survival probability (using the analystical solution), or it is included in ODE form (in dervatives.m):

The IT model first finds the maximum of the scaled damage over time:

And then uses the cumulative distribution of the thresholds to find the survival probability:

**This script:** Long acute toxicity test for guppy (*Poecilia reticulata*) exposed to the insecticide dieldrin. Data from: Bedaux and Kooijman (1994), http://dx.doi.org/10.1007/BF00469427. This scrip deals with the basics only; many options for post-calculations are provided in the byom_guts_extra script.

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

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

% observed number of survivors, time in days, conc. in ug/L DATA{1} = [ -1 0 3.2 5.6 10 18 32 56 100 0 20 20 20 20 20 20 20 20 1 20 20 20 20 18 18 17 5 2 20 20 19 17 15 9 6 0 3 20 20 19 15 9 2 1 0 4 20 20 19 14 4 1 0 0 5 20 20 18 12 4 0 0 0 6 20 19 18 9 3 0 0 0 7 20 18 18 8 2 0 0 0 ]; % Note: For survival data, the weights matrix W has a different meaning % than for continuous data: it can be 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). By default, the weights matrix for survival data is filled with % zeros. % scaled damage, can have no observations DATA{2} = 0;

## Initial values for the state variables

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

X0mat(1,:) = DATA{1}(1,2:end); % scenarios (concentrations) X0mat(2,:) = 1; % initial survival probability X0mat(3,:) = 0; % initial scaled damage

## 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 glo.locD = 2; % location of scaled damage in the state variable list glo.locS = 1; % location of survival probability in the state variable list % syntax: par.name = [startvalue fit(0/1) minval maxval optional:log/normal scale (0/1)]; par.kd = [0.5 1 1e-3 10 1]; % dominant rate constant, d-1 par.mw = [5 1 0 1e6 1]; % median threshold for survival (ug/L) par.hb = [0.01 1 1e-4 1 1]; % background hazard rate (1/d) par.bw = [0.05 1 1e-4 1e6 1]; % killing rate (L/ug/d) (SD only) par.Fs = [3 1 1 100 1]; % fraction spread of threshold distribution (IT only) switch glo.sel % make sure that right parameters are fitted case 1 % for SD ... par.Fs(2) = 0; % never fit the threshold spread! case 2 % for IT ... par.bw(2) = 0; % never fit the killing rate! 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 (',char(181),'g/L)']; % specify the x-axis label (same for all states) glo.xlab = 'time (days)'; glo.leglab1 = 'conc. '; % legend label before the 'scenario' number glo.leglab2 = [char(181),'g/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 function is called that will do the calculation and the plotting. 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.

For the demo, the iterations are turned off (opt_optim.it = 0).

glo.useode = 0; % use the analytical solution in simplefun.m (constant exposure) % par = start_vals(par); % experimental start-value finder; use at your own risk % % Note: in this way, start_vals will overwrite the parameter structure par! 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

Provisional goodness-of-fit measures for each data set 0.9501 ================================================================================= Results of the parameter estimation with BYOM version 4.0 ================================================================================= Filename : byom_guts_start Analysis date : 27-Apr-2017 (12:33) Data entered : data state 1: 8x8, 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 103 simplex iterations used to optimise. Minus log-likelihood has reached the value 161.527 (AIC=331.053). ================================================================================= kd 0.7911 (fit: 1, initial: 0.5) mw 5.205 (fit: 1, initial: 5) hb 0.008352 (fit: 1, initial: 0.01) bw 0.03761 (fit: 1, initial: 0.05) Fs 3 (fit: 0, initial: 3) ================================================================================= Time required: 1.4 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. 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.

Options for the profiling can be set using opt_prof (see prelim_checks). For this demo, the level of detail of the profiling is set to 'course' and no sub-optimisations are used. However, consider these options when working on your own data (e.g., set opt_prof.subopt=10).

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 FOLLOWING LINE(S) TO CALCULATE % % run a profile for selected parameters ... % calc_proflik(par_out,'kd',opt_prof); % calculate a profile % calc_proflik(par_out,'mw',opt_prof); % calculate a profile % or run profiles for all fitted parameters. all_profiles(par_out,opt_prof);

95% confidence interval from the profile ================================================================================= kd interval: 0.493 - 1.317 Time required: 13.8 secs 95% confidence interval from the profile ================================================================================= mw interval: 4.102 - 6.897 Time required: 24.4 secs 95% confidence interval from the profile ================================================================================= hb interval: 0.002146 - 0.02159 Time required: 38.5 secs 95% confidence interval from the profile ================================================================================= bw interval: 0.02638 - 0.05423 Time required: 51.4 secs

## Calculate LCx versus time

Here, the LCx (by default the LC50) is calculated at several time points. When sufficient points are specified, a smooth line for LCx versus time will be produced. LCx values are also printed on screen. If a sample from parameter space is available (either from the slice sampler or the likelihood region), it can be used to calculate confidence bounds.

Options for LCx (with confidence bounds) can be set using opt_lc50 (see prelim_checks).

opt_lc50.type = 0; % use the values from the slice sampler (1) or likelihood region (2) to make intervals Tend = [1:0.2:3 3.5 4:8]; % times at which to calculate LCx, relative to control calc_conf_lc50(par_out,Tend,opt_lc50); % calculates LCx values, CI requires that there is a mat file with sample

Calculate LCx values ================================================================================= LC50 ( 1 days): 75.7432 LC50 ( 1.2 days): 57.3192 LC50 ( 1.4 days): 45.7563 LC50 ( 1.6 days): 37.9606 LC50 ( 1.8 days): 32.4188 LC50 ( 2 days): 28.3156 LC50 ( 2.2 days): 25.1778 LC50 ( 2.4 days): 22.7148 LC50 ( 2.6 days): 20.7392 LC50 ( 2.8 days): 19.1255 LC50 ( 3 days): 17.7868 LC50 ( 3.5 days): 15.2778 LC50 ( 4 days): 13.5454 LC50 ( 5 days): 11.3396 LC50 ( 6 days): 10.0146 LC50 ( 7 days): 9.14071 LC50 ( 8 days): 8.52548 ================================================================================= No confidence intervals calculated Time required: 52.4 secs

## 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_guts1.html byom_guts1.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_guts1.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. S = X(1); % state 1 is the survival probability at previous time point Dw = X(2); % state 2 is the scaled damage 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); % dominant rate constant, d-1 mw = par.mw(1); % median of threshold distribution, ug/L bw = par.bw(1); % killing rate, L/ug/d Fs = par.Fs(1); % fraction spread of threshold distribution (used in call_deri) hb = par.hb(1); % background hazard rate, d-1 %% Calculate the derivatives % This is the actual model, specified as a system of ODEs. % % *Note:* scaled damage and background hazard are calculated here for every % special case. However, survival is done here only for the pure SD model. % Pure IT is calculated in <call_deri.html call_deri.m>. dDw = kd * (c - Dw); % first order damage build-up from c (scaled) switch glo.sel case 1 % stochastic death h = bw * max(0,Dw-mw); % calculate the hazard rate dS = -(h + hb)* S; % change in survival probability (incl. background mort.) case 2 % individual tolerance dS = -hb* S; % only background hazard rate % mortality due to the chemical is included in call_deri! case 3 error('Use the GUTS models in the "full" or "reduced" folder for the IT+SD model!') end dX = [dS;dDw]; % collect both derivatives in one vector

## 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_guts1.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 %% 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 % smaller stepsize is a good idea for pulsed exposures; otherwise, stepsize % may become so large that a pulse is missed completely TE = 0; % dummy for time of events % For the GUTS packages, we need a long time vector for IT and for 'proper % GUTS' (SD and IT combined). For the IT this is needed when damage can % also decrease in time: we have to find the maximum of the damage over % time. For the proper model, we need to numerically integrate the hazard % rate over time here. If the time span of the data set or simulation is % very long, consider increasing the size of t below. t_rem = t; % remember the original time vector if glo.sel ~= 1 && length(t)<100 t = unique([t;(linspace(min(t),max(t),100))']); end 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. % % Stochastic death is calculated in simplefun.m or in derivatives.m % directly. For the individual tolerance calculation, we need a few extra % things here. % Note that for IT, a long time vector will be generated above, which is % needed to catch the maximum dose metric over time during pulsed exposure. if glo.sel == 2 % for individual tolerance ... mw = par.mw(1); % median of thhreshold distribution, ug/L Fs = par.Fs(1); % fraction spread of the threshold distribution Dw = Xout(:,glo.locD); % take the correct state variable as scaled damage % Dw(Dw<=0) = 1e-20; % take care scaled damage <= 0 (probably not needed) maxDw = zeros(size(Dw)); % intitialise a new vector for maximum Dw over time for j = 1:length(t), % this section makes sure that Dw does not decrease in time (dead animals dont become alive) % Note: for IT, the time vector will always be long (see above) maxDw(j,:) = max(Dw(1:j,:)); % so look for max of Dw over time end % Log-logistic distribution is used for GUTS-IT beta = log(39)/log(Fs); % translate Fs to beta S = (1 ./ (1+(maxDw/mw).^beta)) .* Xout(:,glo.locS); % survival probability % the survival due to the chemical is multiplied with the background % survival (calculated in derivatives), assuming logistic distrib. Xout(:,glo.locS) = S; % replace correct state by newly calculated survival prob. [~,loct] = ismember(t_rem,t); % find where the requested time points are in the long Xout Xout = Xout(loct,:); % only keep the ones we asked for end %% Events function % Modify this part of the code if _eventson_=1. % This subfunction catches the 'events': in this case, it looks for the % internal concentration where the threshold is exceeded (only when % useode=1). This only makes sense for SD. % % 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) global glo mw = par.mw(1); % median for thresholds nevents = 1; % number of events that we try to catch value = zeros(nevents,1); % initialise with zeros value(1) = X(glo.locD) - mw; % thing to follow is damage 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