The openGUTS prototype

Table of contents

Contents

About

Main script for the openGUTS prototype. This prototype forms the basis of the standalone openGUTS software that is currently under development. This programme follows the workflow as laid down in the 2018 EFSA Scientific Opinion on TKTD modelling:

Input files for survival data and exposure profiles are read from the respective folders input_data and input_profile. Information from the calibration is saved into the folder output_sample. Figures and log files for the output can be found in the folder [output_report].

% =========================================================================
% Copyright (c) 2018-2019, Tjalling Jager (tjalling@debtox.nl).
% This file is part of the Matlab prototype of openGUTS.
%
% This program is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation, either version 3 of the License, or (at your
% option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.
%
% You should have received a copy of the GNU General Public License along
% with this program. If not, see <https://www.gnu.org/licenses/>.
% =========================================================================

BLOCK 1. Initial things

Make sure the [engine] folder (which holds the actual functions that perform the calculations) is in the Matlab search path. The script [initial_setup] takes care of some general setting up. Also, in this block, you can set a name for the analysis (used for identifying the output files), and set a switch to use a saved sample from a previous calibration rather than loading a file and calibrating it.

pathdefine    % places the [engine] folder in Matlab search path
initial_setup % calls a script that sets things up

rng(2.2) % seed the random number generator with a specific value (this overrides the shuffle in the [initial_setup])
% Fixing the seed is helpful for testing: it ensures that every run of this file will be exactly the same.

% Provide a name for the study to be used for the plots and saved files.
GLO.basenm       = 'propiconazole'; % name used for identifying the analysis
use_saved_sample = 0; % when set to 1, use saved calibration file rather than make a new calibration
GLO.fix_hb_cal   = 1; % when set to 1, fit [hb] to the control data for all calibration data sets, and keep fixed
% Note: for the validation data set, there is a separate switch for fixing [hb] in BLOCK 4.

BLOCK 2. Calibration

Based on the setting of [use_saved_sample], a saved calibration is shown (regenerated from the project MAT file) or a new calibarion is performed. The calibration function [calc_optim] can be used for both purposes. The Matlab GUI [uigetfile] is used to select files.

Note: the second and third input arguments for [calc_optim] are a switch for SD/IT (1/2) and a flag for optimisation (1) or show saved optimisation (2). The [plot_main] also has flags for SD/IT (1/2), and a second argument, which is a vector of two flags: for the stage of the analysis we are in, and for whether to show model curves or not.

if use_saved_sample == 1
    % BLOCK 2.1. Use saved calibration results.
    % Note: Optimised parameter values and sample for confidence region
    % will be read from a saved MAT file (selected by the Matlab GUI).
    % Therefore, there is no need for calibration (as long as the saved
    % file exists).
    %
    % Note: Calling [calc_optim] to show saved optimisation makes [DATAcal]
    % from the saved set global, as well as [GLO.basenm], [GLO.concunit]
    % and [GLO.fix_hb_cal]. These are all associated with the previous
    % calibration, so they should not be changed after calibration has been
    % performed.

    filename = uigetfile('output_sample/*.mat','Select a saved project file (selecting either SD or IT will use both)','MultiSelect','off'); % use Matlab GUI to select MAT file
    if ~iscell(filename) && numel(filename) == 1 && filename(1) == 0 % if cancel is pressed ...
        return % simply stop
    end
    GLO.basenm = filename(1:end-16); % remove the last part (_parspace_SD.mat) of the filename so that the basename remains
    % First, see if there is an output folder for this specific studyname, and otherwise create it.
    savestr = ['output_report',filesep,GLO.basenm];
    if ~exist(savestr,'dir')
        mkdir(savestr); % if not, create it
    end

    % Use [calc_optim] to recreate the plot of parameters space and [plot_main] to plot the fit (standard plot).
    calc_optim([],1,2); % SD just use saved set, without new optimisation
    calc_optim([],2,2); % IT just use saved set, without new optimisation
    plot_main(DATAcal,2,[1 1]) % IT plot of calibration data with the fit (stage 1, model on)
    plot_main(DATAcal,1,[1 1]) % SD plot of calibration data with the fit (stage 1, model on)

else
    % BLOCK 2.2. Load and prepare calibration data.
    % The data set will be loaded and translated to the internal
    % representation, which is returned in [DATAcal]. [DATAcal] is made
    % global (in [initial_setup]) as it is called many times during
    % optimisation.

    [filename,filepath] = uigetfile('input_data/*.txt','Select file(s) for calibration','MultiSelect','on'); % use Matlab GUI to select file(s)
    if ~iscell(filename) && numel(filename) == 1 && filename(1) == 0 % if cancel is pressed ...
        return % simply stop
    end

    [DATAcal,error_flag] = load_data(filename,filepath,1); % load and prepare the data set for stage 1 (returned as structured cell array, so [data_all] can include multiple data sets)
    if error_flag == 1 % see if there are errors spotted in the input file
        return % simply return (stop and do nothing)
    end
    plot_main(DATAcal,[],[1 0]) % inspection plot of data (stage 1, model off)
    % Note: the 1 in the calls above signal that we are loading, preparing
    % and plotting data for stage 1: calibrations.

    % BLOCK 2.3. Ranges for the model parameters. Note: the parameter
    % structure [pmat] will be automatically generated by [startgrid]
    % below. Experts can make changes to the [pmat] after that.
    %
    % WARNING: some settings of [pmat] will lead to poor performance or
    % even complete failure of the calibration algorithm. Therefore, when
    % messing with this meatrix, expertise and care are essential to get
    % reliable results.

    % Make separate parameter matrices for SD and IT.
    pmat_SD = startgrid(1); % create reasonable SD parameter ranges from the data set
    pmat_IT = startgrid(2); % create reasonable IT parameter ranges from the data set
    % At this place, modifications to the parameter matrices can be
    % inserted.

    % BLOCK 2.4. Perform calibration.
    % Note that the starting values in [pmat] are only used for fixed
    % parameters in the analysis; for fitted parameter, only the ranges and
    % other settings are used. [calc_optim] finds the global optimum,
    % confidence intervals on parameters, and the joint confidence region
    % (to be used for predictions). The sample is saved in a MAT file with
    % an extension indicating the analysis (SD or IT), which is
    % automatically used by the post-analyses.

    % Use [calc_optim] to perform the calibration and [plot_main] to plot the fit (standard plot).
    calc_optim(pmat_SD,1,1); % SD, new optimisation
    calc_optim(pmat_IT,2,1); % IT, new optimisation
    plot_main(DATAcal,1,[1 1]) % SD plot of calibration data with the fit (stage 1, model on)
    plot_main(DATAcal,2,[1 1]) % IT plot of calibration data with the fit (stage 1, model on)

end

% return % the return stops the analysis, so we can run this script in parts
Following data sets are loaded for calibration:
    set: 1, file: propiconazole_constant.txt
 
Settings for parameter search ranges:
=====================================================================
kd    bounds:   0.001641 -      143.8 1/d      fit: 1 (log)
mw    bounds:   0.002202 -      35.56 uM       fit: 1 (norm)
hb    bounds:    0.01307 -    0.01307 1/d      fit: 0 (norm)
bw    bounds:  0.0007332 -       9310 1/(uM d) fit: 1 (log)
Fs    bounds:          1 -          1 [-]      fit: 0 (norm)
=====================================================================
Special case: stochastic death (SD) 
Background hazard rate fixed to a value fitted on controls
 
Starting round 1 with initial grid of 2016 parameter sets
  Status: best fit so far is (minloglik) 131.8625
Starting round 2, refining a selection of 200 parameter sets, with 60 tries each
  Status: 16 sets within total CI and 8 within inner. Best fit: 125.8154
Starting round 3, refining a selection of 200 parameter sets, with 40 tries each
  Status: 268 sets within total CI and 89 within inner. Best fit: 125.8154
Starting round 4, refining a selection of 771 parameter sets, with 60 tries each
  Status: 7868 sets within total CI and 2672 within inner. Best fit: 125.8154
  Finished sampling, running a simplex optimisation ...
  Status: 7869 sets within total CI and 2673 within inner. Best fit: 125.8153
Starting round 5, creating the profile likelihoods for each parameter
  Finished profiling, running a simplex optimisation on the best fit set found ...
  Status: 7870 sets within total CI and 2674 within inner. Best fit: 125.8153
 
=================================================================================
Results of the parameter estimation
=================================================================================
   openGUTS prototype: v0.7 of 12 July 2019
   Base name    : propiconazole 
   Analysis date: 12-Jul-2019 (12:13) 
   Following data sets are loaded for calibration:
     set: 1, file: propiconazole_constant.txt
   Sample: 7870 sets in joint CI and 2674 in inner CI.
   Propagation set: 1225 sets will be used for error propagation.
   Minus log-likelihood has reached 125.8153 (AIC=257.63). 
Best estimates and 95% CIs on single parameters
==========================================================================
kd    best:      2.191 (     1.630  -     3.353 ) 1/d       fit: 1 (log)
mw    best:      16.93 (     15.74  -     17.57 ) uM        fit: 1 (norm)
hb    best:    0.01307 (       NaN  -       NaN ) 1/d       fit: 0 (norm)
bw    best:     0.1306 (   0.08627  -    0.1912 ) 1/(uM d)  fit: 1 (log)
Fs    best:      1.000 (       NaN  -       NaN ) [-]       fit: 0 (norm)
==========================================================================
Special case: stochastic death (SD) 
Background hazard rate fixed to a value fitted on controls
 
Extra information from the fit:
==========================================================================
Depuration/repair time (DRT95)    : 1.37 (0.893 - 1.84) days 
Parameter ranges used for the optimisation:
kd    range:   0.001641 -      143.8 
mw    range:   0.002202 -      35.56 
hb    range:    0.01307 -    0.01307 
bw    range:  0.0007332 -       9310 
Fs    range:          1 -          1 
==========================================================================
Time required for optimisation: 16.9 secs
 
 
Settings for parameter search ranges:
=====================================================================
kd    bounds:   0.001641 -      143.8 1/d      fit: 1 (log)
mw    bounds:   0.002202 -      71.85 uM       fit: 1 (norm)
hb    bounds:    0.01307 -    0.01307 1/d      fit: 0 (norm)
bw    bounds:        Inf -        Inf 1/(uM d) fit: 0 (norm)
Fs    bounds:          1 -         20 [-]      fit: 1 (log)
=====================================================================
Special case: individual tolerance (IT) 
Background hazard rate fixed to a value fitted on controls
 
Starting round 1 with initial grid of 2016 parameter sets
  Status: best fit so far is (minloglik) 131.5096
Starting round 2, refining a selection of 200 parameter sets, with 60 tries each
  Status: 17 sets within total CI and 7 within inner. Best fit: 127.9594
Starting round 3, refining a selection of 200 parameter sets, with 40 tries each
  Status: 298 sets within total CI and 106 within inner. Best fit: 127.9594
Starting round 4, refining a selection of 767 parameter sets, with 60 tries each
  Status: 8094 sets within total CI and 2833 within inner. Best fit: 127.9594
  Finished sampling, running a simplex optimisation ...
  Status: 8095 sets within total CI and 2834 within inner. Best fit: 127.9593
Starting round 5, creating the profile likelihoods for each parameter
  Finished profiling, running a simplex optimisation on the best fit set found ...
  Status: 8096 sets within total CI and 2835 within inner. Best fit: 127.9593
 
=================================================================================
Results of the parameter estimation
=================================================================================
   openGUTS prototype: v0.7 of 12 July 2019
   Base name    : propiconazole 
   Analysis date: 12-Jul-2019 (12:14) 
   Following data sets are loaded for calibration:
     set: 1, file: propiconazole_constant.txt
   Sample: 8096 sets in joint CI and 2835 in inner CI.
   Propagation set: 1294 sets will be used for error propagation.
   Minus log-likelihood has reached 127.9593 (AIC=261.92). 
Best estimates and 95% CIs on single parameters
==========================================================================
kd    best:     0.7292 (    0.5458  -    0.9359 ) 1/d       fit: 1 (log)
mw    best:      17.73 (     15.30  -     20.03 ) uM        fit: 1 (norm)
hb    best:    0.01307 (       NaN  -       NaN ) 1/d       fit: 0 (norm)
bw    best:        Inf (       NaN  -       NaN ) 1/(uM d)  fit: 0 (norm)
Fs    best:      1.714 (     1.511  -     2.047 ) [-]       fit: 1 (log)
==========================================================================
Special case: individual tolerance (IT) 
Background hazard rate fixed to a value fitted on controls
 
Extra information from the fit:
==========================================================================
Beta of the threshold distribution: 6.80 (5.11 - 8.88) (-) 
Depuration/repair time (DRT95)    : 4.11 (3.20 - 5.49) days 
Parameter ranges used for the optimisation:
kd    range:   0.001641 -      143.8 
mw    range:   0.002202 -      71.85 
hb    range:    0.01307 -    0.01307 
bw    range:        Inf -        Inf 
Fs    range:          1 -         20 
==========================================================================
Time required for optimisation: 21.2 secs
 
Goodness of fit measures for calibration
  Special case: stochastic death (SD)
=====================================================================
Model efficiency (NSE, r-square)          : 0.9811 
Normalised root-means-square error (NRMSE): 7.949 % 
Survival probability prediction error (SPPE) for each treatment
     Data set     treatment    value
        1         Control     +0.0937 % 
        1              T1     +0.0937 % 
        1              T2     -9.91 % 
        1              T3     -14.9 % 
        1              T4     +3.79 % 
        1              T5     -0.558 % 
        1              T6     -0.690 % 
        1              T7     -0.0297 % 
=====================================================================
Warning: these measures need to be interpreted more qualitatively
  as they, strictly speaking, do not apply to quantal data
 
Goodness of fit measures for calibration
  Special case: individual tolerance (IT)
=====================================================================
Model efficiency (NSE, r-square)          : 0.9605 
Normalised root-means-square error (NRMSE): 11.665 % 
Survival probability prediction error (SPPE) for each treatment
     Data set     treatment    value
        1         Control     +0.0937 % 
        1              T1     +0.395 % 
        1              T2     -5.77 % 
        1              T3     -4.40 % 
        1              T4     +21.1 % 
        1              T5     -9.25 % 
        1              T6     -4.71 % 
        1              T7     -1.12 % 
=====================================================================
Warning: these measures need to be interpreted more qualitatively
  as they, strictly speaking, do not apply to quantal data
 

BLOCK 3. Predictions: LCx versus time

LCx,t values are calculated for several effect percentages and several preset time points (defined in [initial_setup]). Optimised parameter values and sample to create CIs will be read from the saved MAT file.

calc_lcx(16,1); % SD, calculate LCx values with CI (final timepoint t=16)
calc_lcx(16,2); % IT, calculate LCx values with CI (final timepoint t=16)

% return % the return stops the analysis, so we can run this script in parts
Results table for LCx,t (uM), with 95% CI
  Special case: stochastic death (SD)
======================================================================
time (d)       LC50               LC20               LC10               
======================================================================
   1   33.3 ( 29.8 -  37.7)  25.5 ( 23.1 -  28.1)  23.1 ( 20.8 -  25.3) 
   2   22.6 ( 21.5 -  24.0)  19.5 ( 18.6 -  20.3)  18.6 ( 17.5 -  19.3) 
   3   20.1 ( 19.2 -  21.0)  18.2 ( 17.2 -  18.8)  17.7 ( 16.6 -  18.2) 
   4   19.0 ( 18.1 -  19.8)  17.7 ( 16.6 -  18.3)  17.4 ( 16.2 -  18.0) 
   7   17.9 ( 16.9 -  18.6)  17.3 ( 16.1 -  17.9)  17.1 ( 15.9 -  17.7) 
  14   17.4 ( 16.2 -  18.0)  17.1 ( 15.9 -  17.7)  17.0 ( 15.8 -  17.6) 
  21   17.2 ( 16.0 -  17.9)  17.0 ( 15.8 -  17.7)  17.0 ( 15.7 -  17.6) 
  28   17.1 ( 15.9 -  17.8)  17.0 ( 15.8 -  17.6)  17.0 ( 15.7 -  17.6) 
  42   17.1 ( 15.9 -  17.7)  17.0 ( 15.7 -  17.6)  17.0 ( 15.7 -  17.6) 
  50   17.0 ( 15.8 -  17.7)  17.0 ( 15.7 -  17.6)  16.9 ( 15.7 -  17.6) 
 100   17.0 ( 15.8 -  17.6)  16.9 ( 15.7 -  17.6)  16.9 ( 15.7 -  17.6) 
======================================================================
 
Results table for LCx,t (uM), with 95% CI
  Special case: individual tolerance (IT)
======================================================================
time (d)       LC50               LC20               LC10               
======================================================================
   1   34.3 ( 31.4 -  37.9)  27.9 ( 25.4 -  30.5)  24.8 ( 22.1 -  27.3) 
   2   23.1 ( 21.7 -  24.7)  18.8 ( 17.2 -  20.4)  16.7 ( 14.8 -  18.4) 
   3   20.0 ( 18.5 -  21.6)  16.3 ( 14.6 -  18.0)  14.5 ( 12.5 -  16.3) 
   4   18.7 ( 17.0 -  20.7)  15.3 ( 13.4 -  17.2)  13.6 ( 11.5 -  15.5) 
   7   17.8 ( 15.6 -  20.1)  14.6 ( 12.2 -  16.8)  12.9 ( 10.5 -  15.2) 
  14   17.7 ( 15.2 -  20.1)  14.5 ( 12.0 -  16.8)  12.8 ( 10.3 -  15.1) 
  21   17.7 ( 15.2 -  20.1)  14.5 ( 12.0 -  16.8)  12.8 ( 10.3 -  15.1) 
  28   17.7 ( 15.2 -  20.1)  14.5 ( 12.0 -  16.8)  12.8 ( 10.3 -  15.1) 
  42   17.7 ( 15.2 -  20.1)  14.5 ( 12.0 -  16.8)  12.8 ( 10.3 -  15.1) 
  50   17.7 ( 15.2 -  20.1)  14.5 ( 12.0 -  16.8)  12.8 ( 10.3 -  15.1) 
 100   17.7 ( 15.2 -  20.1)  14.5 ( 12.0 -  16.8)  12.8 ( 10.3 -  15.1) 
======================================================================
 

BLOCK 4. Validation: comparison to new data

The function [plot_main] can directly be used to make predictions (with CI) for new scenarios (not used for calibration) and compare them with data. As input, it needs a new data set, in the same way as for the calibration data set. This also implies that it is possible to use multiple data sets (although we allow only 1 for the standalone software). Optimised parameter values and confidence region will be read from a saved MAT file.

The data set is treated in the same manner as the calibration data. It is translated to the internal representation, which is returned in [data_val]. This definition is NOT made global; it is simply passed on to the plotting routine.

GLO.fix_hb_val = 1; % when set to 1, fit [hb] to the control data for all validation data sets, and keep fixed

% BLOCK 4.1. Load and prepare validation data.
[filename,filepath] = uigetfile('input_data/*.txt','Select file(s) for validation','MultiSelect','on'); % use Matlab GUI to select file(s)
if ~iscell(filename) && numel(filename) == 1 && filename(1) == 0 % if cancel is pressed ...
    return % simply stop
end

[data_val,error_flag] = load_data(filename,filepath,2); % load and prepare the data set for stage 2 (returned as structured cell array, so [data_all] can include multiple data sets)
if error_flag == 1 % see if there are errors spotted in the input file
    return % simply return (stop and do nothing)
end
plot_main(data_val,1,[2 0]) % inspection plot of data (stage 2, model off)
% Note: the 2 in the calls above signal that we are loading, preparing and
% plotting data for stage 2: validations.

% BLOCK 4.2. Call the plotting routine for validation plots (with model output).
plot_main(data_val,1,[2 1]) % SD validation: compare to model predictions (stage 2, model on)
plot_main(data_val,2,[2 1]) % IT validation: compare to model predictions (stage 2, model on)

% return % the return stops the analysis, so we can run this script in parts
Following data sets are loaded for validation:
    set: 1, file: propiconazole_pulsed_linear.txt
Goodness of fit measures for validation
  Special case: stochastic death (SD)
=====================================================================
Model efficiency (NSE, r-square)          : 0.5213 
Normalised root-means-square error (NRMSE): 14.822 % 
Survival probability prediction error (SPPE) for each treatment
     Data set     treatment    value
        1         Control     +0.00623 % 
        1      close pulses     +14.8 % 
        1      wide pulses     +11.9 % 
        1        constant     -12.9 % 
=====================================================================
Warning: these measures need to be interpreted more qualitatively
  as they, strictly speaking, do not apply to quantal data
 
Goodness of fit measures for validation
  Special case: individual tolerance (IT)
=====================================================================
Model efficiency (NSE, r-square)          : 0.7504 
Normalised root-means-square error (NRMSE): 10.701 % 
Survival probability prediction error (SPPE) for each treatment
     Data set     treatment    value
        1         Control     +0.00623 % 
        1      close pulses     -11.8 % 
        1      wide pulses     -20.6 % 
        1        constant     -12.8 % 
=====================================================================
Warning: these measures need to be interpreted more qualitatively
  as they, strictly speaking, do not apply to quantal data
 

BLOCK 5. Predictions: LPx

LPx is calculated for various effect levels (defined in [initial_setup]). This is the factor by which a given exposure scenario needs to be multiplied to result in x% mortality at the end of the profile, relative to the control. Optimised parameter values and confidence region will be read from the saved MAT file.

Exposure scenarios need to be entered as text files. However, only one treatment is allowed per data set. Multiple data sets can be used (although we allow only 1 for the standalone software). The scenario is thus a matrix with two columns: time (in days) and exposure concentration (in the same units as the data used for calibration!). This will be linearly interpolated.

Note: the second entry in the call to [calc_lpx] is a flag that is used for faster calculations (no CI) and batch processing. Note that flag setting 2 will probably not be used in the standalone software, but it is handy for testing (as CIs on LPx can take a very long calculation time). flag 1) regular, 2) no CI (very fast), 3) batch processing (no output to screen, no diary, no CIs)

