BYOM function simplefun.m (the model as explicit equations)

Syntax: Xout = simplefun(t,X0,par,c)

This function calculates the output of the model system. It is linked to the script files byom_guts1.m. As input, it gets:

Time t is handed over as a vector, and scenario name c as single number, by call_deri.m (you do not have to use them in this function). Output Xout (as matrix) provides the output for each state at each t.



function Xout = simplefun(t,X0,par,c)
global glo   % allow for global parameters in structure glo (handy for switches)

% This calculation of survival only works when the exposure concentration
% is constant over time! This is the reduced GUTS model.

Unpack initial states

The state variables enter this function in the vector _X_0. However, the analytical solution is (for now) limited to the situation where the initial scaled damage level is zero.

% S0  = X0(1); % survival probability at t=0
% Dw0 = X0(2); % scaled damage (referenced to water) at t=0

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   =;   % median threshold, ug/L
bw   =;   % 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 model output

This is the actual model, specified as explicit function(s):

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

S  = exp(-hb*t); % initialise S with background mortality for survival
% Dw = Dw0*exp(-kd*t)+(1-exp(-kd*t)) * c; % scaled damage over time
Dw = (1-exp(-kd*t)) * c; % scaled damage over time
% Note: the initial damage is excluded as the analytical solution assumes Dw0=0!

if glo.sel == 1 % only for stochastic death model

    % This is the same calculation as in DEBtool fomort
    t0 = -log(max(1e-8,1-mw/max(1e-8,c)))/kd; % no-effect-time
    t1 = ones(length(t),1);                   % column-matrix of ones
    f  = (bw/kd)*max(0,t1*exp(-kd*t0) - exp(-kd*t)).*(t1*c) - ...
        bw*(t1*(max(0,c-mw)')).*max(0,t - t1*t0);

    S = min(1,exp(f)) .* S; % multiply with background mortality


% for individual tolerance, the calculation is done in call_deri and only
% background survival is calculated here

Xout = [S Dw]; % combine them into a matrix