BYOM, ERA-special-Daphnia with simple-compound model: byom_debtox_daphnia.m

Table of contents

Contents

About

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: simple DEBtox model for toxicants, based on DEBkiss and formulated in compound parameters. The model includes flexible modules for toxicokinetics/damage dynamics and toxic effects. The DEBkiss e-book (see http://www.debtox.info/book_debkiss.html) provides a partial description of the model; the publication of Jager in Ecological Modelling contains the full details: https://doi.org/10.1016/j.ecolmodel.2019.108904.

This script: The water flea Daphnia magna exposed to fluoranthene; data set from Jager et al (2010), http://dx.doi.org/10.1007/s10646-009-0417-z. Published as case study in Jager & Zimmer (2012). Simultaneous fit on growth, reproduction and survival. Parameter estimates differ slightly from the ones in Jager & Zimmer due to several model changes and choices in the analysis. Furthermore, here the raw data on individual females is used, such that the special functions for censoring the reproduction data can be used.

Copyright (c) 2012-2021, Tjalling Jager, all rights reserved.
This source code is licensed under the MIT-style license found in the
LICENSE.txt file in the root directory of BYOM.

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
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(1) % set path to the BYOM/engine directory (option 1 uses parallel toolbox)
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:

% NOTE: time MUST be entered in DAYS for the estimation of starting values
% to provide proper search ranges! Controls must use identifiers 0 for the
% true control and 0.1 for the solvent control. If you want to use another
% identifier for the solvent, change parameter id_solvent in
% automatic_runs.

TR   = 1; % transformations for continuous data
opt  = 2; % select an option (opt=1 is recommended for Daphnia)

% Options to deal with repro data set:
% 0) Check whether we can use a single intermoult period for the entire
%    data set. Screen output will show mean intermoult times and brood
%    sizes across the replicates, as function of brood number and treatment.
% 1) Cumulate reproduction, but remove the time points with zero
%    reproduction. Good for clutch-wise reproduction.
% 2) Cumulate reproduction, but don't remove zeros. Good for continuous
%    reproduction or when animals are not followed individually.
% 3) Shift neonate release back to the previous moult. When this option is
%    used, don't shift the model predictions with <glo.Tbp>: the data now
%    represent egg production rather than neonate release.

% Put the position of the various states in globals, to make sure that the
% correct one is selected for extra things (e.g., for plotting in
% plot_tktd, in call_deri for accommodating 'no shrinking', for population
% growth rate).
glo.locD = 1; % location of scaled damage in the state variable list
glo.locL = 2; % location of body size in the state variable list
glo.locR = 3; % location of cumulative reproduction in the state variable list
glo.locS = 4; % location of survival probability in the state variable list

fnames = {'data_FLU_Dmagna'};
read_datafiles(fnames,TR,opt); % helper function defines DATA, W and glo.LabelTable

% Note: optionally, add some info to the MAT filename. The MAT filename
% will already include MoA and feedbacks, but if you want to try other
% things as well (changing opt, calibrating on the validation data, etc)
% it can be helpful to change the name to use this script but not
% overwrite previous MAT files.

% glo.basenm  = [mfilename,'_CAL1']; % remember the filename for THIS file for the plots
save([glo.basenm,'_DATA'],'DATA','W') % save MAT file with data set
% Saving the data set is handy to allow for a simple reconstruction of the
% calibrations, without needing to define the data again, in the same way.
% Note that if you use DATAx and Wx, you need to save them as well!

switch opt
    case 0 % check intermoult duration (and mean brood size)
        return % we need to stop to check the results (no data are created)
    case 3 % then we have shifted the data set, so no shifting of model needed
        glo.Tbp = 0; % time that eggs spend in brood pouch (model output for repro will be shifted)
    otherwise % we need to shift the model output
        glo.Tbp = 3; % time that eggs spend in brood pouch (model output for repro will be shifted)
        % glo.Tbp = 0; % in some cases, effect may be on the eggs in the pouch, so no shift proposed
end

Initial values for the state variables

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

X0mat      = [glo.LabelTable.Scenario]'; % the scenarios (here identifiers)
X0mat(2,:) = 0; % initial values state 1 (scaled damage)
X0mat(3,:) = 0; % initial values state 2 (body length, initial value overwritten by L0)
X0mat(4,:) = 0; % initial values state 3 (cumulative reproduction)
X0mat(5,:) = 1; % initial values state 4 (survival probability)

Initial values for the model parameters

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

