BYOM function simplefun.m (the model as explicit equations)
Syntax: Xout = simplefun(t,X0,par,c)
This function calculates the output of the reduced GUTS model system. Note that the survival probability due to chemical stress is all calculated in call_deri.m. This version can deal with time-varying exposure concentrations, as long as the exposure profile is specified using make_scen as type 2, 3 or 4 (smooth splining requires the ODE version of the model in derivatives.m). If type 1 (interpolation) is used for the exposure profile, make_scen will produce an error. I now also added a make_scen type 4, which applies the analytical solution for scaled damage under linear forcing sequentially (many thanks to Bob Kooi). This allows to use an analytical solution for more complex exposure scenarios as well.
NOTE: watch out when the rate constant (kd and kc) get close to each other. In the analytical solution, there are divisions by the difference between two rate constants, so they are not allowed to be equal. And, there will be numerical problems when they get too close. I try to catch these cases, and modify one or more of the rate constants a tiny bit. However, it is good to be aware of this limitation (when in doubt: also run the ODE version).
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)
Unpack initial states
The state variables enter this function in the vector X0. The initial scaled damage level does not have to be zero. However, a non-zero value would generally be meaningless (look at the full model for examples with non-zero initial body residue). Furthermore, non-zero Dw0 will likely lead to erroneous results for the calculation of LCx and LPx.
% S0 = X0(1); % survival probability at t=0 Dw0 = X0(2); % scaled damage (referenced to external concentrations) 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 % mw = par.mw(1); % median threshold (used in call_deri) % bw = par.bw(1); % killing rate (used in call_deri) % Fs = par.Fs(1); % fraction spread of NEC distribution (used in call_deri) hb = par.hb(1); % background hazard rate
Calculate the model output
This is the actual model, specified as explicit function(s):
S = exp(-hb*t); % fill S with background mortality for survival % The calculation of the actual survival probability is done in call_deri % and only background survival is calculated here. % Do we need to do fast or slow kinetics? if isfield(glo,'fastslow') % does that field exist? than we may need to do something else if glo.fastslow ~= 'o' % if it is not 'off' ... c_v = c * ones(length(t),1); % if no exposure profile is specified, simply copy c over all time points if isfield(glo,'int_scen') % if it exists: use it to derive current external conc. if ismember(c,glo.int_scen) % is c in the scenario range global? c_v = make_scen(-3,c,t); % use make_scen again to derive actual exposure concentration vector % the -3 lets make_scen know we are calling for fast/slow kinetics (and need a conc. vector) end end if glo.fastslow == 'f' Dw = c_v; % fast kinetics: damage equals external concentration Dw(1) = Dw0; % make initial one Dw0 (especially needed for IT) else Dw = cumtrapz(t,c_v); % for slow kinetics, integrate over t end Xout = [S Dw]; % combine them into a matrix return % return as we are done here end end % If regular kinetics ... Tev = [0 c]; % events setting: without anything else, assume it is constant kc = 0; % assume no disappearance of the test chemical if isfield(glo,'int_scen') % if it exists: use it to derive current external conc. if ismember(c,glo.int_scen) % is c in the scenario range global? [Tev,kc] = make_scen(-2,c,t); % use make_scen again to derive actual exposure concentration % the -2 lets make_scen know we are calling from simplefun (and need events and kc) end end if sum(Tev(:,2)) == 0 && Dw0 == 0 % no need to calculate anything if there is only zero concentrations in Tev Dw = zeros(length(t),1); Xout = [S Dw]; % combine them into a matrix return % return as we are done here end % initialise internal concentrations and damage with NaNs Dw = nan(length(t),1); Dw(1) = Dw0; % the first element is the starting concentration diff_rel = 1e-5; % minimum relative difference between the rate constants if size(Tev,2) == 2 % than we have a pulsed or static renewal scenario if abs(1-kd/kc) < diff_rel % the analytical solution does not allow the two rate constants % to be exactly the same (or too close) ... when kc=0, the abs % gives inf, so the part of the code below is not run. kd = kd * (1+diff_rel); % increase kd tiny bit end for i = 1:size(Tev,1) % run through all event periods a = kd * 1 * Tev(i,2); % modify the uptake flux to the new start concentration in this period % compound parameter a is used in the solution below if i<size(Tev,1) ind_t = (t>Tev(i,1) & t<=Tev(i+1,1)); % find the logical indices for the period else ind_t = (t>Tev(i,1)); % find the logical indices for the last period % if the scenario is longer than the data, this is all zeros, % which leads to an empty te, but this does not produce an % error so it is fine. end te = t(ind_t) - Tev(i,1); % find the new part of the time vector, and make it start at zero again % analytical solution Dw(ind_t) = (a/(kd-kc))*exp(-kc * te) + (Dw0 - (a/(kd-kc))) * exp(-kd * te); % find new starting values at exact moment of new period if i < size(Tev,1) % don't do this for last period te = Tev(i+1,1) - Tev(i,1); % start time for new period Dw0 = (a/(kd-kc))*exp(-kc * te) + (Dw0 - (a/(kd-kc))) * exp(-kd * te); end end else % we are using linear interpolation in a forcing series if t(end) > Tev(end,1) % do we ask for more points than in Tev? Tev = [Tev; t(end) 0 0]; % add a dummy time point in Tev (from t) end for i = 1:size(Tev,1)-1 % run through all events (not last one that we added) ind_t = (t>Tev(i,1) & t<=Tev(i+1,1)); % find the logical indices for the period te = t(ind_t) - Tev(i,1); % find the new part of the time vector, and make it start at zero again % Thanks to Bob Kooi for using Maple to obtain the solution a = Tev(i,2); % take as initial concentration the new one in this period b = Tev(i,3); % take as slope the new one in this period Dw(ind_t) = b * te + a - b / kd + exp(-kd * te) * (Dw0 - a + b / kd); % find new starting values at exact moment of new period if i < size(Tev,1)-1 % don't do this for last period te = Tev(i+1,1) - Tev(i,1); % start time for new period Dw0 = b * te + a - b / kd + exp(-kd * te) * (Dw0 - a + b / kd); end end end Xout = [S Dw]; % combine them into a matrix