Ecological drivers of stability and instability in marine ecosystems

Abstract

A stability analysis of the steady state of marine ecosystems is described. The study was motivated by the approximate invariance of biomass in logarithmic size intervals, which is widely observed in marine ecosystems. This invariance is recovered as the steady state of dynamic models of size spectra, which, unlike traditional species-based models of food webs, explicitly account for the mass gained by an individual organism when it eats a prey item. Little is known about the ecological conditions affecting the stability of the steady state, and a new method is developed to examine this. The results show that stability is enhanced by: (a) decreasing the mean predator-to-prey mass ratio (PPMR), (b) increasing the diet breadth of predators, (c) increasing the strength of intrinsic mortality relative to predation mortality, (d) increasing the biomass conversion efficiency. When perturbed from steady state, size spectra develop a wave-like shape, with an average wavelength especially sensitive to the mean PPMR. These waves move from small to large body size at an average speed which depends on the rate of growth of organisms. In contrast to traditional food web models, stability is enhanced as connectance (diet breadth) increases and as food chain length is increased by reducing the PPMR.

Appendix

Linear evolution equation for perturbations

We define the following dimensionless variables:

$$ x=\ln\left(\frac{w}{w_0}\right), \qquad \bar t=\frac{t}{t_0}, \qquad v(x) = e^{(\alpha+\rho+1) x} \frac{\phi(w)}{\phi(w_0)}. $$

If t ₀ is chosen to be $1/(Aw_0^{\alpha+1}\phi(w_0))$, then the dimensionless search rate constant $\bar A$ is equal to 1. In these new variables, and with the constraint that ξ must equal ρ if intrinsic mortality is present (i.e. if μ > 0), the deterministic jump-growth equation (see Eqs. 1 and 2 in the main text) is

$$\begin{array}{rll} \frac{\partial }{\partial \bar t}v(x) &=& e^{-\rho x} \Bigg(-\bar\mu v(x) + \int s(r) \Big[-e^{(\alpha+\rho) r}v(x)v(x-r) \\ &&{\kern50pt}-e^{-\rho r} v(x)v(x+r) \\ &&{\kern50pt}+e^{(\alpha+\rho)r + (\alpha+2\rho)\psi(r)}v(x-\psi(r)) \\ &&{\kern55pt}\times v(x-r-\psi(r))\Big] {\rm d}r \Bigg), \end{array}$$

(11)

where ψ(r) = ln (1 + Ke ^− r) and $\bar\mu = t_0w_0^{-\rho}\mu$ (note that $\bar\mu$ is the same as the dimensionless mortality parameter $\hat\mu/(\hat Au_0)$ that is varied in Fig. 4 in the main text). Henceforth, all variables and parameters will be dimensionless and the overbars will be omitted. Note that v _s(x) = 1 is a steady state of Eq. 11, provided that

$$ 0 \mu\int s(e^r) \left[ -e^{(\alpha+\rho)r}e^{-\rho r} + e^{(\alpha+\rho)r+(\alpha+2\rho)\psi(r)}\right] {\rm d}r. $$

(12)

This equation determines the intrinsic mortality exponent ρ as a function of α, K, μ and the feeding kernel s.

In order to investigate the behaviour of small perturbations to this steady state, we seek solutions of the form v(x, t) = 1 + ϵ(x, t), where |ϵ(x, t)| ≪ 1. Substituting this into Eq. 11, using Eq. 12 to eliminate terms of order ϵ ⁰ and neglecting terms of order ϵ ² and higher, we obtain a linearised equation for ϵ:

$$\begin{array}{rll} \frac{\partial }{\partial t}\epsilon (x) &=& e^{-\rho x} \Bigg(-\mu\epsilon(x) + \int s(e^r)\Big[ -e^{(\alpha+\rho)r}(\epsilon(x)\notag\\ &&{\kern50pt}+\epsilon(x-r)) - e^{-\rho r}(\epsilon(x)+\epsilon(x+r)) \\ && {\kern50pt}+ e^{(\alpha+\rho)r+(\alpha+2\rho)\psi(r)}(\epsilon(x-\psi(r)) \\ &&{\kern50pt}+\epsilon(x-r-\psi(r))) \Big] {\rm d}r. \Bigg) \end{array}$$

(13)

In order to rewrite the right-hand side of this equation in the form of a linear integral operator acting on ϵ, the integral may be written as a sum of three individual integrals and change of integration variable z = f(r) made in each integral so that the integral contains ϵ(z). Recombining the resulting expression into a single integral then leads to

$$\begin{array}{rll} \frac{\partial }{\partial t}\epsilon (x) &=& e^{-\rho x}\left( b_0\epsilon(x) + \int G(x-r)\epsilon(r) {\rm d}r\right) \\ &&-\; M_{\rm s}(x)\epsilon(x), \end{array}$$

(14)

where

