BYOM function derivatives.m (the model in ODEs)

Syntax: dX = derivatives(t,X,par,c)

This function calculates the derivatives for the model system. It is linked to the script files for the snails in the DEBkiss1 package. As input, it gets:

Time t and scenario name c are handed over as single numbers by call_deri.m (you do not have to use them in this function). Output dX (as vector) provides the differentials for each state at t.

Contents

Start

function dX = derivatives(t,X,par,c)
global glo          % allow for global parameters in structure glo

Unpack states

The state variables enter this function in the vector X. Here, we give them a more handy name.

WV = X(1); % state 1 is the structural body mass
cR = X(2); % state 2 is the cumulative reproduction (not used here)
WB = X(3); % state 3 is the egg buffer of assimilates

WV = max(0,WV); % for poor choice of paramaters, WV can become nagtive, which smothers the ODE solver
Not enough input arguments.

Error in derivatives (line 34)
WV = X(1); % state 1 is the structural body mass

Unpack parameters

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.

Tlag  = glo.Tlag;   % lag time for start embryo development
dV    = glo.dV;     % dry weight density of structure

if c==1
    Tlag = 0; % the c=1 scenario for embryos does not include a lag time
    % for the juvenile/adult runs, this does not hurt as Tlag=0 anyway
end

sJAm  = par.sJAm(1);  % specific assimilation rate
sJM   = par.sJM(1);   % specific maintenance costs
WB0   = par.WB0(1);   % initial weight of egg
WVp   = par.WVp(1);   % body mass at puberty
yAV   = par.yAV(1);   % yield of assimilates on volume (starvation)
yBA   = par.yBA(1);   % yield of egg buffer on assimilates
yVA   = par.yVA(1);   % yield of structure on assimilates (growth)
kap   = par.kap(1);   % allocation fraction to soma
f1    = par.f1(1);    % scaled food level, ad libitum
f2    = par.f2(1);    % scaled food level, regime 2
f3    = par.f3(1);    % scaled food level, regime 3
fB    = par.fB(1);    % scaled food level for the embryo (adapted model)
WVf   = par.WVf(1);   % half-saturation size for initial food limitation

Calculate the derivatives

This is the actual model, specified as a system of ODEs. Note: some unused fluxes are calculated because they may be needed for future modules (e.g., respiration and ageing).

L = (WV/dV)^(1/3); % volumetric length from dry weight

% select what to do with maturity maintenance
if glo.mat == 1
    sJJ = sJM * (1-kap)/kap; % add specific maturity maintenance with the suggested value
else
    sJJ = 0; % or ignore it
end

WVb = WB0 *yVA * kap; % predicted body mass at birth

if WB > 0 % if we have an embryo (run from byom_snail_embryo.m)
    if c==2 % for the 'adapted' scenario ...
        f = fB; % assimilation at limited rate
    else % for the 'standard' scenario
        f = 1;  % assimilation at maximum rate
    end
else % we have jubeniles/adults (run from byom_snail_fit.m or byom_snail_starv.m)
    switch c % run through the three feeding regimes to select proper f
        case 1
            f = f1;
        case 2
            f = f2;
            yAV = 7*0.8; % this is used for the adapted starvation simulation
            % no problem if it is also on for the fits, as there is no
            % starvation occurring there anyway!
        case 3
            f = f3;
    end
    if WVf > 0
        f = f / (1+WVf/WV); % hyperbolic relationship for f with body weight
    end
end

% if WV < WVb && WB <= 0
%     error('Your animal is initiated or has shrunk to less than the size at birth! Increase the start weight or check the food levels.')
% end

JA = f * sJAm * L^2;          % assimilation
JM = sJM * L^3;               % somatic maintenance
JV = yVA * (kap*JA-JM);       % growth

if WV < WVp                   % below size at puberty
    JR = 0;                   % no reproduction
    JJ = sJJ * L^3;           % maturity maintenance flux
    % JH = (1-kap) * JA - JJ; % maturation flux (not used!)
else
    JJ = sJJ * (WVp/dV);      % maturity maintenance
    JR = (1-kap) * JA - JJ;   % reproduction flux
end

% starvation rules may override these fluxes
if kap * JA < JM      % allocated flux to soma cannot pay maintenance
    if JA >= JM + JJ  % but still enough total assimilates to pay both maintenances
        JV = 0;       % stop growth
        if WV >= WVp  % for adults ...
            JR = JA - JM - JJ; % repro buffer gets what's left
        else
            % JH = JA - JM - JJ; % maturation gets what's left (not used!)
        end
    elseif JA >= JM   % only enough to pay somatic maintenance
        JV = 0;       % stop growth
        JR = 0;       % stop reproduction
        % JH = 0;     % stop maturation for juveniles (not used)
        JJ = JA - JM; % maturity maintenance flux gets what's left
    else              % we need to shrink
        JR = 0;       % stop reproduction
        JJ = 0;       % stop paying maturity maintenance
        % JH = 0;     % stop maturation for juveniles (not used)
        JV = (JA - JM) / yAV; % shrink; pay somatic maintenance from structure
    end
end

% we do not work with a repro buffer here, so a negative JR does not make
% sense; therefore, minimise to zero
JR = max(0,JR);

% Calcululate the derivatives
dWV = JV;             % change in body mass
dcR = yBA * JR / WB0; % continuous reproduction flux
if WB > 0             % for embryos ...
    dWB = -JA;        % decrease in egg buffer
else
    dWB = 0;          % for juveniles/adults, no change
end
if WV < WVb/4 && dWB == 0 % dont shrink below quarter size at birth
    dWV = 0; % simply stop shrinking ...
end

if t < Tlag
    dX = [0;0;0];
else
    dX = [dWV;dcR;dWB];       % collect both derivatives in one vector
end