BYOM function call_deri.m (calculates the model output)
Syntax: [Xout TE] = call_deri(t,par,X0v)
- 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. 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.
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
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 to use a stiff solver instead of the standard one % Unpack the vector X0v, which is X0mat for one scenario X0 = X0v(2:end); % these are the intitial states for a scenario
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 ode15s instead.
c = X0v(1); % the concentration (or scenario number) t = t(:); % force t to be a row vector (needed when useode=0) % specify options for the ODE solver options = odeset; % start with default options if eventson == 1 options = odeset(options, 'Events',@eventsfun); % add an events function end options = odeset(options, 'RelTol',1e-4,'AbsTol',1e-7); % specify tightened tolerances options = odeset(options,'InitialStep',max(t)/1000,'MaxStep',max(t)/100); % specify smaller stepsize % smaller stepsize is a good idea for pulsed exposures; otherwise, stepsize % may become so large that a pulse is missed completely TE = 0; % dummy for time of events % For the GUTS packages, we need a long time vector for IT and for 'proper % GUTS' (SD and IT combined). For the IT this is needed when damage can % also decrease in time: we have to find the maximum of the damage over % time. For the proper model, we need to numerically integrate the hazard % rate over time here. If the time span of the data set or simulation is % very long, consider increasing the size of t below. t_rem = t; % remember the original time vector if glo.sel ~= 1 && length(t)<100 t = unique([t;(linspace(min(t),max(t),100))']); end if useode == 1 % use the ODE solver to calculate the solution % call the ODE solver (use ode15s for stiff problems, escpecially for pulses) if isempty(options.Events) % if no events function is specified ... 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 switch stiff case 0 [~,Xout,TE,YE,IE] = ode45(@derivatives,t,X0,options,par,c); case 1 [~,Xout,TE,YE,IE] = ode113(@derivatives,t,X0,options,par,c); case 2 [~,Xout,TE,YE,IE] = ode15s(@derivatives,t,X0,options,par,c); end end 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
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.
Stochastic death is calculated in simplefun.m or in derivatives.m directly. For the individual tolerance calculation, we need a few extra things here.
% Note that for IT, a long time vector will be generated above, which is % needed to catch the maximum dose metric over time during pulsed exposure. if glo.sel == 2 % for individual tolerance ... mw = par.mw(1); % median of thhreshold distribution, ug/L Fs = par.Fs(1); % fraction spread of the threshold distribution Dw = Xout(:,glo.locD); % take the correct state variable as scaled damage % Dw(Dw<=0) = 1e-20; % take care scaled damage <= 0 (probably not needed) maxDw = zeros(size(Dw)); % intitialise a new vector for maximum Dw over time for j = 1:length(t), % this section makes sure that Dw does not decrease in time (dead animals dont become alive) % Note: for IT, the time vector will always be long (see above) maxDw(j,:) = max(Dw(1:j,:)); % so look for max of Dw over time end % Log-logistic distribution is used for GUTS-IT beta = log(39)/log(Fs); % translate Fs to beta S = (1 ./ (1+(maxDw/mw).^beta)) .* Xout(:,glo.locS); % survival probability % the survival due to the chemical is multiplied with the background % survival (calculated in derivatives), assuming logistic distrib. 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 end
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.
function [value,isterminal,direction] = eventsfun(t,X,par,c) global glo mw = par.mw(1); % median for thresholds nevents = 1; % number of events that we try to catch value = zeros(nevents,1); % initialise with zeros value(1) = X(glo.locD) - mw; % thing to follow is damage minus threshold isterminal = zeros(nevents,1); % do NOT stop the solver at an event direction = zeros(nevents,1); % catch ALL zero crossing when function is increasing or decreasing