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 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'.
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 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
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
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
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); % threshold for effects 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 dose metric 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