% BLOCK 5.1. Load and prepare prediction data.
[filename,filepath] = uigetfile('input_profile/*.txt','Select file(s) with exposure profiles for standard predictions','MultiSelect','on'); % use Matlab GUI to select file(s)
if ~iscell(filename) && numel(filename) == 1 && filename(1) == 0 % if cancel is pressed ...
    return % simply stop
end

[data_pred,error_flag] = load_data(filename,filepath,3); % load and prepare the data set for stage 3 (returned as structured cell array, so [data_all] can include multiple data sets)
if error_flag == 1 % see if there are errors spotted in the input file
    return % simply return (stop and do nothing)
end
plot_main(data_pred,[],[3 0],[]) % inspection plot of data (stage 3, plot model off)
% Note: the 3 in the calls above signal that we are loading, preparing and
% plotting data for stage 3: predictions.

% BLOCK 5.2. Calculate LPx and plot.
sdit = 1; % select death mechanism: 1) SD 2) IT
LPx  = calc_lpx(data_pred,1,sdit); % calculates LPx values with CI (and makes a plot)
plot_main(data_pred,sdit,[3 1],[ones(size(LPx,1),1) LPx]) % plot of data (stage 3, plot model, without CIs), pass on the calculated LPx values
% Note: setting the flags in [plot_main] from [3 1] to [3 2] creates model
% curves with CIs (terribly slow for FOCUS profiles). Note: the vector
% [ones(size(LPx,1),1) LPx] in [calc_lpx] specifies which multiplication
% factors to use for plotting. I add the [ones] here to also plot the
% survival due to the unmodified exposure profile (as is also done in the
% standalone software).

