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.m, or the explicit function(s) in 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.m or simplefun.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.

This function is for the reduced GUTS model. The external function derivatives or simplefun provides the scaled damage over time and the background mortality. Mortality due to chemical stress is calculated in this function as a form of 'output mapping'.

Author: Tjalling Jager
Date: September 2018
Web support: http://www.debtox.info/byom.html
Back to index walkthrough_guts.html

Start
Initial settings
Calculations
Output mapping
Events function

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.m, or the analytical solution in 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 or 2 to use a stiff solver instead of the standard one
min_t    = 500;          % minimum length of time vector (affects ODE stepsize as well)

% 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.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 ode113 or ode15s instead.

c       = X0v(1); % the concentration (or scenario number)
t       = t(:);   % force t to be a row vector (needed when useode=0)
t_rem   = t;      % remember the original time vector (as we will add to it)

options = odeset; % start with default options for the ODE solver
% options = odeset(options, 'RelTol',1e-5,'AbsTol',1e-8); % specify tightened tolerances

if isfield(glo,'int_scen') && ismember(c,glo.int_scen) % is c in the scenario range global?
    % then we have a time-varying concentration
    Tev = make_scen(-2,c,-1); % use make_scen again to derive actual exposure concentration
    % the -2 lets make_scen know we need events, the -1 that this is also needed for splines!
    min_t = max(min_t,length(Tev(Tev(:,1)<t(end),1))*2);
    % For very long exposure profiles (e.g., FOCUS profiles), we now
    % automatically generate a larger time vector (twice the number of the
    % relevant points in the scenario).
    options = odeset(options,'InitialStep',t_rem(end)/(10*min_t),'MaxStep',t_rem(end)/min_t);
    % For the ODE solver, when we have a time-varying exposure set here, we
    % base minimum step size on min_t. Small stepsize is a good idea for
    % pulsed exposures; otherwise, stepsize may become so large that a
    % concentration change is missed completely.
else
    % options = odeset(options,'InitialStep',max(t)/100,'MaxStep',max(t)/10); % specify smaller stepsize
    % for constant concentrations, we can use a default or we can remove this option alltogether!
    % Matlab uses as default the length of the time vector divided by 10.
end
% For the GUTS cases, we need a long time vector. For the hazard and proper
% model, we need to numerically integrate the hazard rate over time. For IT
% in reduced models, we could make shortcuts, but since we don't know how
% people may modify the models in derivatives or simplefun, it is better to
% be general here. If the time span of the data set or simulation is very
% long (and the number of points in Tev is not), consider increasing the
% size of min_t.
if length(t) < min_t % make sure there are at least min_t points
    t = unique([t;(linspace(t(1),t(end),min_t))']);
end

TE = 0; % dummy for time of events
if useode == 1 % use the ODE solver to calculate the solution
    if eventson == 0 % if no events function is used ...
        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
        options = odeset(options, 'Events',@eventsfun); % add an events function
        switch stiff % note that YE and IE can be added as additional outputs
            case 0
                [~,Xout,TE,~,~] = ode45(@derivatives,t,X0,options,par,c);
            case 1
                [~,Xout,TE,~,~] = ode113(@derivatives,t,X0,options,par,c);
            case 2
                [~,Xout,TE,~,~] = ode15s(@derivatives,t,X0,options,par,c);
        end
    end
    Xout = max(0,Xout); % in extreme cases, states can become ever so slightly negative

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.

This set of files is geared towards the use of the analytical solution for the damage level in simplefun (although it can also be used with the ODEs in derivatives). Therefore, only the damage level and background survival is available at this point, and the survival due to the chemical is added for each death mechanism below.

mw   = par.mw(1);             % median of threshold distribution
bw   = par.bw(1);             % killing rate
Fs   = max(1+1e-6,par.Fs(1)); % fraction spread of the threshold distribution
beta = log(39)/log(Fs);       % shape parameter for logistic from Fs
S    = Xout(:,glo.locS);      % take background survival from the model
Dw   = Xout(:,glo.locD);      % take the correct state variable for scaled damage

switch glo.sel
    case 1 % stochastic death
        haz    = bw * max(0,Dw-mw);          % calculate hazard for each time point
        cumhaz = cumtrapz(t,haz);            % integrate the hazard rate numerically
        S      = S .* min(1,exp(-1*cumhaz)); % calculate survival probability, incl. background

    case 2 % individual tolerance, log-logistic distribution is used
        mw = max(mw,1e-100); % make sure that the threshold is not exactly zero ...
        % New method to make sure that Dw does not decrease over time (dead
        % animals dont become alive). This increases speed, especially when
        % there is very little decrease in Dw.
        maxDw = Dw; % copy the vector Dw to maxDw
        ind   = find([0;diff(maxDw)]<0,1,'first'); % first index to places where Dw has decreased
        while ~isempty(ind) % as long as there is a decrease somewhere ...
            maxDw(ind:end) = max(maxDw(ind:end),maxDw(ind-1)); % replace every later time with max of that and previous point
            ind = find([0;diff(maxDw)]<0,1,'first'); % any decrease left?
        end
        S = S .* (1 ./ (1+(maxDw/mw).^beta)); % survival probability
        % the survival due to the chemical is multiplied with the background survival

    case 3 % mixed model
        % This is a fast way for GUTS proper. Damage is calculated only
        % once, as it is the same for all individuals. Survival for
        % different NECs is calculated below.
        n          = 200; % number of slices from the threshold distribution
        Fs2        = 999^(1/beta); % fraction spread for 99.9% of the distribution
        z_range    = linspace(mw/(1.5*Fs2),mw*Fs2,n); % range of NECs to cover 99.9%
        prob_range = ((beta/mw)*(z_range/mw).^(beta-1)) ./ ((1+(z_range/mw).^beta).^2); % pdf for the log-logistic (Wikipedia)
        prob_range = prob_range / sum(prob_range); % normalise the densities to exactly one
        S1         = zeros(length(t),1); % initialise the survival probability over time with zeros
        for i = 1:n % run through the different thresholds
            haz    = bw * max(0,Dw-z_range(i));  % calculate hazard for each NEC
            cumhaz = cumtrapz(t,haz);            % integrate the hazard rate numerically
            surv   = min(1,exp(-1*cumhaz));      % calculate survival probability
            S1     = S1 + surv * prob_range(i);  % add to S1, weighted for this NECs prob density
        end
        S = S .* S1; % make sure to add background hazard
end

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

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