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 can be modified by the user to decide whether to calculate the model results using the ODEs in derivatives.m, or the analytical solution in simplefun.m. Further, the events function can be turned on here, initial values can be determined by a parameter, and zero-variate data can be calculated.
useode = 1; % set to 1 solves the set of ODEs in derivatives, zero expects % an explicit function in simplefun eventson = 1; % set to 1 uses the events function: also modify the % sub-function at the bottom of this function! % Unpack the vector X0v, which is X0mat for one scenario X0 = X0v(2:end); % these are the intitial states for a scenario % % if needed, extract parameters from par that influence initial states in X0 L0 = par.L0(1); % initial body length (mm) X0(3) = L0; % put this parameter in the correct location of the initial vector % % if needed, calculate model values for zero-variate data from parameter set % if ~isempty(zvd) % zvd.ku(3) = par.Piw(1) * par.ke(1); % add model prediction as third value % end
Error using call_deri (line 46) Not enough input arguments.
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 TE = 0; % dummy for time of events if useode == 1 % use the ODE solver to calculate the solution % call the ODE solver if isempty(options.Events) % if no events function is specified ... [tn, Xout] = ode45(@derivatives, t, X0, options, par, c); else % with an events functions ... [tn, Xout, TE, YE, IE] = ode45(@derivatives, t, X0, options, par, c); % 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 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.
% Xout(:,1) = Xout(:,1).^3; % e.g., do something on first column, like cube it ... % % To obtain the output of the derivatives at each time point. The values in % % dXout might be used to replace values in Xout, if the data to be fitted % % are the changes (rates) instead of the state variable itself. % % dXout = zeros(size(Xout)); % initialise with zeros % for i = 1:length(t) % run through all time points % dXout(i,:) = derivatives(t(i),Xout(i,:),par,c); % % derivatives for each stage at each time % end
Modify this part of the code if eventson=1. This subfunction catches the 'events': in this case, it looks for the scaled internal concentration where the NEC is exceeded, and the start of investment in reproduction (L=Lp). This function should be adapted to the problem you are modelling (this one matches the byom_debtox_... files). You can catch more events by making a vector out of values.
Note that the eventsfun has the same inputs, in the same sequence, as derivatives.m.
function [value,isterminal,direction] = eventsfun(t,X,par,c) Lp = par.Lp(1); % length at puberty (mm) c0 = par.c0(1); % no-effect concentration sub-lethal c0s = par.c0s(1); % no-effect concentration survival nevents = 3; % number of events that we try to catch value = zeros(nevents,1); % initialise with zeros value(1) = X(1) - c0; % follow when scaled internal conc exceed the NEC value(2) = X(3) - Lp; % follow when body length exceeds length at puberty value(3) = X(1) - c0s; % follow when scaled internal conc exceed the NEC 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