% return % the return stops the analysis, so we can run this script in parts
Following exposure profiles are loaded for predictions:
    set: 1, file: profile_monit.txt
Results table for LPx (-) with 95% CI: profile_monit
  Special case: stochastic death (SD)
==================================================================
LPx     best       CI
==================================================================
LP10: 1.105e+05 (1.046e+05 - 1.148e+05)   
LP50: 1.382e+05 (1.315e+05 - 1.457e+05)   
==================================================================
Time required for the LPx calculations: 1 min, 28.7 secs
 

BLOCK 6. Batch predictions for LPx

Note that [calc_batch_lpx] returns the filename for the scenario that leads to the lowest (=worst) LP50 value. Note that the call to [calc_batch_lpx] includes a correction factor [corr_factor] that can be used to correct for different units of the exposure profile, if needed. We decided not to use that for the software.

% BLOCK 6.1. Select files with exposure profiles for batch prediction.
[filename,filepath] = uigetfile('input_profile/*.txt','Select file(s) with exposure profiles for batch predictions','MultiSelect','on'); % use Matlab GUI to select file(s)
if ~iscell(filename) && numel(filename) == 1 && filename(1) == 0 % if cancel is pressed ...
    return % simply stop
end

% Note: the EFSA opinion also uses the propiconazole data set as a case
% study. However, there is some confusion regarding the concentration
% units. The calibration data are reported as nmol/ml (or uM) whereas the
% exposure profiles are in ug/L. I intend to distribute these profiles with
% the prototype unchanged, to allow the users to reproduce the same numbers
% as in the opinion. Below, we can set a correction factor. However, these
% files are examples only, and not meant to be used for research or
% regulatory application.