% global parameters as part of the structure glo
glo.FBV    = 0.02;    % dry weight egg as fraction of structural body weight (-) (for losses with repro; approx. for Daphnia magna)
glo.KRV    = 1;       % part. coeff. repro buffer and structure (kg/kg) (for losses with reproduction)
glo.kap    = 0.8;     % approximation for kappa (for starvation response)
glo.yP     = 0.8*0.8; % product of yVA and yAV (for starvation response)
glo.Lm_ref = 5;       % reference max length for scaling rate constants
glo.len    = 2;       % switch to fit length 1) with shrinking, 2) without shrinking (used in call_deri.m)
% NOTE: the settings above are species specific! Make sure to use the same
% settings for validation and prediction as for calibration!
%
% NOTE: For arthropods, one would generally want to fit the model without
% shrinking (since the animals won't shrink in length).

% syntax: par.name = [startvalue fit(0/1) minval maxval optional:log/normal scale (0/1)];
par.L0   = [1     1 0.5  1.5  1]; % body length at start experiment (mm)
par.Lp   = [1.8   1 0.5  4.0  1]; % body length at puberty (mm)
par.Lm   = [3     1 2    7    1]; % maximum body length (mm)
par.rB   = [0.15  1 0.02 0.3  1]; % von Bertalanffy growth rate constant (1/d)
par.Rm   = [10    1 5    30   1]; % maximum reproduction rate (#/d)
par.f    = [1     0 0    2    1]; % scaled functional response
par.hb   = [0.01  1 1e-3 0.07 0]; % background hazard rate (d-1)
% - Note 1: it does not matter whether length measures are entered as actual
% length or as volumetric length (as long as the same measure is used
% consistently).
% - Note 2: hb is fitted on log-scale. This is especially helpful for
% fit_tox(1)=-2 when hb is fitted along with the other parameters. The
% simplex fitting has trouble when one parameter is much smaller than
% others.
% - Note 3: the min-max ranges in par are appropriate for Daphnia magna. For
% other species, these ranges and starting values need to be modified.

glo.names_sep = {}; % no parameters can differ between data sets
% glo.names_sep = {'f';'L0'}; % names of parameters that can differ between data sets
% par.L01   = [0.9 1 0.5  1.5  1]; % body length at start experiment (mm)
% par.f1    = [0.9 1 0    2    1]; % scaled functional response

% Note: using separate parameters for separate data sets requires using
% specific identifiers, and using exposure scenarios with make_scen (also
% for constant exposure and for the controls)!

% extra parameters for special situations
par.Lf   = [0 0 0 1e6 1]; % body length at half-saturation feeding (mm)
par.Lj   = [0 0 0 1e6 1];  % body length at end acceleration (mm)
par.Tlag = [0 0 0 1e6 1];  % lag time for start development

ind_tox = length(fieldnames(par))+1; % index where tox parameters start
% the parameters below this line are all treated as toxicity parameters!

% When using the parameter-space explorer, start values and ranges are not
% needed (filled later by startgrid_debtox); only the fit/fix mark is
% relevant here. But make sure that the value in the first column is within
% the bounds (and not zero for log-scale parameters).
par.kd   = [0.08   1 0.01  10 0]; % dominant rate constant (d-1)
par.zb   = [0.1    1 0    1e6 1]; % effect threshold energy budget ([C])
par.bb   = [75     1 1e-6 1e6 0]; % effect strength energy-budget effect (1/[C])
par.zs   = [0.23   1 0    1e6 1]; % effect threshold survival ([C])
par.bs   = [0.7    1 1e-6 1e6 0]; % effect strength survival (1/([C] d))

% After optimisation, copy-paste relevant lines (fitted parameters) from screen below!

Time vector and labels for plots

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

% specify the y-axis labels for each state variable
glo.ylab{1} = ['scaled damage (',char(181),'g/L)'];
glo.ylab{2} = 'body length (mm)';
if isfield(glo,'Tbp') && glo.Tbp > 0
    glo.ylab{3} = ['cumul. repro. (shift ',num2str(glo.Tbp),'d)'];
else
    glo.ylab{3} = 'cumul. repro. (no shift)';
end
glo.ylab{4} = 'survival fraction (-)';

% 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),'M']; % legend label after the 'scenario' number
% Note: these legend labels will not be used when we make a glo.LabelTable

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.

NOTE: for this package, the options useode and eventson in glo will not be functional: the ODE solver is always used, and the events function as well.

% -------------------------------------------------------------------------
% Configurations for the ODE solver
glo.stiff = [0 3]; % ODE solver 0) ode45 (standard), 1) ode113 (moderately stiff), 2) ode15s (stiff)
% Second argument is for default sloppy (0), normally tight (1), tighter
% (2), or very tight (3) tolerances. Use 1 for quick analyses, but check
% with 3 to see if there is a difference! Especially for time-varying
% exposure, there can be large differences between the settings!
glo.break_time = 0; % break time vector up for ODE solver (1) or don't (0)
% Note: breaking the time vector is a good idea when the exposure scenario
% contains discontinuities. Don't use for continuous splines (type 1) as it
% will be much slower. For FOCUS scenarios (high time resolution), breaking
% up is not efficient and does not appear to be necessary.
% -------------------------------------------------------------------------

opt_optim.fit    = 1; % fit the parameters (1), or don't (0)
opt_plot.bw      = 0; % if set to 1, plots in black and white with different plot symbols
opt_plot.annot   = 2; % annotations in multiplot for fits: 1) box with parameter estimates 2) single legend
opt_plot.repls   = 0; % set to 1 to plot replicates, 0 to plot mean responses
basenm_rem       = glo.basenm; % remember basename as automatic_runs may modify it!

