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:
- t is the time vector
- X0 is a vector with the initial values for states
- par is the parameter structure
- c is the external concentration (or scenario number)
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
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 threshold, 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 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 end % 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