$$\begin{array}{rll} b_0 &=& -\mu - \int s(e^r) \left(e^{(\alpha+\rho)r}+e^{-\rho r}\right){\rm d}r, \\ G(r) &=& -s(e^r)e^{(\alpha+\rho)r} - s(e^{-r})e^{\rho r} \notag\\&&+\; s\left(\frac{K}{e^r-1}\right)\frac{K^{\alpha+\rho}e^{(\alpha+2\rho+1)r}}{(e^r-1)^{\alpha+\rho+1}} \nonumber\\ && +\; s\left(e^r-K\right) \frac{ e^{(\alpha+2\rho+1)r}}{(e^r-K)^{\rho+1}} ,\\ M_{\rm s}(x) &=& \left\{ \begin{array}{ll} \mu_{\rm s} e^{p(x-x_{\rm s})} & x\ge x_{\rm s}, \\ 0 & x<x_{\rm s}, \end{array}\right. . \end{array}$$

Here the perturbation damping term M _s(x), introduced to reduce gradually the amplitude of the perturbations for x > x _s, has also been included. Note that the integral term in the definition of b ₀ corresponds to the constant C in Eq. 7 in the main text.

The McKendrick–von Foerster approximation can be obtained by carrying out a Taylor expansion in Eq. 13 up to terms of order K ²:

$$\begin{array}{rll} \frac{\partial }{\partial t}\epsilon (x) &=& -M_{\rm s}(x)\epsilon(x) + e^{-\rho x} \notag\\ &&\times\Bigg(\! -\mu \epsilon(x) \int s(e^r)\Big[\!-e^{-\rho r}\left(\epsilon(x)\epsilon(xr)\right) \\ &&{\kern14pt}\;\; + K e^{(\alpha+\rho-1) r}\left((\alpha+2\rho)(\epsilon(x)+\epsilon(x-r))\right.\notag\\ &&{\kern14pt}\;\; -\left.\epsilon^\prime(x)-\epsilon'(x-r)\right) + \frac{K^2}{2} e^{(\alpha+\rho-2)r}\nonumber \\ &&{\kern14pt}\;\; \times\Big((\alpha+2\rho)(\alpha+2\rho-1)(\epsilon(x)+\epsilon(x-r)) \\ &&{\kern14pt}\;\; + (1-2\alpha-4\rho)(\epsilon^\prime(x)+\epsilon^\prime(x-r)) \\ &&{\kern14pt}\;\; + \epsilon^{\prime\prime}(x)+\epsilon^{\prime\prime}(x-r)\Big) \Big] {\rm d}r \Bigg), \end{array}$$

(15)

where $\epsilon^\prime(x) = (\partial/\partial x)\epsilon(x)$. As before, suitable changes of integration variables can be used to rewrite this equation in the form

$$\begin{array}{rll} \frac{\partial }{\partial t}\epsilon (x) &=& e^{-\rho x}\Bigg( b_0\epsilon(x) + b_1\epsilon^\prime(x) + b_2\epsilon^{\prime}(x) \nonumber \\ &&{\kern24pt}\quad + \int G_0(x-r)\epsilon(r) {\rm d}r\notag\\ &&{\kern24pt}\quad + \int G_1(x-r)\epsilon^{\prime}(r) {\rm d}r \notag\\ &&{\kern24pt}\quad +\int G_2(x-r)\epsilon^{\prime}(r) {\rm d}r \Bigg). \end{array}$$

As discussed below, numerical studies require the assumption that ϵ(x) is zero for x outside some finite interval, say 0 ≤ x ≤ a. In practice, a is taken to be sufficiently large that the damping term reduces ϵ(a) to a negligibly small value. Under the assumption that ϵ(a) = 0, integration by parts of the above equation gives

$$\begin{array}{rll} \frac{\partial }{\partial t}\epsilon (x) &=& e^{-\rho x}\Big( b_0\epsilon(x) + b_1\epsilon^\prime(x) + b_2\epsilon^{\prime\prime}(x) \notag\\ &&{\kern24pt}\quad -(G_1(x)+G_2^\prime(x))\epsilon(0) -G_2(x)\epsilon^\prime(0) \nonumber \\ &&{\kern24pt}\quad +\int G(x-r)\epsilon(r) {\rm d}r \Big) \\ &&-\, M_{\rm s}(x)\epsilon(x), \end{array}$$

(16)

where

$$\begin{array}{rll} b_0 &=& -\mu + \int s(e^r) \left( -e^{-\rho r} + (\alpha+2\rho)e^{(\alpha+\rho-1)r}\vphantom{\frac{K^2}{2}}\right.\notag\\&& {\kern23mm}\quad\quad\; \times \left.\left( K \frac{K^2}{2}(\alpha\rho)e^{-r}\right)\right) {\rm d}r,\\ b_1 &=& \int s(e^r)e^{(\alpha+\rho-1)r} \left(K \frac{K^2}{2}(1\alpha\rho)e^{-r}\right) {\rm d}r, \\ b_2 &=& \int s(e^r) \frac{K^2}{2} e^{(\alpha+\rho-2)r} {\rm d}r, \\ G(r) & = & \Big( -s(e^{-r})e^{\rho r} + K e^{(\alpha+\rho-1)r} \nonumber \\ &&{\kern8pt} \times\left((\rho+1)s(e^r)-e^r s^\prime(e^r)\right)+ \frac{K^2}{2} e^{(\alpha+\rho-2)r}\notag\\ &&{\kern8pt} \times \left( (\rho^2+3\rho+2) s(e^r) -2(\rho+1) \right.\\&& {\kern8pt} \times\left.e^r s'(e^r) + e^{2r} s''(e^r) \right) \Big). \end{array}$$

This is the second-order (K ²) approximation to the jump-growth equation. If all terms of order K ² are also neglected, then the McKendrick–von Foerster (first-order) approximation to the jump-growth equation is obtained.

Both Eq. 14 and its approximation, Eq. 16, are in the form of a linear evolution equation:

$$ \frac{\partial \epsilon }{\partial t} = L\epsilon, $$

(17)

where L is a linear operator.

Numerical solution of the linear evolution equation

To study Eq. 17 numerically, the variable x is restricted to a finite range 0 ≤ x ≤ a. Biologically, this means that perturbations to the steady state are only allowed within this specified size range and the system is assumed to be at steady state outside this range. Furthermore, the perturbations in the range 0 ≤ x ≤ a are assumed to have no effect on the system outside this size range, which therefore remains at steady state. The variable x is also discretised using a fixed step size h. Thus, ϵ(x) is approximated by a finite vector ${\boldsymbol \epsilon}$ and Eq. 17 becomes a system of linear ordinary differential equations, $\dot{\boldsymbol \epsilon} = M{\boldsymbol \epsilon}$, where M is a matrix defined by:

$$ M_{j,k} = \frac{\partial}{\partial {\boldsymbol \epsilon}_k} (L {\boldsymbol \epsilon})_j . $$

The solution to this system of differential equations is

$$ {\boldsymbol \epsilon}(t) = \sum\limits_{i=1}^N a_i {\boldsymbol v}_i e^{\lambda_i t}, $$

(18)

where λ _i and ${\boldsymbol v}_i$ are the eigenvalues and eigenvectors of M and the a _i are the coefficients of expanding the initial perturbation in terms of eigenvectors: $\boldsymbol \epsilon(0) =\sum a_i {\boldsymbol v}_i$. If Re(λ _i ) < 0 for all i, then small perturbations will always decay to zero in the long term. Conversely, if at least one eigenvalue has a positive real part, perturbations will grow.

The eigenvalues and eigenvectors of M can be computed by standard techniques. In the case of the full jump-growth equation with a realistic predator-to-prey mass ratio, taking a step size that is sufficiently small to resolve accurately the kernel G makes the problem computationally infeasible. However, the McKendrick–von Foerster approximations, with or without the order K ² terms, are more tractable.

Calculating approximate wavelength and wave speed of eigenfunctions

The eigenfunctions of the linear operator L are not, in general, plane waves e ^ikx (unless ρ = 0) but are oscillatory complex functions of x. An approximate wavelength may be obtained from the average angular velocity of the eigenfunction in the complex plane. For a plane wave ϵ(x) = e ^ikx, the angular velocity dθ/dx (where $\theta=\arg(\epsilon(x))$) is constant and equal to k. The wavelength is equal to 2π/k. The wavelength of a general eigenfunction, at a particular value of x, may therefore be approximated by

$$ \textrm{wavelength} = 2\pi\left(\frac{{\rm d}\theta}{{\rm d}x}\right)^{-1}, $$

The average wavelength of the eigenfunction over the range 0 ≤ x ≤ a can be approximated by

$$ \textrm{average wavelength} = \frac{2\pi\Delta x}{\Delta \theta}, $$

where Δθ and Δx = a are the net change in $\arg(\epsilon(x))$ and x, respectively, between x = 0 and x = a.

An initial perturbation that is equal to the eigenfunction ϵ(x) leads, after a period of time t, to a perturbation given by ϵ(x, t) = ϵ(x) e ^λt, where λ is the associated eigenvalue. In general, this is a complex-valued solution (both ϵ(x) and λ are complex-valued). A real-valued solution may be obtained by taking the sum of this solution with the corresponding solution for the complex conjugate eigenfunction and eigenvalue (this real-valued solution is the function that is plotted in Fig. 2.) This gives

$$ \epsilon(x,t) = \left( \epsilon_1(x)\cos(\omega t)+ \epsilon_2(x)\sin(\omega t)\right) e^{qt}, $$

where ϵ(x) = ϵ ₁(x) + iϵ ₂(x) and λ = q + iω. From this, it may be seen that the perturbations returns periodically to its original shape (with amplitude multiplied by e ^qt), with period T = 2π/ω. Therefore, the average wave speed is given by

$$ \textrm{average wave speed} = \frac{\textrm{average wavelength}}{\textrm{period}} = \frac{\omega \Delta x}{\Delta\theta}. $$