corr_factor = 1; % to reconstruct the results in the EFSA opinion (and for the software, we don't use this correction factor)
% corr_factor = 1 / 342.2; % correction factor for ug/L (profiles) -> uM (calibration data)

% BLOCK 6.2. Run batch-calculation of all profiles and return worst one.
sdit       = 1; % select death mechanism: 1) SD 2) IT
scen_worst = calc_batch_lpx(filename,filepath,sdit,corr_factor); % call [calc_batch_lpx] to do the calculations, and return name of worst profile

% BLOCK 6.3. Detailed analysis for the worst one from the batch calculations.
% Note: this is identical to the procedure in Block 5, but using the
% filename in [scen_worst] rather than selecting files with the Matlab GUI.
[data_pred,error_flag] = load_data(scen_worst,filepath,3,corr_factor); % load and prepare the data set for stage 3 (returned as structured cell array, so [data_all] can include multiple data sets)
if error_flag == 1 % see if there are errors spotted in the input file
    return % simply return (stop and do nothing)
end

LPx = calc_lpx(data_pred,2,sdit); % calculates LPx values (1 with CI, 2 without) (and makes a plot)
plot_main(data_pred,sdit,[3 1],LPx) % plot of data (stage 3, plot model without CIs), pass on the calculated LPx values