% Select what to fit with fit_tox (this is a 3-element vector). NOTE: use
% identifier c=0 for regular control, and c=0.1 for solvent control. When
% entering more than one data set, use 100 and 100.1 for the controls of
% the second data set (and 101, 102 ... for the treatments), 200 and 200.1
% for the controls of the third data set. etc.
%
% First element of fit_tox is which part of the data set to use:
%   fit_tox(1) = -2  comparison between control and solvent control (c=0 and c=0.1)
%   fit_tox(1) = -1  control survival (c=0) only
%   fit_tox(1) = 0   controls for growth/repro (c=0) only, but not for survival
%   fit_tox(1) = 1   all treatments, but, when fitting, keep all control parameters fixed;
%               run through all elements in MOA and FEEDB sequentially and
%               provide a table at the end (plots are made incl. control)
%
% Second element of fit_tox is whether to fit or only to plot:
%   fit_tox(2) = 0   don't fit; for standard optimisations, plot results for
%               parameter values in [par], for parspace optimisations, use saved mat file.
%   fit_tox(2) = 1   fit parameters
%
% Third element is what to use as control (fitted for fit_tox(1) = -1 or 0)
% (if this element is not present, only regular control will be used)
%   fit_tox(3) = 1   use regular control only (identifier 0)
%   fit_tox(3) = 2   use solvent control only (identifier 0.1)
%   fit_tox(3) = 3   use both regular and control (identifier 0 and 0.1)
%
% The strategy in this script is to fit hb to the control data first. Next,
% fit the basic parameters to the control data. Finally, fit the toxicity
% parameter to the complete data set (keeping basic parameters fixed. The
% code below automatically keeps the parameters fixed that need to be
% fixed.
%
% MOA: Mode of action of toxicant as set of switches
% [assimilation/feeding, maintenance costs (somatic and maturity), growth costs, repro costs]
% [1 0 0 0 0]   assimilation/feeding
% [0 1 0 0 0]   costs for maintenance
% [0 0 1 1 0]   costs for growth and reproduction
% [0 0 0 1 0]   costs for reproduction
% [0 0 0 0 1]   hazard for reproduction
%
% FEEDB: Feedbacks to use on damage dynamics as set of switches
% [surface:volume on uptake, surface:volume on elimination, growth dilution, losses with reproduction]
% [1 1 1 1]     all feedbacks
% [1 1 1 0]     classic DEBtox (no losses with repro)
% [0 0 1 0]     damage that is diluted by growth
% [0 0 0 0]     damage that is not diluted by growth

% These are the MoA's and feedback configurations that will be run
% automatically when fit_tox(1) = 1. This needs to be defined here, even
% for control fits when it is not used.

MOA   = [0 0 0 1 0];
FEEDB = [0 0 0 0;1 1 1 0]; % you may want to try more feedbacks ...
% MOA   = [0 0 0 1 0; 0 0 0 0 1]; % these are costs for repro and repro hazards
% FEEDB = allcomb([0 1],[0 1],[0 1],[0 1]); % or test ALL feedback configurations

% For fit_tox(1)~=1, the setting of MOA and FEEDB has no impact; they must
% be defined to prevent errors. For fit_tox(1)=1, setting is relevant. Note
% that if you run multiple configurations, the last one will remain in the
% memory, unless glo.moa and glo.feedb are redefined.

% Note: when using the parspace explorer, the automatic_runs uses the start
% ranges as produced by startgrid_debtox. These are preliminary! It may be
% needed to tweak these, but that should then be done in startgrid_debtox
% or in automatic_runs. So watch out when the analysis runs into min-max
% bounds.

% ===== FITTING CONTROLS ==================================================
% Simplex fitting works fine for control data
opt_optim.type = 1; % optimisation method: 1) default simplex, 4) parspace explorer

% opt_optim.type     = 4; % optimisation method 1) simplex, 4 parameter-space explorer
% opt_optim.ps_plots = 0; % when set to 1, makes intermediate plots of parameter space to monitor progress
% opt_optim.ps_rough = 1; % set to 1 for rough settings of parameter-space explorer, 0 for settings as in openGUTS (2 for extra rough)

% % Compare controls in data set
% fit_tox = [-2 1 3];
% automatic_runs_debtox2019(fit_tox,par,ind_tox,[],MOA,FEEDB,opt_optim,opt_plot);
% % script to run the calculations and plot, automatically
% % NOTE: I think the likelihood-ratio test is often too strict, and that
% % using both controls should be the default situation.

% Fit hb in data set
fit_tox = [-1 1 3]; % use both controls
par_out = automatic_runs_debtox2019(fit_tox,par,ind_tox,[],MOA,FEEDB,opt_optim,opt_plot);
% par_out = automatic_runs_debtox2019(fit_tox,par,ind_tox,[],MOA,FEEDB,opt_optim,opt_plot,opt_prof);
% script to run the calculations and plot, automatically
par = copy_par(par,par_out,1); % copy fitted parameters into par, and keep fit mark in par

% Fit other control parameters (not hb) in data set
fit_tox = [0 1 3]; % use both controls
par_out = automatic_runs_debtox2019(fit_tox,par,ind_tox,[],MOA,FEEDB,opt_optim,opt_plot);
% par_out = automatic_runs_debtox2019(fit_tox,par,ind_tox,[],MOA,FEEDB,opt_optim,opt_plot,opt_prof);
% script to run the calculations and plot, automatically
par = copy_par(par,par_out,1); % copy fitted parameters into par, and keep fit mark in par

% return

% ===== FITTING TOX DATA ==================================================
% Use the parameter-space explorer for fitting the treatments. Note that,
% with fit_tox=1, the tox parameters in par are replaced (when fitted) with
% estimates based on the data set using startgrid_debtox (called in
% automatic_runs_debtox2019).
%
% Note: you can use the rough settings to find better ranges and restart
% with refined settings. With opt_optim.ps_profs = 0, the algorithm will
% provide new search ranges on screen that can be directly copied-pasted
% into this script. You may comment out the fitting of the controls above.
% Make sure that skip_sg is then set to 1 to avoid startgrid_debtox to
% overwrite par again.

