%% get_pars_9 % Obtains 9 DEB parameters from 9 data points at abundant food %% function par = get_pars_9(data, fixed_par, chem_par) % created 2013/07/08 by Bas Kooijman; modified 2015/01/17, 2020/12/24 %% Syntax % par = <../get_pars_9.m *get_pars_9*> (data, fixed_par, chem_par) %% Description % Obtains 9 DEB parameters from 9 data points at abundant food. % Accounts for type M acceleration % % Input % % data: 9-vector with zero-variate data % % d_V: g/cm^3 specific density of structure % a_b: d, age at birth % a_p: d, age at puberty % a_m: d, age at death due to ageing % W_b: g, wet weight at birth % W_j: g, wet weight at metamorphosis % W_p: g, wet weight at puberty % W_m: g, maximum wet weight % R_m: #/d, maximum reproduction rate % % * fixed_par: optional vector with k_J, s_G, kap_R, kap_G % * chem_par: optional 4 vector with w_V, w_E, mu_V, mu_E % % Output % % * par: 9-vector with DEB parameters % % p_Am: J/d.cm^2, {p_Am}, max specific assimilation rate % v: cm/d, energy conductance % kap: -, allocation fraction to soma % p_M: J/d.cm^3, [p_M], specific somatic maintenance costs % E_G: J/cm^3, [E_G] specific cost for structure % E_Hb: J, E_H^b, maturity at birth % E_Hj: J, E_H^j, maturity at metamorphosis % E_Hp: J, E_H^p, maturity at puberty % h_a: 1/d^2, ageing acceleration %% Remarks % The theory behind this mapping is discussed in % . % See for the inverse mapping and % , % , % , % , % , % , % , % , % . %% Example of use % See <../mydata_get_par_9.m *mydata_get_par_9*> % assumptions: % abundant food (f = 1) % a_b, a_p, a_m and R_m are temp-corrected to T_ref = 293 K % p_T = 0; % J/d.cm^2, {p_T}, surf-spec som maint if exist('fixed_par','var') == 0 || isempty(fixed_par) k_J = 0.002; % 1/d, maturity maintenance rate coefficient s_G = 1e-4; % -, Gompertz stress coefficient (= small) kap_R = 0.95; % -, reproduction efficiency kap_G = 0.80; % -, growth efficiency else k_J = fixed_par(1); % 1/d, maturity maintenance rate coefficient s_G = fixed_par(2); % -, Gopertz stress coefficient (= small) kap_R = fixed_par(3); % -, reproduction efficiency kap_G = fixed_par(4); % -, growth efficiency end if exist('chem_par', 'var') == 0 || isempty(chem_par) % C:H:O:N = 1:1.8:0.5:0.15 w_V = 23.9; % g/C-mol, molecular weight of structure w_E = 23.9; % g/C-mol, molecular weight of reserve mu_V = 5E5; % J/C-mol, chemical potential of structure mu_E = 5.5E5; % J/C-mol, chemical potential of reserve else w_V = chem_par(1); w_E = chem_par(2); mu_V = chem_par(3); mu_E = chem_par(4); end % unpack data d_V = data(1); % g/cm^3 specific density of structure a_b = data(2); % d, age at birth a_p = data(3); % d, age at puberty a_m = data(4); % d, age at death due to ageing W_b = data(5); % g, wet weight at birth W_j = data(6); % g, wet weight at metamorphosis W_p = data(7); % g, wet weight at puberty W_m = data(8); % g, maximum wet weight R_m = data(9); % #/d, maximum reproduction rate s_M = (W_j/ W_b)^(1/3); % -, acceleration factor l_b = s_M * (W_b/ W_m)^(1/3); % -, scaled length at birth l_j = s_M * (W_j/ W_m)^(1/3); % -, scaled length at metamorphosis l_p = s_M * (W_p/ W_m)^(1/3); % -, scaled length at puberty s = (1/ l_b - 1) * log((s_M - l_j)/ (s_M - l_p)); % -, help var a_j = (a_p * log(s_M) + a_b * s)/ (log(s_M) + s); % d, age at metam E_G = d_V * mu_V/ kap_G/ w_V; % J/cm^3, [E_G] costs for structure r_B = log((s_M - l_j)/(s_M - l_p))/ (a_p - a_j); % 1/d, von Bertalanffy growth rate t_Em = a_b/ 3.7/ l_b; % d, max reserve residence time p_M = 3 * r_B * E_G/ (1 - 3 * r_B * t_Em); % J/d.cm^3, spec somatic maintenance k_M = p_M/ E_G; % 1/d, somatic maintenance rate coefficient k_M = fzero(@fnget_kM, k_M, [], l_b, a_b, r_B); % 1/d, somatic maintenance rate coefficient p_M = E_G * k_M; % J/d.cm^3, spec somatic maintenance g = r_B/ (k_M/ 3 - r_B); % -, energy investment ratio t_Em = 1/ g/ k_M; % d, max reserve residence time kap = 1 - R_m/ s_M^3 * l_b^3 * (1.75/ g + 1)/ (0.95 * k_M); % -, allocation fraction to soma kap = fzero(@fnget_kap, kap, [], w_E, d_V, mu_E, g, E_G, l_b, l_j, l_p, W_m, t_Em, p_M, k_J, kap_R, R_m); E_m = E_G/ kap/ g; % J/cm^3, (max) reserve capacity w = E_m * w_E/ d_V/ mu_E;% -, contribution of reserve to weight L_m = (W_m/ (1 + w))^(1/3)/ s_M; % cm, maximum structural length v = L_m/ t_Em; % cm/d, energy conductance p_Am = v * E_m; % J/d.cm^2, max spec assimilation rate L_b = l_b * L_m; % cm, structural length at birth L_j = l_j * L_m; % cm, structural length at metamorphosis L_p = l_p * L_m; % cm, structural length at puberty x_b = g/ (1 + g); % -, see Tab 2.1 of DEB3 alpha_b = 3 * g * x_b^(1/3)/ l_b; % -, see Tab 2.1 of DEB3 uE0 = (3 * g/ (alpha_b - beta0(0, x_b)))^3; % -, see (2.42) of DEB3 E_0 = uE0 * L_m^3 * E_G/ kap; % J, initial reserve options = odeset('RelTol', 1e-10); [L HE] = ode45(@dget_HE, [1e-8; L_b; L_j; L_p], [0; E_0], options, L_b, L_j, L_m, s_M, E_m, v, p_M, E_G, kap, k_J); E_Hb = HE(2,1); E_Hj = HE(3,1); E_Hp = HE(4,1); % J, maturity levels h_a = 4.27/ a_m^3/ k_M/ g; % 1/d^2 ageing acceleration % pack par par = [p_Am; v; kap; p_M; E_G; E_Hb; E_Hj; E_Hp; h_a]; end %% subfunctions function f = fnget_kM(k_M, l_b, a_b, r_B) k_M = max(k_M, 3 * r_B + 1e-10); g = r_B/ (k_M/ 3 - r_B); % -, energy investment ratio f = a_b - get_tb(g, 1, l_b)/ k_M; % set f to zero end function f = fnget_kap(kap, w_E, d_V, mu_E, g, E_G, l_b, l_j, l_p, W_m, t_Em, p_M, k_J, kap_R, R_m) s_M = l_j/ l_b; % -, acceleration factor k_M = p_M/ E_G; % 1/d, somatatic maintenance rate coefficient E_m = E_G/ kap/ g; % J/cm^3, (max) reserve capacity w = E_m * w_E/ d_V/ mu_E;% -, contribution of reserve to weight L_m = (W_m/ (1 + w))^(1/3)/ s_M; % cm, maximum structural length v = L_m/ t_Em; % cm/d, energy conductance L_b = l_b * L_m; % cm, structural length at birth L_j = l_j * L_m; % cm, structural length at metamorphosis L_p = l_p * L_m; % cm, structural length at puberty x_b = g/ (1 + g); % -, see Tab 2.1 of DEB3 alpha_b = 3 * g * x_b^(1/3)/ l_b; % -, see Tab 2.1 of DEB3 uE0 = (3 * g/ (alpha_b - beta0(0, x_b)))^3; % -, see (2.42) of DEB3 E_0 = uE0 * L_m^3 * E_G/ kap; % J, initial reserve options = odeset('RelTol', 1e-10); [L HE] = ode45(@dget_HE, [1e-8; L_b; L_j; L_p], [0; E_0], options, L_b, L_j, L_m, s_M, E_m, v, p_M, E_G, kap, k_J); E_Hp = HE(4,1); % J, maturity at puberty % set f to zero f = R_m - kap_R * (s_M^3 * L_m^3 * k_M * E_G * (1 - kap)/ kap - k_J * E_Hp)/ E_0; end function dHE = dget_HE(L, EH, L_b, L_j, L_m, s_M, E_m, v, p_M, E_G, kap, k_J) % forward integration of E_H and E from L = 0 to L_p, with acceleration % unpack states E_H = EH(1); % J, maturity E = EH(2); % J, reserve V = L^3; % cm^3, structural volume k_M = p_M/ E_G; % 1/d, som maint rate coeff g = E_G/ kap/ E_m; % -, energy investment ratio if L < L_b % embryo r = (E * v/ L - p_M * V/ kap)/ (E + E_G * V/ kap); % 1/d, specific growth rate p_C = E * (v/ L - r); % J/d, mobilisation rate dE_H = (1 - kap) * p_C - k_J * E_H; % J/d, d/dt E_H, change in maturity dE = - p_C; % J/d, d/dt E, change in reserve dL = L * r/ 3; % cm/d, d/dt L, change in length elseif L < L_j % acceleration r = k_M * (L_m/ L_b - 1)/ (1 + 1/ g); % 1/d, specific growth rate p_C = E * (v/ L_b - r); % J/d, mobilisation rate dE_H = (1 - kap) * p_C - k_J * E_H; % J/d, d/dt E_H, change in maturity dL = L * r/3; % cm/d, d/dt L, change in length dE = E_m * 3 * L^2 * dL; % J/d, d/dt E, change in reserve else % juvenile, post metam r = (E * s_M * v/ L - p_M * V/ kap)/ (E + E_G * V/ kap); % 1/d, specific growth rate p_C = E * (s_M * v/ L - r); % J/d, mobilisation rate dE_H = (1 - kap) * p_C - k_J * E_H;% J/d, d/dt E_H, change in maturity dL = L * r/ 3; % cm/d, d/dt L, change in length dE = E_m * 3 * L^2 * dL; % J/d, d/dt E, change in reserve end dHE = [dE_H; dE]/ dL; % J/cm, change in maturity and reserve end