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:
- t is the time point, provided by the ODE solver
- X is a vector with the previous value of the states
- par is the parameter structure
- c is the external concentration (or scenario number)
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.
function dX = derivatives(t,X,par,c)
global glo % allow for global parameters in structure glo
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
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