ecostate (Thorson et al. 2025,
n.d.) is an R package that implements a nonlinear state-space
model representing equilibrium conditions from Ecopath (Polovina 1984), biomass dynamics from Ecosim
(CARL Walters, Christensen, and Pauly
1997), and age/stage-structured dynamics from Ecosim version-2
(Carl Walters et al. 2000). The
parameterization is heavily inspired by RPath (Lucey, Gaichas, and Aydin 2020). However, it
also adds features that are not available in EwE or Rpath including:
Estimating equilibrium and dynamical parameters simultaneously: ecostate uses maximum likelihood estimation to identify the values (and standard errors) for parameters representing both ecosystem equilibrium (e.g., equilibrium biomass, production/consumption per biomass), ecosystem dynamics (e.g., ratio of initial to equilibrium biomass), and measurement parameters (e.g., catchability coefficient);
Process errors: ecostate can estimate errors in both biomass index and catch measurements, as well as unexplained variation in dynamics (termed “process errors”). It then estimates the variance of process errors using hierarchical modelling techniques;
Priors and likelihood penalties: ecostate allows the user to specify a Bayesian prior or likelihood penalty on parameters, rather than fixing them a priori.
These features ensure that anyone running ecostate can rapidly update the model using new data or ecosystem assumptions. Results can then be exactly reproduced, allowing model performance to be explored via simulation modelling.
The model starts by defining equilibrium biomass \(\bar{\beta}_s\) for each taxon \(s \in \{1,2,...,S\}\), and using \(i\) and \(j\) to refer to taxa when they are prey or predator. Every taxon is defined as a heterotroph (consumer), autotroph (producer), or detritus pool, which affects dynamics as explained later. This equilibrium depends upon prey production per biomass \(p_i\), the proportion of prey production that is explained in the model (termed “ecotrophic efficiency”) \(e_i\), diet proportion \(d_{i,j}\) for each predator \(j\) and prey \(i\), and predator consumption per biomass \(w_j\):
\[ \underbrace{\bar{\beta}_i}_{\text{Equilibrium biomass as prey}} \times \underbrace{p_i}_{\text{Prey production per biomass}} \times \underbrace{e_i}_{\text{Ecotrophic efficiency}} = \\ \sum_{j=1}^S{\left( \underbrace{d_{i,j}}_{\text{Proportion of diet for predator }j \text{ by prey }i} \times \underbrace{\bar{\beta}_j}_{\text{Equilibrium biomass for predator }j} \times \underbrace{w_j}_{\text{Predator consumption per biomass}} \right)} \] Autotroph and detritus taxa \(j\) have no consumption so \(d_{i,j}\) for these taxa. We can then solve for the \(S \times S\) matrix of equilibrium consumption \(\bar{\mathbf{C}}\) as:
\[ \bar{\mathbf{C}} = \mathbf{D \odot ( 1 (\bar{\beta} \odot w)}^T ) \]
Given these equilibrium values, we then define a differential equation for biomass dynamics which we integrate over time to calculate biomass at times \(t \in \{ 1, 2, ..., T \}\). In parallel, we also integrate catches \(\eta_s(t)\) over time, while resetting \(\eta_s(t)=0\) at the beginning of each time-interval and recording the integrated catch at the end of each interval:
\[ \begin{gather} \frac{\text{d}}{\text{d}t} \mathbf{\beta}(t) = & \left( \underbrace{\mathbf{g}(t)}_{\text{Growth rate}} - \underbrace{\mathbf{m}(t)}_{\text{Natural mortality rate}} - \underbrace{\mathbf{f}(t)}_{\text{Fishing mortality rate}} + \underbrace{\mathbf{\epsilon}(t)}_{\text{Process error}} \right) & \odot \mathbf{\beta}(t) \\ \frac{\text{d}}{\text{d}t} \mathbf{\eta}(t) = & \mathbf{f}(t) & \odot \mathbf{\beta}(t) \end{gather} \]
where biomass growth rate \(g_s(t)\) and natural mortality rate \(m_s(t)\) are calculated from predicted consumption \(\mathbf{C}(t)\) representing the mass \(c_{i,j}(t)\) of prey \(i\) consumed by predator \(j\), and fishing mortality \(f_s(t)\) is an estimated parameter that is informed by catch data.
Consumption \(c_{i,j}(t)\) varies around the equilibrium consumption:
\[ c_{i,j}(t) = \underbrace{\bar{c}_{i,j}(t)}_{\text{equilibrium consumption rate}} \times \underbrace{\frac{x_{i,j}\frac{\beta_j(t)}{\bar{\beta}_j}}{x_{i,j} - 1 + \frac{\beta_j(t)}{\bar{\beta}_j}}}_{\text{predator functional response}} \times \underbrace{\frac{\beta_i(t)}{\bar{\beta}_i}}_{\text{prey functional response}} \] where \(\mathbf{X}\) is the matrix of predator-prey vulnerability parameters \(x_{i,j}\).
Natural mortality is then calculated as:
\[ m_s(t) = \underbrace{\frac{\sum_{j=1}^S c_{s,j}(t)}{\beta_s(t)}}_{\text{Predation rate}} + \begin{cases} \underbrace{p_s(1-e_s)}_{\text{Residual natural mortality rate}} & \text{ if } s \text{ is autotroph or heterotroph} \\ \underbrace{\nu_s}_{\text{Detritus loss rate}} & {\text{ if } s \text{ is detritus}} \end{cases} \] where residual natural mortality \(p_s(1-e_s)\) accounts for predation by unmodeled taxa, senescence and disease.
The residual mortality rate is defined differently for detritus than other taxa, and this detritus loss rate \(\nu_s\) is defined to ensure that net detritus accumulation matches net consumption plus export at equilibrium:
\[ \underbrace{\bar{\beta}_s \nu_s}_{\text{Equil. detritus loss}} = \underbrace{\sum_{i=1}^S \sum_{j=1}^S \mu_j \bar{c}_{i,j} + \sum_{j=1}^S \bar{\beta}_j p_s (1-e_s)}_{\text{Equil. detritus accumulation}} - \underbrace{\sum_{j=1}^S \bar{c}_{s,j}}_{\text{Equil. detritus consumption}} \] where \(\mu_j\) is the proportion of consumption that is not assimilated for predator \(j\) such that \(\sum_{i=1}^S \sum_{j=1}^S \mu_j \bar{c}_{i,j}\) flows to detritus via unassimilated consumption. Similarly, \(p_j(1-e_j)\) is unexplained mortality which we assume flows to consumption with total \(\sum_{j=1}^S \bar{\beta}_j p_s (1-e_s)\).
Biomass growth rate \(g_s(t)\) is also calculated from consumption: \[ g_s(t) = \begin{cases} \frac{p_s}{w_s} \times \frac{\sum_{i=1}^S c_{i,s}(t)}{\beta_s(t)} & \text{ if } s \text{ is heterotroph} \\ \frac{p_s \bar{\beta}_s}{\beta_s(t)} \times \frac{x_{s,s}\frac{\beta_s(t)}{\bar{\beta}_s}}{x_{s,s} - 1 + \frac{\beta_s(t)}{\bar{\beta}_s}} & \text{ if } x \text{ is autotroph} \\ \frac{\sum_{i=1}^S \sum_{j=1}^S \mu_j c_{i,j}(t) + \sum_{j=1}^S \beta_j(t)p_j(1-e_j)}{\beta_s(t)} & \text{ if } s \text{ is detritus} \end{cases} \] This expression states that:
Heterotrophs: The growth rate for heterotrophs is consumption per biomass times the ratio of production and consumption per biomass;
Autotrophs: The autotroph growth rate is calculated so that production is constant but adjusted by a Type-2 functional response using vulnerability \(x_{s,s}\);
Detritus: The detritus growth rate is total detritus gain per detritus biomass.
ecostate also includes features that allow age- and size-structured dynamics (Thorson et al. n.d.) that are inspired by the parameterization from Rpath (Lucey, Gaichas, and Aydin 2020). This requires defining one or more age-structured populations (“multi-stanza groups”) \(g2\), and each multi-stanza group has a growth coefficient \(K_{g2}\), allometric scaling of consumption to weight \(d_{g2}\), age at maturity \(A_{\text{mat},g2}\) or weight at maturity \(W_{\text{mat},g2}\), and density dependence in larval survival \(X_{\text{spawn},g2}\). Each multi-stanza group then has one or more stanzas \(s2[g2]\), where each stanza has a minimum age \(a_{\text{min},s2}\) and therefore corresponds to ages from that minimum to the minimum for the next stanza (or a high upper limit for the oldest stanza). Finally each stanza corresponds to a single “taxon” \(s[s2[g2]]\), such that the biomass for stanza \(s2[g2]\) of age-structured population \(g2\) is \(\beta_{s[s2[g2]]}\).
The link between biomass dynamics for taxa \(s\) and stanza \(s2[g2]\) is established via the generalized von Bertalanffy growth equation:
\[ \frac{\text{d}}{\text{d}t} \omega = \underbrace{H\omega^d}_{\text{Anabolism}} - \underbrace{K\omega}_{\text{Catabolism}} \] Integrating this expression over time where individuals start at zero mass then results in the generalized von Bertalanffy growth function:
\[ \omega(a) = \omega_{\infty} (1 - e^{K(1-d)a})^{\frac{1}{1-d}} \]
We subsequently track abundance-at-age \(\nu_{a^*}\) for \(n_{\Delta}\) fractional ages \(a^*\) within each integer age \(a\), and calculate equilibrium weight-at-age \(\omega_{a^*}\) for those fractional ages. We start by using a first-order explicit Euler approximation to growth:
\[ \bar{\omega}_{a^*+1} = \bar{\omega}_{a^*} + \frac{H \bar{\omega}_{a^*}^d}{n_{\Delta}} - \frac{K \bar{\omega}_{a^*}}{n_{\Delta}} \] We then assume that anabolism will vary linearly with consumption per metabolic demand. After some algebra, this results in:
\[ \bar{\alpha}_{a^*} \underbrace{\left( \frac{\bar{Q}_{s2}}{\sum_{a^{*'} \in s2} \nu_{a^{*'}} \omega_{a^{a'}}^d)} \right)}_{\text{Equilibrium consumption per demand}} = \underbrace{\bar{\omega}_{a^*+1} - \bar{\omega_{a^*}} + \frac{K}{n_{\Delta}} \bar{\omega}_{a^*}}_{\text{Equilibrium anabolism}} \] where the right-hand-side is known from equilibrium weight-at-age, the second-term of the left-hand-side is calculated from equilibrium biomass-dynamics, and we can use this to solve for \(\bar{\alpha}_{a^*}\), which converts the ratio of consumption and metabolim demand to anabolism for a given fractional age.
After calculating \(\bar{\alpha}_{a^*}\) at equilibrium, we can then use non-equilibrium consumption \(Q_{s2}\) to calculate non-equilibrium growth:
\[ \omega_{a^*+1} = \omega_{a^*} + \bar{\alpha}_{a^*} \left( \frac{Q_{s2}}{\sum_{a^{*'} \in s2} \nu_{a^{*'}} \omega_{a^{*'}}^d} \right) - \frac{K}{n_{\Delta}} \omega_{a^*} \] Similarly, we use the Euler approximation to update abundance at age:
\[ \nu_{a^*+1} = \nu_{a^*} e^{-f_{s^*} - m_{s^*}} \] where \(s^*\) is the taxon associated with fractional age \(a^*\) within a given stanza \(s2[g2]\), and \(f_{s^*}\) and \(m_{s^*}\) match the corresponding rates in the differential equation for biomass dynamics. We additionally include a Beverton-Holt stock recruit model based on mature biomass to calculate \(\nu_{a^*}\) for the initial fractional age \(a^* = 1\).
To couple the biomass-dynamics and age- and size-structured submodels, we simply use an Ordinary Differential Equation (ODE) solver to approximate the integral of biomass dynamics \(\beta\) over one time-step. Using the average consumption for each taxon over this interval, we then use the forwards Euler approximation to calculate \(n_{\Delta}\) updates in abundance-at-age and weight-at-age. After these updates, we then calculate total biomass for each taxon, and plug those values into the biomass \(\beta\). These two will closely match in the absence of age-transitions (where animals transition between stanzas as they age), and using the age-structured update instead ensures that biomass \(\beta\) follows the age-structured dynamics while also closely matching their biomass-dynamic interactions.