% return
 
Starting batch processing ...
 
LPx values, scenarios ordered by descending risk (LP50)
  Special case: stochastic death (SD)
file analysed             LP10       LP50 
==================================================================
     cereal_D1_ditch      2.293      2.648 
      cereal_D5_pond      11.04      11.65 
      cereal_D4_pond      11.07      11.70 
     cereal_D3_ditch      8.595      12.25 
       apple_R1_pond      15.90      16.71 
    cereal_D1_stream      13.62      20.23 
    cereal_R4_stream      30.57      38.41 
     apple_R2_stream      31.84      43.40 
    cereal_D5_stream      113.6      193.3 
    cereal_D4_stream      119.9      203.9 
==================================================================
Following exposure profiles are loaded for predictions:
    set: 1, file: cereal_D1_ditch.txt
 
Results table for LPx (-) with 95% CI: cereal_D1_ditch
  Special case: stochastic death (SD)
==================================================================
LP10 at end of profile: 2.293 
LP50 at end of profile: 2.648 
==================================================================

BLOCK 7. Demonstration of predictions for test design

The ability to make predictions can also be used for efficient test design. Here, an example is given. The user should prepare one or more scenarios of choice, e.g., with various pulses. With MF, the user can specify multiplication factors to plot, and [plot_main] will plot the associated damage and survival over time (here with CIs).

For testing, after calibration with propiconazole data for constant exposure, try the two test profiles in the directory [pulses_test_design]. These also work well with the dieldrin case study.

[filename,filepath] = uigetfile('input_profile/*.txt','Select file(s) with exposure profiles for test-design predictions','MultiSelect','on'); % use Matlab GUI to select file(s)
if ~iscell(filename) && numel(filename) == 1 && filename(1) == 0 % if cancel is pressed ...
    return % simply stop
end

[data_pred,error_flag] = load_data(filename,filepath,3); % load and prepare the data set for stage 3 (returned as structured cell array, so [data_all] can include multiple data sets)
if error_flag == 1 % see if there are errors spotted in the input file
    return % simply return (stop and do nothing)
end

sdit = 1; % select death mechanism: 1) SD 2) IT
MF   = [1 2 2.5 5]; % these are the multiplication factors to plot
MF   = ones(data_pred.nr_sets,1) * MF; % make sure there is the same MF for all profiles that are entered

plot_main(data_pred,sdit,[3 2],MF) % plot of data (stage 3, plot model, with CIs), with the different MF values
Following exposure profiles are loaded for predictions:
    set: 1, file: test_close_pulse.txt
    set: 2, file: test_wide_pulse.txt