opt_optim.type     = 4; % optimisation method 1) simplex, 4 parameter-space explorer
% Note: to use saved set for parameter-space explorer, use fit_tox(2)=0!
opt_optim.ps_plots = 0; % when set to 1, makes intermediate plots of parameter space to monitor progress
opt_optim.ps_rough = 1; % set to 1 for rough settings of parameter-space explorer, 0 for settings as in openGUTS (2 for extra rough)
% -----------------------------------------------------------------------
% THIS BLOCK: SETTINGS FOR ROUGHLY GOING THROUGH MANY CONFIGURATIONS
% opt_optim.ps_rough = 2; % set to 1 for rough settings of parameter-space explorer, 0 for settings as in openGUTS (2 for extra rough)
opt_optim.ps_profs = 0; % when set to 1, makes profiles and additional sampling for parameter-space explorer
% glo.stiff          = [2 1]; % for quick exploration of many options, could try stiff solver with sloppy tolerances if default setting is slow/stuck
glo.stiff          = [0 3]; % use ode45 with very strict tolerances
skip_sg            = 0; % set to 1 to skip startgrid completely (use ranges in <par> structure)
% -----------------------------------------------------------------------
% % THIS BLOCK: SETTINGS TO REFINE FROM RANGES COPIES FROM SCREEN INTO THIS SCRIPT
% opt_optim.ps_profs = 1; % when set to 1, makes profiles and additional sampling for parameter-space explorer
% glo.stiff          = [0 3]; % use ode45 with very strict tolerances
% skip_sg            = 1; % set to 1 to skip startgrid completely (use ranges in <par> structure)
% -----------------------------------------------------------------------

fit_tox = [1 1 3]; % use both controls, and fit tox data
[par_out,best_MoaFb] = automatic_runs_debtox2019(fit_tox,par,ind_tox,skip_sg,MOA,FEEDB,opt_optim,opt_plot);
% script to run the calculations and plot, automatically
% Note: automatic_runs will return the BEST parameter set in par_out.
disp_settings_debtox2019 % display a bit more info on the settings on screen
print_par(par_out) % print the complete best parameter vector that can be copied-pasted
% This includes the control parameters (the parspace explorer itself will
% only plot the results for the fitted parameters).
Goodness-of-fit measures for each state and each data set (R-square)
Based on the individual replicates, accounting for transformations, and including t=0
state 1.    NaN
state 2.    NaN
state 3.    NaN
state 4.  -0.34
Warning: R-square is not appropriate for survival data so interpret
qualitatively! 
Warning: (R-square not calculated when missing/removed animals; use plot_tktd) 
 
=================================================================================
Results of the parameter estimation with BYOM (par-engine) version 6.0 BETA 7
=================================================================================
   Filename      : byom_debtox_daphnia_hb 
   Analysis date : 30-Nov-2021 (18:36) 
   Data entered  :
     data state 1: 0x0, no data.
     data state 2: 0x0, no data.
     data state 3: 0x0, no data.
     data state 4: 12x5, survival data.
   Search method: Nelder-Mead simplex direct search, 2 rounds. 
     The optimisation routine has converged to a solution
     Total 19 simplex iterations used to optimise. 
     Minus log-likelihood has reached the value 11.98 (AIC=25.96). 
=================================================================================
L0              1 (fit: 0, initial:  1) 
Lp            1.8 (fit: 0, initial: 1.8) 
Lm              3 (fit: 0, initial:  3) 
rB           0.15 (fit: 0, initial: 0.15) 
Rm             10 (fit: 0, initial: 10) 
f               1 (fit: 0, initial:  1) 
hb       0.004813 (fit: 1, initial: 0.01) 
Lf              0 (fit: 0, initial:  0) 
Lj              0 (fit: 0, initial:  0) 
Tlag            0 (fit: 0, initial:  0) 
kd           0.08 (fit: 0, initial: 0.08) 
zb            0.1 (fit: 0, initial: 0.1) 
bb             75 (fit: 0, initial: 75) 
zs           0.23 (fit: 0, initial: 0.23) 
bs            0.7 (fit: 0, initial: 0.7) 
=================================================================================
Parameter hb is fitted on log-scale.
=================================================================================
Mode of action used: 0  0  0  1  0
Feedbacks used     : 0  0  0  0
=================================================================================
Time required: 5.6 secs
Plots result from the optimised parameter values. 
par.hb     = [ 0.0048133  1      0.001       0.07 0]; 
Goodness-of-fit measures for each state and each data set (R-square)
Based on the individual replicates, accounting for transformations, and including t=0
state 1.    NaN
state 2.   0.95
state 3.   0.73
state 4.    NaN
=================================================================================
Results of the parameter estimation with BYOM (par-engine) version 6.0 BETA 7
=================================================================================
   Filename      : byom_debtox_daphnia_ctrl 
   Analysis date : 30-Nov-2021 (18:37) 
   Data entered  :
     data state 1: 0x0, no data.
     data state 2: 12x25, continuous data, no transf.
     data state 3: 12x50, continuous data, no transf.
     data state 4: 0x0, no data.
   Search method: Nelder-Mead simplex direct search, 2 rounds. 
     The optimisation routine has converged to a solution
     Total 184 simplex iterations used to optimise. 
     Minus log-likelihood has reached the value 1572.57 (AIC=3155.13). 
