BYOM function call_deri.m (calculates the model output)

Syntax: [Xout TE] = call_deri(t,par,X0v)

This function calls the ODE solver to solve the system of differential equations specified in derivatives.m, or the explicit function(s) in simplefun.m. It is linked to the model in the directory 'compound_model'. As input, it gets:

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.

Copyright (c) 2012-2019, Tjalling Jager, all rights reserved.
This source code is licensed under the MIT-style license found in the
LICENSE.txt file in the root directory of BYOM.



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

Initial settings

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

% if needed, extract parameters from par that influence initial states in X0
% start from specified initial size in a model parameter
L0M = par.L0M(1); % initial body length (mm)
X0(2) = glo.dV*((L0M*glo.delM)^3); % recalculate to dwt, and put this parameter in the correct location of the initial vector

% % start from the value provided in X0mat
% if glo.len ~= 0 % only if the switch is set to 1 or 2!
%     % Assume that X0mat includes initial length; if it includes mass, no correction is needed here
%     % The model is set up in masses, so we need to translate the initial length
%     % (in mm body length) to body mass using delM. Assume that size is second state.
%     X0(2) = glo.dV*((X0(2) * glo.delM)^3); % initial length to initial dwt
% end

% % start at hatching, so calculate initial length from egg weight
% this is possible here as well, but requires a lot of recalculations
% (which can be copied from derivatives)

% % if needed, calculate model values for zero-variate data from parameter set
% if ~isempty(zvd)
%     zvd.LmM(3) = (par.sJAm(1) * par.kap(1)/par.sJM(1))/glo.delM; % add model prediction as third value
% end


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
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 (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);
    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);
else % alternatively, use an explicit function provided in simplefun!
    Xout = simplefun(t,X0,par,c);

if isempty(TE) || all(TE == 0) % if there is no event caught
    TE = +inf; % return infinity

Output mapping

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.

% here, translate from the state body weight to output in length, if needed
if glo.len ~= 0 % only if the switch is set to 1 or 2!
    % Here, I translate the model output in body mass to estimated physical body length
    W = Xout(:,2); % assume that dry body weight is the second state
    LM = (W/glo.dV).^(1/3)/glo.delM; % estimated physical body length
    if glo.len == 2 % when animal cannot shrink in length (but does on volume!)
        maxLM = 0;%zeros(size(Lw)); % remember max size achieved
        for i = 1:length(LM) % run through length data
            maxLM = max([maxLM;LM(i)]); % new max is max of old max and previous size
            LM(i) = maxLM; % replace value in output vector
    Xout(:,2) = LM; % replace body weight by physical body length

% % 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

Events function

Modify this part of the code if eventson=1. This subfunction catches the 'events': in this case, this one catches when the scaled internal concentration exceeds one of the NECs, and catches the switch at puberty.

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
LpM = par.LpM(1); % length at puberty (mm)
zb  = par.zb(1);  % effect threshold for the energy budget
zs  = par.zs(1);  % effect threshold for survival

WVp = glo.dV * ((LpM * glo.delM)^3); % translate physical length to dry weight

nevents = 3; % number of events that we try to catch

value    = zeros(nevents,1); % initialise with zeros
value(1) = X(1) - zb;        % follow when scaled damage exceed the threshold for energy-budget effects
value(2) = X(1) - zs;        % follow when scaled damage exceed the threshold for survival
value(3) = X(2) - WVp;       % follow when body weight exceeds length at puberty

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