Introduction & Philosophy
A key design feature of the mizer package is that a multi-species size-spectrum model can be set up using only a small amount of information. This is possible because mizer uses allometric scaling relations and size-based feeding rules to choose sensible default values for parameters that are not explicitly provided by the user.
For example, to set up a model using the newMultispeciesParams() function, the user only needs to provide the maximum sizes (\(w_{\text{max}}\)) for the species in the community. Almost all other physiological, ecological, and mortality parameters can be automatically filled in using defaults.
This vignette details the philosophy, the relationships, and the mathematical derivations behind the default parameter values calculated by mizer. For a description of the model equations themselves, see the vignette on the model description.
Mathematical Derivations
1. Maximum Intake Rate Coefficient (\(h\))
The maximum intake rate coefficient \(h_i\) for species \(i\) is calculated by get_h_default(). The derivation is based on the idea that an average individual must grow from egg size \(w_{\text{min}}\) to maturity size \(w_{\text{mat}}\) in a given time \(\text{age}_{\text{mat}}\) when feeding at a constant feeding level \(f_0\) (which defaults to \(0.6\)).
By definition, the rate of somatic growth is: \[\frac{dw}{dt} = g(w)\]
Prior to maturation, individuals invest nothing in reproduction (\(\psi(w) = 0\)), so growth is equal to the energy available for growth and reproduction \(E_r(w)\): \[g(w) = \max(0, \alpha f(w) h w^n - \text{metab}(w))\]
Under the default assumptions: * The feeding level is constant: \(f(w) = f_0\). * The activity cost is zero: \(k = 0\), so \(\text{metab}(w) = k_s w^p\). * The metabolic exponent equals the intake exponent: \(p = n\).
Substituting these in gives: \[\frac{dw}{dt} = \alpha f_0 h w^n - k_s w^n = (\alpha f_0 h - k_s) w^n\]
We also define the critical feeding level \(f_c\) as the feeding level at which assimilation exactly balances metabolic costs, i.e., \(g(w) = 0\). This implies: \[\alpha f_c h w^n = k_s w^n \implies k_s = \alpha f_c h\]
Using this expression for \(k_s\) in the growth equation yields: \[\frac{dw}{dt} = (\alpha f_0 h - \alpha f_c h) w^n = \alpha h (f_0 - f_c) w^n\]
We can solve this differential equation by separating variables and integrating from birth (\(t = 0\), \(w = w_{\text{min}}\)) to maturity (\(t = \text{age}_{\text{mat}}\), \(w = w_{\text{mat}}\)): \[\int_{w_{\text{min}}}^{w_{\text{mat}}} w^{-n} \, dw = \alpha h (f_0 - f_c) \int_{0}^{\text{age}_{\text{mat}}} dt\]
Evaluating the integrals: \[\frac{w_{\text{mat}}^{1-n} - w_{\text{min}}^{1-n}}{1-n} = \alpha h (f_0 - f_c) \text{age}_{\text{mat}}\]
Solving for \(h\) gives the formula implemented in get_h_default(): \[h = \frac{w_{\text{mat}}^{1-n} - w_{\text{min}}^{1-n}}{\text{age}_{\text{mat}} (1-n) \alpha (f_0 - f_c)}\]
If the age at maturity is not known, mizer attempts to estimate it from von Bertalanffy parameters using: \[\text{age}_{\text{mat}} = - \frac{\ln\left(1 - \left(\frac{w_{\text{mat}}}{w_{\infty}}\right)^{1/b}\right)}{k_{\text{vb}}} + t_0\] If no growth information is provided at all, mizer defaults to \(h = 30\).
2. Search Volume Coefficient (\(\gamma\))
The search volume coefficient \(\gamma_i\) is calculated by get_gamma_default(). It is selected such that if the prey abundance were described by a power-law resource spectrum: \[N_R(w_p) = \kappa w_p^{-\lambda}\] then the resulting encounter rate would lead to the target feeding level \(f_0\) (default \(0.6\)).
The encounter rate \(E_i(w)\) of predator species \(i\) at size \(w\) is: \[E_i(w) = \gamma_i w^q \int \theta_{iR} N_R(w_p) \phi_i(w, w_p) w_p \, dw_p\]
Substituting the power-law resource spectrum: \[E_i(w) = \gamma_i w^q \theta_{iR} \kappa \int w_p^{1-\lambda} \phi_i(w, w_p) \, dw_p\]
If the predation kernel depends only on the predator/prey mass ratio \(x = w/w_p\), then we can substitute \(w_p = w/x\) and \(dw_p = -w/x^2 \, dx\): \[\int w_p^{1-\lambda} \phi_i(w/w_p) \, dw_p = w^{2-\lambda} \int_0^\infty x^{\lambda - 2} \phi_i(x) \, dx\]
Thus, the encounter rate scales with size as: \[E_i(w) = \gamma_i w^{q + 2 - \lambda} \theta_{iR} \kappa \int_0^\infty x^{\lambda - 2} \phi_i(x) \, dx\]
Let us define the size-independent component of the encounter rate as: \[A_{\text{avail}, i} = \theta_{iR} \kappa \int_0^\infty x^{\lambda - 2} \phi_i(x) \, dx\] so that \(E_i(w) = \gamma_i w^{q + 2 - \lambda} A_{\text{avail}, i}\).
The feeding level is given by: \[f_i(w) = \frac{E_i(w)}{E_i(w) + h_i w^n}\]
For the feeding level to be independent of size and equal to \(f_0\), we must require the size scaling exponents of encounter and maximum intake to match: \[q + 2 - \lambda = n \implies q = n + \lambda - 2\] When this holds, the size \(w\) cancels out: \[f_0 = \frac{\gamma_i A_{\text{avail}, i}}{\gamma_i A_{\text{avail}, i} + h_i}\]
Solving this equation for \(\gamma_i\) yields: \[\gamma_i = \frac{h_i}{A_{\text{avail}, i}} \frac{f_0}{1 - f_0}\]
In the code, \(A_{\text{avail}, i}\) is calculated numerically by setting \(\gamma_i = 1\) and evaluating the encounter rate at the maximum size grid point divided by the appropriate power of \(w\): \[\text{avail_energy}_i = \frac{E_{i}|_{\gamma_i=1}(w_{\text{max}})}{w_{\text{max}}^{q + 2 - \lambda}}\] Then, the default \(\gamma\) is: \[\gamma_i = \frac{h_i}{\text{avail_energy}_i} \frac{f_0}{1 - f_0}\]
The standard metabolic rate coefficient \(k_{s,i}\) is calculated by get_ks_default(). It is chosen so that the critical feeding level \(f_c\) (default \(0.2\)) is exactly the feeding level required for assimilation to balance metabolism at maturity size \(w_{\text{mat}}\).
At the critical feeding level \(f_c\), the assimilation rate equals metabolic losses. At size \(w_{\text{mat}}\): \[\alpha_i f_c h_i w_{\text{mat}}^n = \text{metab}_i(w_{\text{mat}})\]
Using the default metabolic cost equation \(\text{metab}_i(w) = k_{s,i} w^p\) (with activity \(k = 0\)): \[\alpha_i f_c h_i w_{\text{mat}}^n = k_{s,i} w_{\text{mat}}^p\]
Solving for \(k_{s,i}\) gives the formula: \[k_{s,i} = f_c \alpha_i h_i w_{\text{mat}}^{n - p}\]
If \(n = p\), this simplifies to: \[k_{s,i} = f_c \alpha_i h_i\]
4. Target Feeding Level (\(f_0\))
For models where the search volume coefficient \(\gamma_i\) is provided by the user but the target feeding level \(f_0\) is not, get_f0_default() calculates the inverse of the search volume default.
Using the relationship derived above: \[f_{0,i} = \frac{\gamma_i A_{\text{avail}, i}}{\gamma_i A_{\text{avail}, i} + h_i}\] Since the encounter rate calculated with the given \(\gamma_i\) is \(E_i(w) = \gamma_i w^{q+2-\lambda} A_{\text{avail}, i}\), we evaluate: \[\text{avail_energy}_i = \frac{E_i(w_{\text{max}})}{w_{\text{max}}^{q+2-\lambda}}\] And then calculate: \[f_{0,i} = \frac{1}{h_i / \text{avail_energy}_i + 1}\]