=================================================================================
L0         0.8696 (fit: 1, initial:  1) 
Lp          1.798 (fit: 1, initial: 1.8) 
Lm          3.113 (fit: 1, initial:  3) 
rB         0.1437 (fit: 1, initial: 0.15) 
Rm          10.02 (fit: 1, initial: 10) 
f               1 (fit: 0, initial:  1) 
hb       0.004813 (fit: 0, initial: 0.004813) 
Lf              0 (fit: 0, initial:  0) 
Lj              0 (fit: 0, initial:  0) 
Tlag            0 (fit: 0, initial:  0) 
kd           0.08 (fit: 0, initial: 0.08) 
zb            0.1 (fit: 0, initial: 0.1) 
bb             75 (fit: 0, initial: 75) 
zs           0.23 (fit: 0, initial: 0.23) 
bs            0.7 (fit: 0, initial: 0.7) 
=================================================================================
Mode of action used: 0  0  0  1  0
Feedbacks used     : 0  0  0  0
=================================================================================
Time required: 43.9 secs
Plots result from the optimised parameter values. 
par.L0     = [    0.8696  1        0.5        1.5 1]; 
par.Lp     = [    1.7976  1        0.5          4 1]; 
par.Lm     = [    3.1134  1          2          7 1]; 
par.rB     = [   0.14375  1       0.02        0.3 1]; 
par.Rm     = [    10.021  1          5         30 1]; 
Warning: The parameter-space explorer is not thoroughly tested for more than 4
free parameters! 
 
Settings for parameter search ranges:
=====================================================================
L0    fixed :     0.8696              fit: 0
Lp    fixed :      1.798              fit: 0
Lm    fixed :      3.113              fit: 0
rB    fixed :     0.1437              fit: 0
Rm    fixed :      10.02              fit: 0
f     fixed :          1              fit: 0
hb    fixed :   0.004813              fit: 0
Lf    fixed :          0              fit: 0
Lj    fixed :          0              fit: 0
Tlag  fixed :          0              fit: 0
kd    bounds:       0.01 -         10 fit: 1 (log)
zb    bounds:  0.0003547 -     0.8445 fit: 1 (norm)
bb    bounds:     0.5862 -  1.216e+04 fit: 1 (log)
zs    bounds:  0.0003547 -     0.8445 fit: 1 (norm)
bs    bounds:   0.005766 -      495.8 fit: 1 (log)
=====================================================================
Mode of action: 0  0  0  1  0
Feedbacks     : 0  0  0  0
Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 8).

ans = 

 ProcessPool with properties: 

            Connected: true
           NumWorkers: 8
              Cluster: local
        AttachedFiles: {}
    AutoAddClientPath: true
          IdleTimeout: 30 minutes (30 minutes remaining)
          SpmdEnabled: true

Using latin-hypercube sampling in first round, instead of regular grid.
 
Starting round 1 with initial grid of 10000 parameter sets
  Status: best fit so far is (minloglik) 4024.7335
Starting round 2, refining a selection of 400 parameter sets, with 25 tries each
  Status: 13 sets within total CI and 8 within inner. Best fit: 4021.8097
Starting round 3, refining a selection of 400 parameter sets, with 25 tries each
  Status: 53 sets within total CI and 17 within inner. Best fit: 4021.7478
Starting round 4, refining a selection of 400 parameter sets, with 25 tries each
  Status: 269 sets within total CI and 43 within inner. Best fit: 4021.7372
Starting round 5, refining a selection of 467 parameter sets, with 21 tries each
  Status: 1499 sets within total CI and 186 within inner. Best fit: 4021.7372
  Next round will focus on inner rim (outer rim has enough points)
Starting round 6, refining a selection of 400 parameter sets, with 10 tries each
  Status: 3030 sets within total CI and 523 within inner. Best fit: 4021.7372
  Finished sampling, running a simplex optimisation ...
  Status: 3030 sets within total CI and 524 within inner. Best fit: 4021.7364
 
Results from the limited parspace exploration (initial sampling only).
This can be copied-pasted into your BYOM script.
Warning: for oddly-shaped parameter spaces, this may miss local minima!
(in extreme cases, even the global optimum can be missed)
Suggestions for min-max search bounds returned in par_out.
=======================================================================
par.kd     = [  0.070606  1   0.025296    0.38538 0]; 
par.zb     = [     0.127  1   0.054317    0.24126 1]; 
par.bb     = [     74.33  1     11.924     1931.9 0]; 
par.zs     = [   0.29148  1   0.090385    0.84447 1]; 
par.bs     = [   0.48807  1    0.05987     3.1768 0]; 
=======================================================================
 
=================================================================================
Results of the parameter-space exploration
=================================================================================
   Sample: 3030 sets in joint CI and 524 in inner CI.
   Propagation set: 298 sets will be used for error propagation.
Best estimates and 95% CIs on fitted parameters
Rough settings were used.
  In almost all cases, this will be sufficient for reliable results.
Profiling and additional sampling steps were skipped.
  CIs below are just the edges of the inner cloud (underestimating the true CIs).
  Advised search ranges above are derived from approx. CIs +- 20%.
==========================================================================
kd    best:    0.07061 (   0.03366  -    0.1887 ) fit: 1 (log)
zb    best:     0.1270 (   0.07613  -    0.1829 ) fit: 1 (norm)
bb    best:      74.33 (     22.62  -     232.2 ) fit: 1 (log)
zs    best:     0.2915 (    0.1776  -    0.6494 ) fit: 1 (norm)
bs    best:     0.4881 (    0.1620  -     1.445 ) fit: 1 (log)
==========================================================================
 
==========================================================================
Depuration/repair time (DRT95)    : 42.4 (15.9 - 89.0) days 
Note that this depuration time does NOT account for any feedbacks from
life history traits back to TK/damage dynamics (in DEBtox analyses).
==========================================================================
 
Goodness-of-fit measures for each state and each data set (R-square)
Based on the individual replicates, accounting for transformations, and including t=0
state 1.    NaN
state 2.   0.95
state 3.   0.75
state 4.   0.92
Warning: R-square is not appropriate for survival data so interpret
qualitatively! 
Warning: (R-square not calculated when missing/removed animals; use plot_tktd) 
 
=================================================================================
Results of the parameter estimation with BYOM (par-engine) version 6.0 BETA 7
=================================================================================
   Filename      : byom_debtox_daphnia_moa00010_feedb0000 
   Analysis date : 30-Nov-2021 (19:29) 
   Data entered  :
     data state 1: 0x0, no data.
     data state 2: 12x25, continuous data, no transf.
     data state 3: 12x50, continuous data, no transf.
     data state 4: 12x5, survival data.
   Search method: Parameter-space explorer (see CIs above). 
     The optimisation routine has converged to a solution
     Minus log-likelihood has reached the value 4021.74 (AIC=8053.47). 
=================================================================================
L0         0.8696 (fit: 0, initial: NaN) 
Lp          1.798 (fit: 0, initial: NaN) 
Lm          3.113 (fit: 0, initial: NaN) 
rB         0.1437 (fit: 0, initial: NaN) 
Rm          10.02 (fit: 0, initial: NaN) 
f               1 (fit: 0, initial: NaN) 
hb       0.004813 (fit: 0, initial: NaN) 
Lf              0 (fit: 0, initial: NaN) 
Lj              0 (fit: 0, initial: NaN) 
Tlag            0 (fit: 0, initial: NaN) 
kd        0.07061 (fit: 1, initial: NaN) 
zb         0.1270 (fit: 1, initial: NaN) 
bb          74.33 (fit: 1, initial: NaN) 
zs         0.2915 (fit: 1, initial: NaN) 
bs         0.4881 (fit: 1, initial: NaN) 
=================================================================================
Parameters kd, bb and bs are fitted on log-scale.
=================================================================================
Mode of action used: 0  0  0  1  0
Feedbacks used     : 0  0  0  0
=================================================================================
Time required: 52 mins, 28.6 secs
Plots result from the optimised parameter values. 
Warning: The parameter-space explorer is not thoroughly tested for more than 4
free parameters! 
 
Settings for parameter search ranges:
=====================================================================
L0    fixed :     0.8696              fit: 0
Lp    fixed :      1.798              fit: 0
Lm    fixed :      3.113              fit: 0
rB    fixed :     0.1437              fit: 0
Rm    fixed :      10.02              fit: 0
f     fixed :          1              fit: 0
hb    fixed :   0.004813              fit: 0
Lf    fixed :          0              fit: 0
Lj    fixed :          0              fit: 0
Tlag  fixed :          0              fit: 0
kd    bounds:       0.01 -         10 fit: 1 (log)
zb    bounds:  0.0003547 -     0.8445 fit: 1 (norm)
bb    bounds:     0.5862 -  1.216e+04 fit: 1 (log)
zs    bounds:  0.0003547 -     0.8445 fit: 1 (norm)
bs    bounds:   0.005766 -      495.8 fit: 1 (log)
=====================================================================
Mode of action: 0  0  0  1  0
Feedbacks     : 1  1  1  0
Using latin-hypercube sampling in first round, instead of regular grid.
 
Starting round 1 with initial grid of 10000 parameter sets
  Status: best fit so far is (minloglik) 4023.1954
Starting round 2, refining a selection of 400 parameter sets, with 25 tries each
  Status: 11 sets within total CI and 7 within inner. Best fit: 4020.77
 
Slow kinetics indicated: parameter-space explorer is restarting with these settings:
=====================================================================
L0    fixed :     0.8696              fit: 0
Lp    fixed :      1.798              fit: 0
Lm    fixed :      3.113              fit: 0
rB    fixed :     0.1437              fit: 0
Rm    fixed :      10.02              fit: 0
f     fixed :          1              fit: 0
hb    fixed :   0.004813              fit: 0
Lf    fixed :          0              fit: 0
Lj    fixed :          0              fit: 0
Tlag  fixed :          0              fit: 0
kd    bounds:       0.01 -         10 fit: 1 (log)
zb    bounds:  0.0003547 -     0.8445 fit: 1 (log)
bb    bounds:     0.5862 -  1.216e+04 fit: 1 (log)
zs    bounds:  0.0003547 -     0.8445 fit: 1 (log)
bs    bounds:   0.005766 -      495.8 fit: 1 (log)
=====================================================================
Using latin-hypercube sampling in first round, instead of regular grid.
 
Starting round 1 with initial grid of 10000 parameter sets
  Status: best fit so far is (minloglik) 4026.2447
Starting round 2, refining a selection of 400 parameter sets, with 25 tries each
  Status: 9 sets within total CI and 7 within inner. Best fit: 4020.8345
Starting round 3, refining a selection of 400 parameter sets, with 25 tries each
  Status: 32 sets within total CI and 16 within inner. Best fit: 4020.753
Starting round 4, refining a selection of 400 parameter sets, with 25 tries each
  Status: 161 sets within total CI and 35 within inner. Best fit: 4020.7487
Starting round 5, refining a selection of 400 parameter sets, with 25 tries each
  Status: 1072 sets within total CI and 152 within inner. Best fit: 4020.7484
  Next round will focus on inner rim (outer rim has enough points)
Starting round 6, refining a selection of 400 parameter sets, with 10 tries each
  Status: 2337 sets within total CI and 443 within inner. Best fit: 4020.7483
  Next round will focus on inner rim (outer rim has enough points)
Starting round 7, refining a selection of 802 parameter sets, with 5 tries each
  Status: 4413 sets within total CI and 1030 within inner. Best fit: 4020.7481
  Finished sampling, running a simplex optimisation ...
  Status: 4414 sets within total CI and 1024 within inner. Best fit: 4020.745
 
Results from the limited parspace exploration (initial sampling only).
This can be copied-pasted into your BYOM script.
Warning: for oddly-shaped parameter spaces, this may miss local minima!
(in extreme cases, even the global optimum can be missed)
Suggestions for min-max search bounds returned in par_out.
=======================================================================
par.kd     = [    0.0312  1       0.01    0.38005 0]; 
par.zb     = [  0.071205  1   0.020773    0.21594 0]; 
par.bb     = [    98.173  1      18.54     5115.1 0]; 
par.zs     = [   0.17687  1   0.023931    0.65258 0]; 
par.bs     = [   0.91838  1    0.11909     7.2792 0]; 
=======================================================================
 
=================================================================================
Results of the parameter-space exploration
=================================================================================
   Sample: 4414 sets in joint CI and 1024 in inner CI.
   Propagation set: 550 sets will be used for error propagation.
Best estimates and 95% CIs on fitted parameters
Rough settings were used.
  In almost all cases, this will be sufficient for reliable results.
Profiling and additional sampling steps were skipped.
  CIs below are just the edges of the inner cloud (underestimating the true CIs).
  Advised search ranges above are derived from approx. CIs +- 20%.
==========================================================================
kd    best:    0.03120 (   0.01000  -   0.08518 ) fit: 1 (log)
zb    best:    0.07120 (   0.02606  -    0.1257 ) fit: 1 (log)
bb    best:      98.17 (     42.00  -     567.8 ) fit: 1 (log)
zs    best:     0.1769 (   0.05497  -    0.3678 ) fit: 1 (log)
bs    best:     0.9184 (    0.2812  -     4.243 ) fit: 1 (log)
==========================================================================
 
==========================================================================
Depuration/repair time (DRT95)    : 96.0 (35.2 - 300.) days 
Note that this depuration time does NOT account for any feedbacks from
life history traits back to TK/damage dynamics (in DEBtox analyses).
==========================================================================
 
Goodness-of-fit measures for each state and each data set (R-square)
Based on the individual replicates, accounting for transformations, and including t=0
state 1.    NaN
state 2.   0.95
state 3.   0.75
state 4.   0.92
Warning: R-square is not appropriate for survival data so interpret
qualitatively! 
Warning: (R-square not calculated when missing/removed animals; use plot_tktd) 
 
=================================================================================
Results of the parameter estimation with BYOM (par-engine) version 6.0 BETA 7
=================================================================================
   Filename      : byom_debtox_daphnia_moa00010_feedb1110 
   Analysis date : 30-Nov-2021 (21:13) 
   Data entered  :
     data state 1: 0x0, no data.
     data state 2: 12x25, continuous data, no transf.
     data state 3: 12x50, continuous data, no transf.
     data state 4: 12x5, survival data.
   Search method: Parameter-space explorer (see CIs above). 
     The optimisation routine has converged to a solution
     Minus log-likelihood has reached the value 4020.74 (AIC=8051.49). 
=================================================================================
L0         0.8696 (fit: 0, initial: NaN) 
Lp          1.798 (fit: 0, initial: NaN) 
Lm          3.113 (fit: 0, initial: NaN) 
rB         0.1437 (fit: 0, initial: NaN) 
Rm          10.02 (fit: 0, initial: NaN) 
f               1 (fit: 0, initial: NaN) 
hb       0.004813 (fit: 0, initial: NaN) 
Lf              0 (fit: 0, initial: NaN) 
Lj              0 (fit: 0, initial: NaN) 
Tlag            0 (fit: 0, initial: NaN) 
kd        0.03120 (fit: 1, initial: NaN) 
zb        0.07120 (fit: 1, initial: NaN) 
bb          98.17 (fit: 1, initial: NaN) 
zs         0.1769 (fit: 1, initial: NaN) 
bs         0.9184 (fit: 1, initial: NaN) 
=================================================================================
Parameters kd, zb, bb, zs and bs are fitted on log-scale.
=================================================================================
Mode of action used: 0  0  0  1  0
Feedbacks used     : 1  1  1  0
=================================================================================
Time required: 1 hour, 44 mins, 19.5 secs
Plots result from the optimised parameter values. 
 
MoA  feedbacks   MLL     delta-AIC     prob.    best
====================================================================
00010  0000    4021.74       1.98       0.37 
00010  1110    4020.74       0.00       1.00      *
====================================================================
 
========================================================================
Extra information on the model fit (DEBtox2019, stdDEB or extended GUTS)
Mode of action        : 0  0  0  1  0
Feedback configuration: 1  1  1  0
   these are only the settings for the LAST run
 
Model and ODE solver ----------------------------------------------
   Breaking time vector in call_deri: 0
   Settings for ODE solver call_deri: 0  3
   Brood-pouch delay for model (d)  : 3
Options for parspace explorer -------------------------------------
   Rough settings (recommended: 1)  : 1
   Make profile likelihoods         : 0
Options for startgrid ---------------------------------------------
   Skip automatic startgrid ranges  : 0
===================================================================
par.L0     = [    0.8696  0        0.5        1.5 1]; 
par.Lp     = [    1.7976  0        0.5          4 1]; 
par.Lm     = [    3.1134  0          2          7 1]; 
par.rB     = [   0.14375  0       0.02        0.3 1]; 
par.Rm     = [    10.021  0          5         30 1]; 
par.f      = [         1  0          0          2 1]; 
par.hb     = [ 0.0048133  0      0.001       0.07 0]; 
par.Lf     = [         0  0          0      1e+06 1]; 
par.Lj     = [         0  0          0      1e+06 1]; 
par.Tlag   = [         0  0          0      1e+06 1]; 
par.kd     = [    0.0312  1       0.01    0.38005 0]; 
par.zb     = [  0.071205  1   0.020773    0.21594 0]; 
par.bb     = [    98.173  1      18.54     5115.1 0]; 
par.zs     = [   0.17687  1   0.023931    0.65258 0]; 
par.bs     = [   0.91838  1    0.11909     7.2792 0]; 

Plot results with confidence intervals

The following code can be used to make plots with confidence intervals. Options for confidence bounds on model curves can be set using opt_conf (see prelim_checks). The plot_tktd function makes multiplots for the effects data, which are more readable when plotting with various intervals.

Note: the set up below plots the best fit only.

opt_conf.type    = 3; % make intervals from 1) slice sampler, 2)likelihood region, 3) parspace explorer
opt_conf.lim_set = 2; % use limited set of n_lim points (1) or outer hull (2, not for Bayes) to create CIs
opt_tktd.repls   = 0; % plot individual replicates (1) or means (0)
opt_tktd.transf  = 1; % set to 1 to calculate means and SEs including transformations
opt_tktd.obspred = 1; % plot predicted-observed plots (1) or not (0)
opt_tktd.max_exp = 0; % set to 1 to maximise exposure/damage plots on exposure rather than damage
opt_tktd.sppe    = 1; % set to 1 to calculate SPPEs (relative error at end of test)

% change X0mat to avoid plotting controls not used for fitting (assumes
% that regular control has identifier 0 and solvent control 0.1)
switch fit_tox(3)
    case 1 % remove solvent control
        X0mat(:,X0mat(1,:)==0.1) = [];
    case 2 % remove regular control
        X0mat(:,X0mat(1,:)==0) = [];
    case 3
        % leave all in
end

% Below some tricks to allow plotting results for selected configurations.
glo.moa    = MOA(best_MoaFb(1),:); % change global for MoA
glo.feedb  = FEEDB(best_MoaFb(2),:); % change global for feedback configuration
% glo.moa    = MOA(1,:);   % change global for MoA to first one
% glo.feedb  = FEEDB(1,:); % change global for feedback configuration to first one

glo.mat_nm = [basenm_rem,'_moa',sprintf('%d',glo.moa),'_feedb',sprintf('%d',glo.feedb)];
% The glo.mat_nm specifies the filename to read the parameters and sample
% from, as needed for plot_tktd.
glo.basenm = basenm_rem; % return the basenm to this filename
% The plots saved will get their filename based on the name of THIS
% script (so without the specification of MoA and feedbacks).

plot_tktd([],opt_tktd,opt_conf);
% leave the options opt_conf empty to suppress all CIs for these plots
% leave par_out (first input) empty to read it from saved mat file. If we
% ran more configurations, we should watch out when using par_out
% (automatic_run returns the best-fitting one).
 
Amount of 0.3764 added to the chi-square criterion for inner rim
Slightly more is taken to provide a generally conservative estimate of the CIs.
  Outer hull of 549 sets used.
Calculating CIs on model predictions, using sample from parspace explorer ... please be patient.
 
Plot_tktd uses complete parameter set from saved file.
 
Amount of 0.3764 added to the chi-square criterion for inner rim
Slightly more is taken to provide a generally conservative estimate of the CIs.
  Outer hull of 549 sets used.
Calculating CIs on model predictions, using sample from parspace explorer ... please be patient.
 
Time required: 1 hour, 45 mins, 3.3 secs
 
Goodness-of-fit metrics (treat with caution!)
 Values below are calculated on means, incl. t=0, from the residuals
 Residuals are used without accounting for any transformations
 However, means are calculated with transformation (if specified in DATA)
 
            Model efficiency (r2)   NRMSE
================================================================================
  survival    : 0.9154              0.0572 
  body length : 0.9754              0.0445 
  reproduction: 0.9878              0.1640 
================================================================================
 
Relative error percentage at the end of the test, per treatment (as in SPPE)
This is a somewhat curious goodness of fit measure, so interpret with care!
  - for survival it is 100x(observed prob-predicted prob), at end of test
  - for sub-lethal, it is 100x(observed-predicted)/observed, at end of test
  - for repro data, interpret with care when zero-observations have been removed
================================================================================
           Treatment    survival    length   reprod.
================================================================================
             control       9.6       7.9      -7.6
     solvent control     -10.4       1.2      -8.1
                  T1       9.6       1.8      -5.3
                  T2       4.8      -0.1      17.7
                  T3       2.8      -5.8      -Inf
================================================================================