Persistence in Stochastic Lotka--Volterra food chains with intraspecific competition

TL;DR

This paper analyzes the persistence and extinction of species in a stochastic Lotka-Volterra food chain with intraspecific competition, providing criteria based on invasion rates and extending previous work to include predator intraspecific competition.

Contribution

It introduces a comprehensive analysis of stochastic food chains with intraspecific competition, classifying species persistence using invasion rates and extending results to degenerate noise cases.

Findings

Classified species extinction and persistence based on invasion rates.

Provided criteria for food chains of length up to 3.

Extended analysis to both non-degenerate and degenerate noise.

Abstract

This paper is devoted to the analysis of a simple Lotka-Volterra food chain evolving in a stochastic environment. It can be seen as the companion paper of Hening and Nguyen (J. of Math. Biol. `18) where we have characterized the persistence and extinction of such a food chain under the assumption that there is no intraspecific competition among predators. In the current paper we focus on the case when all the species experience intracompetition. The food chain we analyze consists of one prey and $n-1$ predators. The $j$th predator eats the $j-1$st species and is eaten by the $j+1$st predator; this way each species only interacts with at most two other species - the ones that are immediately above or below it in the trophic chain. We show that one can classify, based on the invasion rates of the predators (which we can determine from the interaction coefficients of the system via an algorithm), which species go extinct and which converge to their unique invariant probability measure. We obtain stronger results than in the case with no intraspecific competition because in this setting we can make use of the general results of Hening and Nguyen (Ann. of Appl. Probab.). Unlike most of the results available in the literature, we provide an in depth analysis for both non-degenerate and degenerate noise. We exhibit our general results by analysing trophic cascades in a plant--herbivore--predator system and providing persistence/extinction criteria for food chains of length $n\leq 3$.

Equations

$$\begin{array}{rl}d{x}_1(t)& =x_1(t)(a_{10}-a_{11}{x}_1(t)-a_{12}{x}_2(t))dt\\ d{x}_2(t)& =x_2(t)(-a_{20}+a_{21}{x}_1(t)-a_{22}{x}_2(t)-a_{23}{x}_3(t))dt\\ & \text{makebox}[0.0pt]⋮\\ d{x}_{n-1}(t)& =x_{n-1}(t)(-a_{n-1,0}+a_{n-1,n-2}{x}_{n-2}(t)-a_{n-1,n-1}{x}_{n-1}(t)-a_{n-1,n}{x}_n)dt\\ d{x}_n(t)& =x_n(t)(-a_{n0}+a_{n,n-1}{x}_{n-1}(t)-a_{n,n}{x}_n(t))dt.\end{array}$$

$$\underset{t\to \infty }{\mathrm{limsup}}{x}_i(t)>0,i=1,\dots ,n.$$

$$\underset{t\to \infty }{\lim }{x}_i(t)=0.$$

$$\begin{array}{rl}d{X}_1(t)& =X_1(t)(a_{10}-a_{11}{X}_1(t)-a_{12}{X}_2(t))dt+X_1(t)d{E}_1(t)\\ d{X}_2(t)& =X_2(t)(-a_{20}+a_{21}{X}_1(t)-a_{22}{X}_2(t)-a_{23}{X}_3(t))dt+X_2(t)d{E}_2(t)\\ & \text{makebox}[0.0pt]⋮\\ d{X}_{n-1}(t)& =X_{n-1}(t)(-a_{n-1,0}+a_{n-1,n-2}{X}_{n-2}(t)-a_{n-1,n-1}{X}_{n-1}(t)-a_{n-1,n}{X}_n)dt\\ & +X_{n-1}(t)d{E}_{n-1}(t)\\ d{X}_n(t)& =X_n(t)(-a_{n0}+a_{n,n-1}{X}_{n-1}(t)-a_{nn}{X}_n(t))dt+X_n(t)d{E}_n(t)\end{array}$$

$$\begin{array}{rl}-a_{11}{x}_1-a_{12}{x}_2& =-{\tilde{a}}_{10}\\ a_{21}{x}_1-a_{22}{x}_2-a_{23}{x}_3& ={\tilde{a}}_{20}\\ & \text{makebox}[0.0pt]⋮\\ a_{j-1,j-2}{x}_{j-2}-a_{j-1,j-1}{x}_{j-1}-a_{j-1,j}{x}_j& ={\tilde{a}}_{j-1,0}\\ a_{j,j-1}{x}_{j-1}-a_{jj}{x}_j& ={\tilde{a}}_{j0}.\end{array}$$

$${\mathcal{I}}_{j+1}=-{\tilde{a}}_{j+1,0}+a_{j+1,j}{x}_j^{(j)}.$$

$${\mathbb{R}}_{+}^{(j)}:=\{\mathbf{x}=(x_1,\dots ,x_n)\in {\mathbb{R}}_{+}^n:x_k=0\text{ for }j<k\le n\}\subset {\mathbb{R}}_{+}^n,$$

$${\mathbb{R}}_{+}^{(j),\circ }:=\{\mathbf{x}=(x_1,\dots ,x_n)\in {\mathbb{R}}_{+}^n:x_k>0\text{ for }k\le j;x_k=0\text{ for }j<k\le n\}.$$

$$\Vert \mu ,\nu {\Vert }_{\text{TV}}:=\underset{A\in \mathcal{B}({\mathbb{R}}_{+}^n)}{\sup }|\mu (A)-\nu (A)|.$$

$$P_{\mathbf{X}}(t,\mathbf{x},\cdot )\Rightarrow {\pi }^{\ast },\text{as}t\to \infty ,\mathbf{x}\in {\mathbb{R}}_{+}^{n,\circ }$$

$$\underset{t\to \infty }{\lim }{\mathbb{E}}_{{\mathbf{x}}_0}f(\mathbf{X}(t))={\int }_{{\mathbb{R}}_{+}^n}f({\mathbf{x}}^{\prime }){\pi }_{j^{\ast }}(d{\mathbf{x}}^{\prime }).$$

$$\underset{t\to \infty }{\lim }\Vert P_{\mathbf{X}}(t,\mathbf{x},\cdot )-{\pi }^{\ast }(\cdot ){\Vert }_{\text{TV}}=0,$$

$$P_{\mathbf{X}}(t,\mathbf{x},\cdot )\Rightarrow {\pi }^{\ast },\text{as}t\to \infty$$

$${\mathbb{P}}_{\mathbf{x}}\left\{\underset{t\to \infty }{\lim }{X}_i(t)=0\right\}=1.$$

$$\underset{t\to \infty }{\mathrm{liminf}}{\mathbb{P}}_{\mathbf{x}}\left\{(X_1(s),\dots ,X_k(s))\in K_{\epsilon }\right\}ds\ge 1-\epsilon ,\text{ for any }\mathbf{x}\in {\mathbb{R}}_{+}^{n,\circ },$$

$${\mathbb{P}}_{\mathbf{x}}\left\{\underset{t\to \infty }{\lim }\frac{1}{t}{\int }_0^t{X}_k(s)ds=x_k^{(n)}>0,k=1,\dots ,n\right\}=1$$

$$\underset{t\to \infty }{\lim }\frac{1}{t}{\int }_0^t{X}_i(s)ds=\{\begin{array}{ll}{x}_i^{(j^{\ast })}& \text{ if }i=1,\dots ,j^{\ast },\\ 0& \text{ if }i=j^{\ast }+1,\dots ,n.\end{array}$$

$$\underset{t\to \infty }{\lim }\left(\underset{\mathbf{x}\in K}{\sup }\left|{\int }_{{\mathbb{R}}_{+}^{(j^{\ast }),\circ }}f({\mathbf{x}}^{\prime }){\pi }_{j^{\ast }}(d{\mathbf{x}}^{\prime })-{\mathbb{E}}_{\mathbf{x}}f(\mathbf{X}(t))\right|\right)=0.$$

$$d{X}_i(t)=X_i(t)f_i(\mathbf{X}(t))dt+X_i(t)d{E}_i(t),i=1,\dots ,n$$

$${\tilde{\Pi }}_t(B):=\frac{1}{t}{\int }_0^t{1}_{\{\mathbf{X}(s)\in \cdot \}}ds,t>0,B\in \mathcal{B}({\mathbb{R}}_{+}^n)$$

$$\Big{\{}(x_{1},\dots,x_{n})\in\mathbb{R}^{n}_{+}:x_{i}>0\text{ for }i\in\{\tilde{n}_{1},\dots,\tilde{n}_{k}\};\text{ and }x_{i}=0\text{ if }i\notin\{\tilde{n}_{1},\dots,\tilde{n}_{k}\}\Big{\}}$$

$${\mathbb{R}}_{+}^{\mu }:=\{(x_1,\dots ,x_n)\in {\mathbb{R}}_{+}^n:x_i=0\text{ if }i\in I_{\mu }^c\}$$

$${\mathbb{R}}_{+}^{\mu ,\circ }:=\{(x_1,\dots ,x_n)\in {\mathbb{R}}_{+}^n:x_i=0\text{ if }i\in I_{\mu }^c\text{ and }{x}_i>0\text{ if }{x}_i\in I_{\mu }\},$$

$${\lambda }_j(\mu ):=$$

$$=$$

$$\frac{\ln X_i(t)}{t}=\frac{\ln X_i(0)}{t}+\frac{1}{t}{\int }_0^t\left[f_i(\mathbf{X}(s))-\frac{{\sigma }_{ii}}{2}\right]ds+\frac{1}{t}{\int }_0^{t}d{E}_i(s).$$

$$\frac{1}{t}{\int }_0^t\left[f_i(\mathbf{X}(s))-\frac{{\sigma }_{ii}}{2}\right]ds$$

$${\lambda }_i(\mu )={\int }_{\partial {\mathbb{R}}_{+}^n}\left(f_i(\mathbf{x})-\frac{{\sigma }_{ii}}{2}\right)\mu (d\mathbf{x}).$$

$$\frac{\ln X_i(0)}{t}+\frac{E_i(t)}{t}$$

$${\mathbb{P}}_{\mathbf{x}}\left\{\underset{t\to \infty }{\lim }\left(\frac{\ln X_i(0)}{t}+\frac{E_i(t)}{t}\right)=0\right\}=1.$$

$$\mathbf{c}=(c_1,\dots ,c_n)\in {\mathbb{R}}_{+}^{n,\circ },\text{ where }{c}_1=1,c_i:=\prod _{j=2}^i\frac{a_{k-1,k}}{a_{k,k-1}},i\ge 2.$$

Full text

We exhibit our general results by analysing trophic cascades in a plant–herbivore–predator system and providing persistence/extinction criteria for food chains of length $n\leq 3$.

D. Nguyen was in part supported by the National Science Foundation under grant DMS-1207667.

1. 1 Introduction

2. 2 Mathematical framework

3. 3 Properties of the invasion rates

4. 3.1 Negative invasion rates

5. 3.2 Case study: equal death, competition and predation rates

6. 3.3 Trophic cascades in a plant–herbivore–predator system

7. 4 Discussion

8. A Proofs

1. Introduction

Biological populations usually do not evolve in isolation and as such it is fundamentally important to determine which species persist and which go extinct in a given ecosystem. The fluctuations of the environment make the dynamics of populations inherently stochastic. Consequently, one needs to account for the combined effects of biotic interactions and environmental fluctuations when trying to determine species richness. Sometimes biotic effects can result in species going extinct. However, if one adds the effects of a random environment, extinction might be reversed into coexistence. In other instances, deterministic systems that coexist become extinct once one takes into account environmental fluctuations. A successful method for studying this interplay is to model the populations as discrete or continuous time Markov processes and study the long-term behavior of these processes ([Che00, ERSS13, EHS15, LES03, SLS09, SBA11, BEM07, BS09, BHS08, CM10, CCL*+*09]).

Even though ecological systems are often more complex than a linear food chain, understanding food chain dynamics has a very long theoretical and empirical history and is a topic that is covered extensively in introductory biology ([OB71, WP93, PDJ*+*92, Han92, OPO95, VZSLR99, PCG00]). In nature the way species interact changes at least seasonally and is extremely complicated. Food chains are simplified ‘caricatures’ of the real world but still offer interesting information about various biological features. Most of the time it is more realistic to model a system as a food web, consisting of an interconnection of food chains. However, in certain instances the model can be simplified to a single food chain and one can glean relevant information by analyzing the properties of the linear food chain analytically. For example, if one has a system with three species in which one is a plant, the second one is a herbivore and the third is a predator then there is no reason to have a direct link between the plant and the predator – in this setting one would get a linear food chain (see [MdR12]).

All ecologists are aware that the world is complex. There is a divide between those who are skeptical of theories based on simplified models (see, for example, [Pol91]) and those who think that ‘simple models can be used like a surgeon’s knife, cutting deftly through the cloying of fat of complicating detail to get at the essential sinews of ecological reality’ ([THE*+*10]). Both from practical and analytical perspectives, theoretical models in ecology must significantly simplify the natural complexities. For food webs one of the main simplifications is the community module, that is represented by a food chain. One of the simplest models of food chains is the Lotka–Volterra one. Even though this model is imperfect and does not describe the behavior of any actual ecosystem, it nevertheless captures some key features which carry over to more realistic and analytically intractable models.

In this paper we study models of food chains of arbitrary length. We assume that there is only one species at each trophic level and that each species eats only the one on the adjacent lower trophic level. Furthermore, the ecosystem is supposed to have no immigration or emigration.

Many of the food-chain models studied in the literature are deterministic and of Lotka-Volterra type. Criteria for persistence and extinction have been studied by [GH79, Gar80, FS85] while the global stability of nonnegative equilibrium points was studied by [So79, Har79].

Usually, individuals of the same species have similar requirements for survival. Sometimes their combined demand for a resource is higher than the supply. As such the individuals have to compete for the resource and during this competition some of the individuals naturally become deprived and are therefore less likely to survive. Ecologists usually term this competition among individuals of the same species as intraspecific competition. Arguably all populations experience some form of intraspecific competition. In particular, all vertebrate top predators in terrestrial ecosystems, with the possible exception of some reptiles, have strong intraspecific competition due to direct aggresion or territoriality ([THE*+*10]).

In most situations, competing individuals do not interact with one another directly. One such situation occurs due to exploitation where individuals are affected by the resource that is left after it has been consumed by others. Another form of indirect competition is due to interference. This happens when one individual will prevent another from exploiting the resource within an area of the habitat. Yet another possible type of interaction between individuals of the same species is intraspecific predation. This is the process of both killing and eating an individual of the same species. These observations suggest that one cannot always ignore intraspecific competition, as has been done in previous work.

One example of a deterministic Lotka-Volterra food chain is given by the system:

[TABLE]

The quantities $(x_{1}(t),\dots,x_{n}(t))$ represent the densities of the $n$ species at time $t\geq 0$. In this model $x_{1}$ describes a prey species, which is at the bottom of the food chain. The next $n-1$ species are predators. Species $1$ has a per-capita growth rate $a_{10}>0$ and its members compete for resources according to the intracompetition rate $a_{11}>0$. Predator species $j$ has a death rate $-a_{j0}<0$, preys upon species $j-1$ at rate $a_{j,j-1}>0$, competes with its own members at rate $a_{jj}>0$ and is preyed upon by predator $j+1$ at rate $a_{j,j+1}>0$. The last species, $x_{n}$, is considered to be the apex predator of the food chain.

In the deterministic setting one says that the system (1.1) is persistent if each solution of $\mathbf{x}(t)=(x_{1}(t),\dots,x_{n}(t))$ with $\mathbf{x}(0)\in\mathbb{R}_{+}^{n,\circ}:=((y_{1},\dots,y_{n}):y_{i}>0,i=1,\dots,n)$ satisfies

[TABLE]

We say that species $i$ goes extinct if

[TABLE]

It is natural to analyze the coexistence of species by looking at the average per-capita growth rate of a popualtion when it is rare. Intuitively, if this growth (or invasion) rate is positive, the respective population increases when rare and can invade, while if the growth is negative, the population decreases and goes extinct. If there are only two populations, coexistence is ensured if each population can invade when it is rare and the other population is stationary ([Tur77, CE89, EHS15]).

There is a general theory for coexistence for deterministic models ([Hof81, Hut84, HS89]). It is shown that a sufficient condition for persistence is the existence of a fixed set of weights associated with the interacting populations such that this weighted combination of the populations’s invasion rates is positive for any invariant measure supported by the boundary (i.e. associated to a sub-collection of populations) - see [Hof81].

In order to take into account environmental fluctuations and their effect on the persistence or extinction of species, one approach is to study systems that have random environmental perturbations. One way to do this is by analysing stochastic differential equations that arise by adding noise to ordinary differential equations. For compact state spaces there are results for persistence in [SBA11]. These results have been generalized in [HN18a] where the authors show how, under some natural assumptions, one can characterize the coexistence and extinction of species living on non-compact state spaces. Some of these results hold not only for stochastic differential equations but also for stochastic difference equations (see [SBA11]), piecewise deterministic Markov processes (see [BL16, HS17]), and for general Markov processes (see [Ben16]).

One stochastic version of (1.1) is the process $\mathbf{X}:=(\mathbf{X}(t))_{t\geq 0}=(X_{1}(t),\dots,X_{n}(t))_{t\geq 0}$ defined by the system of stochastic differential equations

[TABLE]

where $\mathbf{E}(t)=(E_{1}(t),\dots,E_{n}(t))^{T}=\Gamma^{\top}\mathbf{B}(t)$ for an $n\times n$ matrix $\Gamma$ such that $\Gamma^{\top}\Gamma=\Sigma=(\sigma_{ij})_{n\times n}$ and $\mathbf{B}(t)=(B_{1}(t),\dots,B_{n}(t))$ is a vector of independent standard Brownian motions living on the probability space $(\Omega,\mathcal{F},\{\mathcal{F}_{t}\}_{t\geq 0},\mathbb{P})$ with a filtration $\{\mathcal{F}_{t}\}_{t\geq 0}$ satisfying the usual conditions. We denote by $\mathbb{P}_{\mathbf{x}}$ (respectively $\mathbb{E}_{\mathbf{x}}$) the probability measure (respectively the expected value) conditioned on $\mathbf{X}(0)=(X_{1}(0),\dots,X_{n}(0))=\mathbf{x}\in\mathbb{R}_{+}^{n}$.

Remark 1.1**.**

There are a few different ways to add stochastic noise to deterministic population dynamics. We assume that the environment mainly affects the growth/death rates of the populations. This way, the growth/death rates in an ODE (ordinary differential equation) model are replaced by their average values plus random noise fluctuation terms. See [Tur77, Bra02, Gar88, HNY18, EHS15, ERSS13, SBA11, HN18a, Gar84] for more details.

Define the stochastic growth rate $\tilde{a}_{10}:=a_{10}-\frac{\sigma_{11}}{2}$ and the stochastic death rates $\tilde{a}_{j0}:=a_{j0}+\frac{\sigma_{jj}}{2},j=1,\dots,n$. For fixed $j\in\{1,\dots,n\}$ write down the system

[TABLE]

It is easy to show that (1.3) has a unique solution, say $(x^{(j)}_{1},\dots,x^{(j)}_{j})$. Define

[TABLE]

We will show that, if (1.3) has a strictly positive solution $(x^{(j)}_{1},\dots,x^{(j)}_{j})$, the invasion rate of predator $X_{j+1}$ in the habitat of $(X_{1},\dots,X_{j})$ is given by (1.4). The invasion rate of predator $X_{j+1}$ is the asymptotic logarithmic growth $\lim_{t\to\infty}\frac{\log X_{j+1}(t)}{t}$ when $X_{j+1}$ is introduced at a low density in $(X_{1},\dots,X_{j})$. We also set $\mathcal{I}_{1}:=\tilde{a}_{10}$ to be the stochastic growth rate of the prey - this can be seen as the invasion rate of the prey into the habitat, when it is introduced at low densities.

Throughout the paper we define $\mathbb{R}^{n}_{+}=[0,\infty)^{n}$ and for $j=1,\dots,n$

[TABLE]

and

[TABLE]

Definition 1.1**.**

One can define a distance on the space of probability measures living on the space $(\mathbb{R}_{+}^{n},\mathcal{B}(\mathbb{R}_{+}^{n}))$ i.e. the Borel measurable subsets of $\mathbb{R}_{+}^{n}$. This is done by defining $\|\cdot,\cdot\|_{\text{TV}}$, the total variation norm, via

[TABLE]

There are different ways one can define the persistence and extinction of species. We review some of these definitions below.

For a system to be strongly stochastically persistent we require that there exists a unique invariant measure $\pi^{*}$ that does not put any mass on the extinction set $S:=\{\mathbf{x}\in\mathbb{R}_{+}^{n}:\Pi_{i=1}^{n}x_{i}=0\}$ and that the distribution of $\mathbf{X}$ converges in some sense to $\pi^{*}$.

Definition 1.2**.**

The process $\mathbf{X}$ is strongly stochastically persistent if it has a unique invariant probability measure $\pi^{*}$ on $\mathbb{R}^{n,\circ}_{+}$ and converges weakly to $\pi^{*}$, that is

[TABLE]

where $P_{\mathbf{X}}(t,\mathbf{x},\cdot)$ is the transition probability of $\mathbf{X}$. This means that for any continuous function $f:\mathbb{R}^{n}_{+}\mapsto\mathbb{R}$ with $\sup_{\mathbf{x}\in\mathbb{R}^{n}_{+}}|f(\mathbf{x})|\leq 1$ and any $\mathbf{x}_{0}\in\mathbb{R}^{n,\circ}_{+}$

[TABLE]

Remark 1.2**.**

We note that if

[TABLE]

then

[TABLE]

so that convergence in total variation implies weak convergence.

Definition 1.3**.**

The species $X_{i}$ goes extinct if for all $\mathbf{x}\in\mathbb{R}^{n,\circ}_{+}$

[TABLE]

Definition 1.4**.**

The species $(X_{1},\dots,X_{j^{*}})$ are ** persistent in probability* if for any $\varepsilon>0$, there exists a compact set $K_{\varepsilon}\subset\mathbb{R}^{(j^{*}),\circ}_{+}$ such that*

[TABLE]

where $\left(x^{(j^{*})}_{1},\dots,x^{(j^{*})}_{j^{*}}\right)\in\mathbb{R}_{+}^{(j^{*}),\circ}$is the unique solution to (1.3) with $j=j^{*}$

We refer the reader to [Sch12] for a discussion of various forms of persistence. With the above concepts in hand we can formulate our main result.

Theorem 1.1**.**

Suppose $n\geq 2$, and $\mathbf{X}(0)=\mathbf{x}\in\mathbb{R}_{+}^{n,\circ}$. We have the following classification.

• (i)

If $\mathcal{I}_{n}>0$ then $(X_{1},\dots,X_{n})$ is persistent in probability. Moreover,

[TABLE]

*where *$\left(x^{(n)}_{1},\dots,x^{(n)}_{n}\right)\in\mathbb{R}_{+}^{(n),\circ}$is the unique solution of (1.3) with $j=n$.

If $\Sigma$ is positive definite, making the noise non-degenerate, then the food chain $\mathbf{X}$ is strongly stochastically persistent and its transition probability converges to its unique invariant probability measure $\pi^{(n)}$ on $\mathbb{R}_{+}^{n,\circ}$ exponentially fast in total variation.

• (ii)

If there exists $0\leq j^{*}<n$ such that $\mathcal{I}_{j^{*}}>0$ and $\mathcal{I}_{j^{*}+1}<0$ then $X_{j^{*}+1},\dots,X_{n}$ go extinct almost surely exponentially fast, as $t\to\infty$, with rates $\mathcal{I}_{j^{*}+1},-\tilde{a}_{j^{*}+2,0},\dots,-\tilde{a}_{n0}$ respectively. Furthermore, $(X_{1},\dots,X_{j^{*}})$ is persistent in probability and with probability $1$

[TABLE]

*where *$\left(x^{(j^{*})}_{1},\dots,x^{(j^{*})}_{j^{*}}\right)\in\mathbb{R}_{+}^{(j^{*}),\circ}$is the unique solution of (1.3) with $j=j^{*}$.

• (iii)

Suppose that $\mathcal{I}_{j^{*}}>0$ and $\mathcal{I}_{j^{*}+1}<0$ for some $j^{*}<n$. Suppose further that there exists a unique invariant probability measure $\pi_{j^{*}}$ on $\mathbb{R}^{(j^{*}),\circ}_{+}$ such that the transition probability measure of $\mathbf{X}$ restricted on $\mathbb{R}^{(j^{*}),\circ}_{+}$ converges weakly uniformly on each compact set to $\pi_{j^{*}}$. By this we mean that for any continuous and bounded function: $f:\mathbb{R}^{(j^{*}),\circ}_{+}\mapsto\mathbb{R}_{+}$ and for any compact set $K\subset\mathbb{R}^{(j^{*}),\circ}_{+}$, we have

[TABLE]

Then for any $\mathbf{x}\in\mathbb{R}^{n,\circ}_{+}$, the transition probability measure $P(t,\mathbf{x},\cdot)$ of $\mathbf{X}$ converges weakly to $\pi_{j^{*}}$ and as a result $(X_{1},\dots,X_{j^{*}})$ is strongly stochastically persistent.

• (iv)

Suppose that $\mathcal{I}_{j^{*}}>0,\mathcal{I}_{j^{*}+1}<0$ for some $j^{*}<n$ and $\Sigma_{j^{*}}$, the principal submatrix of $\Sigma$ obtained by removing the $j^{*}+1$-th,…, $n$-th rows and columns of $\Sigma$, is positive definite. Then for any $\mathbf{x}\in\mathbb{R}^{n,\circ}_{+}$, the transition probability measure $P(t,\mathbf{x},\cdot)$ of $\mathbf{X}$ converges weakly to $\pi_{j^{*}}$ and as a result $(X_{1},\dots,X_{j^{*}})$ is strongly stochastically persistent.

Remark 1.3**.**

If $\Sigma$ is positive definite then all principal submatrices are positive definite, so in particular $\Sigma_{j^{*}}$ from Theorem 1.1 part (iv) is positive definite.

Remark 1.4**.**

We note that by Theorem 1.1 the food chain persists when $\mathcal{I}_{n}>0$ and goes extinct when $I_{j^{*}+1}<1$ for some $j^{*}\leq n-1$. It is key to note that $\mathcal{I}_{j}$ is independent of the coefficients $(a_{lm}),l>j$.

As such, if we add one extra predator at the top of the food chain the quantities $\mathcal{I}_{j}>0,j=2,\dots,n$ remain unchanged and we get one extra invasion rate $\mathcal{I}_{n+1}$. In this setting, when we have $n$ predators, the system persists if $\mathcal{I}_{n+1}>0$ and goes extinct if $I_{j^{*}+1}<1$ for some $j^{*}\leq n$. This means that the introduction of an apex predator makes extinction more likely.

Remark 1.5**.**

The persistence or extinction of species evolving according to system (1.1) when the intraspecies competition for predators is zero (i.e. $a_{ii}=0,i\geq 2$) has been studied by [GH79]. [HN18b] generalized the results from [GH79] to a stochastic setting. The current paper tackles the case when intraspecies competition is nonzero. We get stronger results than in the case without intracompetition because we are able to make use of the general results from [HN18a]. From a technical point of view, strictly positive intracompetition rates make the process return to compact sets exponentially fast. This fact can then be used to prove exponential convergence to an invariant probability measure or extinction.

Most of the results for stochastic food chains only consider chains of length two. We note that our results are new even in the case of food chains of length three.

Theorem 1.1 extends previous results on stochastic Lotka-Volterra systems in two dimensions (see [LB16, HN18a, Rud03]) to an $n$ dimensional setting. We also generalize the work by [Gar84] where the author gives sufficient conditions for persistence of stochastic Lotka-Volterra type food web models in bounded regions of state space. We note that the main results of [Gar84] only say something about persistence until the first exit time of the process from a compact rectangular region $R_{\gamma}\subset\mathbb{R}_{+}^{n,\circ}$. Once the process exits the region, one cannot say whether the species persist or not. Partial results for the existence of invariant probability measures for stochastic Lotka-Volterra systems have been found in [Pol79]. However, these conditions are quite restrictive and impose artificial constraints on the interaction coefficients. In contrast, our results for persistence and extinction are sufficient and (almost) necessary. Moreover, based on which conditions are satisfied, we can say exactly which species persist and which go extinct.

The paper is organized as follows. In Section 2 we present the mathematical framework from [HN18a] and explain how we can apply it in the current context. The proof of Theorem 1.1 is presented in Appendix A. General properties regarding the invasion rates and algorithms for how one can compute these invasion rates appear in Section 3. In Section 3.3 we study a plant–herbivore–predaftor food chain and look at the trophic cascade effect the predator has. Finally, Section 4 is devoted to discussing our results and comparing them to the literature.

2. Mathematical framework

We rewrite (1.2) as

[TABLE]

where $\mathbf{X}(t):=(X_{1}(t),\dots,X_{n}(t))$. This is a stochastic process that takes values in $\mathbb{R}_{+}^{n}:=[0,\infty)^{n}$ and defined on a complete probability space $(\Omega,\mathcal{F},\{\mathcal{F}_{t}\}_{t\geq 0},\mathbb{P})$ with a filtration $\{\mathcal{F}_{t}\}_{t\geq 0}$ satisfying the usual conditions. We mainly focus on the process $\mathbf{X}$ starting at $\mathbf{x}\in\mathbb{R}^{n,\circ}_{+}=(0,\infty)^{n}$. The random normalized occupation measures are defined as

[TABLE]

where $\mathcal{B}(\mathbb{R}_{+}^{n})$ are the Borel measurable subsets of $\mathbb{R}_{+}^{n}$. Note that $\widetilde{\Pi}_{t}(B)$ tells us the fraction of time the process $\mathbf{X}$ spends in the set $B$ during the duration $[0,t]$.

Let $\mathcal{M}$ be the set of ergodic invariant probability measures of $\mathbf{X}$ supported on the boundary $\partial\mathbb{R}^{n}_{+}:=\mathbb{R}_{+}^{n}\setminus\mathbb{R}_{+}^{n,\circ}$. For a subset $\widetilde{\mathcal{M}}\subset\mathcal{M}$, denote by $\operatorname{Conv}(\widetilde{\mathcal{M}})$ the convex hull of $\widetilde{\mathcal{M}}$, that is the set of probability measures $\pi$ of the form $\pi(\cdot)=\sum_{\mu\in\widetilde{\mathcal{M}}}p_{\mu}\mu(\cdot)$ with $p_{\mu}>0,\sum_{\mu\in\widetilde{\mathcal{M}}}p_{\mu}=1$.

Note that each subspace of $\mathbb{R}^{n}_{+}$ of the form

[TABLE]

for some $\tilde{n}_{1},\dots,\tilde{n}_{k}\in\mathbb{N}$ satisfying $0<\tilde{n}_{1}<\dots<\tilde{n}_{k}\leq n$ is an invariant set for the process $\mathbf{X}$. Thus, any ergodic measure $\mu\in\mathcal{M}$ must be supported in such a subspace, that is, there exist $0<n_{1}<\dots<n_{k}\leq n$ (if $k=0$, there are no $n_{1},\dots,n_{k}$) such that $\mu(\mathbb{R}^{\mu,\circ}_{+})=1$ where

[TABLE]

for $I_{\mu}:=\{n_{1},\dots,n_{k}\}$, $I_{\mu}^{c}:=\{1,\dots,n\}\setminus\{n_{1},\dots,n_{k}\}$,

[TABLE]

and $\partial\mathbb{R}_{+}^{\mu}:=\mathbb{R}_{+}^{n}\setminus\mathbb{R}_{+}^{\mu,\circ}$. For the Dirac measure $\boldsymbol{\delta}^{*}$ concentrated at the origin [math], we have $I_{\boldsymbol{\delta}^{*}}=\emptyset$

Remark 2.1**.**

Note that $\operatorname{Conv}(\mathcal{M})$ is exactly the set of invariant probability measures of the process $\mathbf{X}$ supported on the boundary $\partial\mathbb{R}_{+}^{n}$.

For a probability measure $\mu$ on $\mathbb{R}^{n}_{+}$ we define the $i$th Lyapunov exponent (when it exists) via

[TABLE]

Remark 2.2**.**

To determine the Lyapunov exponents of an ergodic invariant probability measure $\mu\in\mathcal{M}$, one can look at the equation for $\ln X_{i}(t)$. An application of Itô’s Lemma yields that

[TABLE]

If $\mathbf{X}$ is close to the support of an ergodic invariant measure $\mu$ for a long time $t\gg 1$, then

[TABLE]

can be approximated by the average with respect to $\mu$

[TABLE]

On the other hand, the term

[TABLE]

is negligible for large $t$ since

[TABLE]

This implies that $\lambda_{i}(\mu),i=1,\dots,n$ are the Lyapunov exponents of $\mu$.

For $\mathbf{x}=(x_{1},\dots,x_{n})\in\mathbb{R}^{n}$, we define the norm $\|\mathbf{x}\|=\max_{i=1}^{n}|x_{i}|$. Let

[TABLE]

One can easily check that there exists $\gamma_{b}>0$ such that

[TABLE]

Then parts (2) and (3) of Assumption 1.1 in [HN18a] are satisfied and one gets the existence and uniqueness of strong solutions to (2.1). Moreover, if $\mathbf{X}(0)=\mathbf{x}\in\mathbb{R}_{+}^{n,\circ}$ then

[TABLE]

In view of [HN18a, Lemma 2.3], for $\mu\in\mathcal{M}$, $\lambda_{i}(\mu)$ is well-defined and

[TABLE]

The intuition behind equation (2.5) is the following: if we are inside the support of an ergodic invariant measure $\mu$ then we are at an ‘equilibrium’ and the process does not tend to grow or decay. If $\mu$ is an invariant probability measure satisfying $\mu(\mathbb{R}_{+}^{(j),\circ})=1$ then we derive from (2.5) that

[TABLE]

That is, the solution of (1.3) is the vector $(\mathbb{E}_{\mu}X_{1},\dots,\mathbb{E}_{\mu}X_{j})$ of the expected values of $(X_{1},\dots,X_{j})$ at stationarity.

The following assumption is shown in [HN18a] to imply strong stochastic persistence

Assumption 2.1**.**

For any $\mu\in\operatorname{Conv}(\mathcal{M})$ one has

[TABLE]

Extinction is ensured by the following two assumptions.

Assumption 2.2**.**

There exists $\mu\in\mathcal{M}$ such that

[TABLE]

If $\\ R_{+}^{n}\neq\{\boldsymbol{0}\}$, suppose further that for any $\nu\in\operatorname{Conv}(\mathcal{M}_{\mu})$ , we have

[TABLE]

where $\mathcal{M}_{\mu}:=\{\nu^{\prime}\in\mathcal{M}:\operatorname{supp}(\nu^{\prime})\subset\partial\mathbb{R}_{+}^{n}\}.$

Define

[TABLE]

and

[TABLE]

Assumption 2.3**.**

Suppose that one of the following is true

• •

$\mathcal{M}^{2}=\emptyset$**

• •

For any $\nu\in\operatorname{Conv}(\mathcal{M}^{2})$, $\max_{\{i=1,\dots,n\}}\left\{\lambda_{i}(\nu)\right\}>0.$

Remark 2.3**.**

We refer the reader to [HN18a] for a detailed discussion of the above assumptions. In short

• •

From a dynamical point of view, the solution in the interior domain $\mathbb{R}^{n,\circ}_{+}$ is persistent if every invariant probability measure on the boundary is a “repeller”. In a deterministic setting, an equilibrium is a repeller if it has a positive Lyapunov exponent (or the eigenvalue of the Jacobian). In a stochastic model, the ergodic invariant measures $\mu\in\mathcal{M}$ play the same role. The $\lambda_{i}(\mu),i=1,\dots,n$ are the Lyapunov exponents of $\mu$ (it can also be seen that $\lambda_{i}(\mu)$ gives the long-term growth rate of $X_{i}(t)$ if $\mathbf{X}$ is close to the support of $\mu$). As a result, if $\max_{i=1}^{n}\{\lambda_{i}(\mu)\}>0$, then the invariant measure $\mu$ is a “repeller”. Therefore, Assumption 2.1 guarantees the persistence of the population.

• •

If an ergodic invariant measure $\mu$ with support on the boundary is an “attractor”, it will attract solutions starting nearby. Intuitively, condition (2.7) forces $X_{i}(t),i\in I_{\mu}^{c}$ to get close to [math] if the solution starts close to $\mathbb{R}^{\mu,\circ}_{+}$.

• •

In order to characterize extinction we need the additional Assumption 2.3 which ensures that apart from those in $\operatorname{Conv}(\mathcal{M}^{1})$, invariant probability measures are “repellers”.

Remark 2.4**.**

*The quantity $\lambda_{i}(\mu)$ can be interpreted as the stochastic growth rate of species $X_{i}$ when introduced at a low density in the habitat consisting of species $\{X_{j},j\in I_{\mu}\}$. Since $\mu$ is a invariant probability measure, the growth rate of any $X_{j},j\in I_{\mu}$ is [math]. *

Example 2.1*.*

Let us start by analyzing the one-dimensional equation for the prey

[TABLE]

In this case $\mathcal{M}=\{\boldsymbol{\delta}^{*}\}$. One can than easily check that

[TABLE]

According to [HN18a, Example 6.2] if $\tilde{a}_{10}>0$ there exists a unique invariant probability measure $\pi^{(1)}$ on $\mathbb{R}_{+}^{\circ}$ and $Z_{1}$ converges exponentially fast to $\pi^{(1)}$. If $\tilde{a}_{10}<0$ then $Z_{1}$ goes extinct. Next, assume one has the prey and one predator

[TABLE]

In view of the analysis from [HN18a, Example 6.2] if $\mathcal{I}_{1}=\lambda_{1}(\boldsymbol{\delta}^{*})=\tilde{a}_{10}<0$ then $X_{1}(t),X_{2}(t)$ converge to 0 almost surely with the exponential rates $\mathcal{I}_{1}=\tilde{a}_{10}$ and $\lambda_{2}(\boldsymbol{\delta}^{*})=-a_{20}-0.5\sigma_{22}$ respectively.

If $\mathcal{I}_{1}>0$ then there exists an invariant measure $\mu_{1}$ on $\mathbb{R}^{\circ}_{1+}:=\{(x_{1},0):x_{1}>0\}$ and

[TABLE]

Now one can compute

[TABLE]

If $\mathcal{I}_{1}>0,\mathcal{I}_{2}<0$ then $Z_{2}$ converges to [math] almost surely with the exponential rate $\lambda_{2}(\mu_{1})$ and the occupation measure of the process $(Z_{1},Z_{2})$ converges to $\mu_{1}$.

If $\mathcal{I}_{1}>0,\mathcal{I}_{2}>0$ the transition probability of $(Z_{1}(t),Z_{2}(t))$ on $\mathbb{R}^{\circ}_{12+}$ converges to an invariant probability measure in total variation with an exponential rate. The case with two predators is treated in [HN18a, Example 6.2].

3. Properties of the invasion rates

We want to say more about the invasion rates $\mathcal{I}_{n+1}$. For this we note by (1.4) that we have to analyze the system (1.3). This can be written in matrix form as

[TABLE]

where $\mathbf{x}^{(n)}=\left(x_{1}^{(n)},\dots,x_{n}^{(n)}\right)^{T}$, $\mathbf{a}=(-\tilde{a}_{10},\tilde{a}_{20},\tilde{a}_{30},\dots,\tilde{a}_{n0})^{T}$ and

[TABLE]

is a tridiagonal $n\times n$ matrix.

It is well-known that the solution can be obtained by a forward sweep that is a special case of Gaussian elimination (see [Mal01]). To simplify notation we let

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Define new coefficients $(c_{1}^{\prime},\dots,c_{j-1}^{\prime}),(d_{1}^{\prime},\dots,d_{j}^{\prime})$ recursively as follows

[TABLE]

and

[TABLE]

Having defined these coefficients the solution to (3.1) can be written as

[TABLE]

Since $\mathcal{I}_{n+1}$ only depends on $x_{n}$, all one needs to do is solve for $d_{n}^{\prime}$. In particular, we can compute directly $\mathcal{I}_{n}$ for $n\leq 4$ as follows:

[TABLE]

Remark 3.1**.**

Since the persistence conditions are $\mathcal{I}_{1},\dots,\mathcal{I}_{4}>0$ we note that they are more likely to be met for the top predator if the growth rate of the prey increases. This agrees with the prediction that the length of a food chain should be an increasing function of the prey growth rate. The fact that the length of a food chain should increase with increasing prey growth rates is a general feature of a multitude of models.

Remark 3.2**.**

The invasion rates are functions of the variances $(\sigma_{ii})_{i=1,\dots,n}$. We note that (at least for $j\leq 4$) $\mathcal{I}_{j}(\sigma_{1},\dots,\sigma_{j})$ is strictly decreasing in each variable $\sigma_{i},1\leq i\leq 4$. Even though we were not able to give explicit formulas for $\mathcal{I}_{n}$ one can see from Proposition 3.1 that for any $1\leq j\leq n$ the quantity

[TABLE]

is strictly decreasing in the variable $\sigma_{u}$ for $1\leq u\leq j$ and independent of $\sigma_{u}$ for $u>j$. As a result environmental stochasticity is seen to increase the risk of extinction.

In the limit of no noise (i.e. $\sigma_{ii}\downarrow 0$ for $1\leq 1\leq 4$) the invasion rates from (3.5) converge to $\widehat{\mathcal{I}}_{i}$, that is $\mathcal{I}_{i}\uparrow\widehat{\mathcal{I}}_{i}$ as $\sigma_{ii}\downarrow 0$, where

[TABLE]

Since we do not assume $\Gamma$ is positive definite we note that our method also works in the deterministic setting. The expressions for $\widehat{\mathcal{I}}_{1},\dots,\widehat{\mathcal{I}}_{4}$ give, correctly, the deterministic invasion rates.

3.1. Negative invasion rates

For fixed $j\in\{1,\dots,n\}$ write down the system

[TABLE]

and suppose it has a strictly positive solution $(x^{(j)}_{1},\dots,x^{(j)}_{j})$. The invasion rate of predator $X_{j+1}$ in the habitat of $(X_{1},\dots,X_{j})$ is given by

[TABLE]

Now, one can look at the system

[TABLE]

and its solution $(x^{(j+1)}_{1},\dots,x^{(j+1)}_{j+1})$.

Proposition 3.1**.**

The following holds

[TABLE]

where $c_{j}^{\prime}$ is defined in (3.2). In particular, $x_{j+1}^{(j+1)}>0$ if and only if $\mathcal{I}_{j+1}>0$.

Proof.

Using (3.1)

[TABLE]

and noting that $d_{j+1}^{\prime}=x_{j+1}^{(j+1)}$, $d_{j}^{\prime}=x_{j}^{(j)}$ we get

[TABLE]

This, together with the expression

[TABLE]

implies that

[TABLE]

Using (3.2) one can easily see that $c^{\prime}_{j}\geq 0$. ∎

Proposition 3.2**.**

If there exists $j^{*}\geq 1$ such that $\mathcal{I}_{j^{*}+1}<0$ then there exists no solution in $\mathbb{R}_{+}^{m,\circ}$ for the system (3.7) with $j=m\in\{j^{*}+1,\dots,n\}$.

Proof.

By Proposition 3.1 we note that

[TABLE]

As a result

[TABLE]

and by Proposition 3.1 $x_{j^{*}+2}^{(j^{*}+2)}<0$. By repeating this argument we see that $x_{m}^{(m)}<0$

∎

3.2. Case study: equal death, competition and predation rates

Consider a simplified setting where $a_{ii}=\alpha,i=1,\dots,n$, $a_{i,i-1}=\beta,i=2,\dots,n$, $a_{i,i+1}=\beta,i=1,\dots,n-1$ and $\tilde{a}_{10}=\delta,\tilde{a}_{i0}=\gamma,i=2,\dots,n$. In this case we want to solve

[TABLE]

where $\mathbf{x}=\left(x_{1},\dots,x_{n}\right)^{T}$, $\mathbf{a}=(-\delta,\gamma,\gamma,\dots,\gamma)^{T}$ and

[TABLE]

We use the technique of [Dow08][Section 3.1] to find the inverse $A^{-1}$ of $A$. The quadratic equation

[TABLE]

hast the distinct roots

[TABLE]

One can write the $n$th row of the inverse matrix $A^{-1}$ as

[TABLE]

Therefore the solution to (3.10) satisfies

[TABLE]

In this case one can write down explicitly the formula for the invasion rate.

[TABLE]

3.3. Trophic cascades in a plant–herbivore–predator system

Let us explore a food chain with two or three species. We will assume $X_{1}$ is a plant, $X_{2}$ is a herbivore eating the plant, and $X_{3}$ is a predator that preys on the herbivore.

We can compute the expected abundances of different species at stationarity using the linear system (3.7). We do this to glean information regarding how these expected abundances are changed by intraspecific competition and environmental stochasticity. According to our notation from Section 3.1 the quantity $x_{i}^{(j)}$ will denote the abundance of species $i$ at stationarity when there are $j$ species present. In the example with three species we are looking at, $x_{1}^{(3)}$ will be the abundance at stationarity of the plant when we have the plant, the herbivore, and the predator present. In contrast, $x_{1}^{(2)}$ is the abundance of the plant at stationarity when we only have the plant and the herbivore present.

Solving the $3\times 3$ or $2\times 2$ system (3.7) directly yields:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Let us explore the effect the introduction of a predator $X_{3}$ has on the expected density of the plant and the herbivore at stationarity. Using the formulas above one can show that, as long as the predator $X_{3}$ persists, i.e. $x_{3}^{(3)}>0$, we will always have $x_{1}^{(3)}-x_{1}^{(2)}>0$. If the predator $X_{3}$ goes extinct, i.e. $x_{3}^{(3)}=0$, then $x_{1}^{(3)}-x_{1}^{(2)}=0$. Similarly, one can show that if $x_{3}^{(3)}>0$ then $x_{2}^{(3)}-x_{2}^{(2)}<0$ and if $x_{3}^{(3)}=0$ then $x_{2}^{(3)}-x_{2}^{(2)}=0$. One can also note that the abundance of the plant species is decreasing as we increase the death rate of the predator.

In order to get more information, we graph $x_{1}^{(3)}-x_{1}^{(2)}$ and $x_{2}^{(3)}-x_{2}^{(2)}$ as functions of the predation rate of the predator on the herbivore,$a_{32}$, and the intracompetition rate of the predator, $a_{33}$. See Figures 3 and 5. Similarly, in Figure 4 (respectively Figure 6) we graph $x_{1}^{(3)}-x_{1}^{(2)}$ (respectively $x_{1}^{(3)}-x_{1}^{(2)}$) as a function of predation rate $a_{32}$, the death rate of the predator $\tilde{a}_{33}$ or the intracompetition rate of the predator $a_{33}$. For Figures 1, 3, 4, 5 and 6 we have set $\tilde{a}_{10}=4$ and all the other coefficients (other than the ones being varied) equal to $1$. For Figure 2 we have set all constant coefficients equal to $1$.

We note that the introduction of the predator is always beneficial to the plant and detrimental to the herbivore. The predator will decrease the population size of the herbivore, which will lead to an increase in the plant population size. The density of the plant is seen to increase as we increase the predation rate $a_{32}$ of the predator on the herbivore, and as we decrease the intraspecific competition rate, $a_{33}$, among predators. Plant density will also increase if the stochastic death rate $\tilde{a}_{30}$ of the predator decreases. The areas of the graphs from Figures 3 and 4 where $x_{1}^{(3)}-x_{1}^{(2)}=0$ (or from Figures 5 and 6 where $x_{2}^{(3)}-x_{2}^{(2)}=0$) are those where the predator $X_{3}$ goes extinct. It turns out that anything that is detrimental to the predator (higher intracompetition rate or higher death rate), is also detrimental to the plant. Similarly, factors that are helping the predator survive (higher predation rate $a_{32}$) increase the density of the plant. If one looks at the herbivore $X_{2}$ then its abundance at stationarity will always suffer by the introduction of the predator.

4. Discussion

Even though environmental stochasticity is often said to be a key factor in the study of the persistence of species, its effect on persistence has not been investigated until recently. [BHS08] showed that if one adds a small diffusion term to a persistent deterministic system then the corresponding differential equation has a positive stationary distribution concentrated on the positive global attractor of the deterministic system. For many systems the random perturbations might not be small. For populations living in a compact state space [SBA11] give sufficient conditions for persistence that extend the results from deterministic systems to randomly forced nonlinear systems. This has been further extended by [HN18a] to non-compact states spaces. [HN18a] are able to give, under some mild assumptions, sufficient and necessary conditions for persistence and extinction for stochastic Kolmogorov systems of the form

[TABLE]

Most results which give sharp, tractable conditions for persistence and extinction of populations in stochastic environments usually treat models with only two species (see [EHS15, Rud03]). However, it has been shown that it makes more sense to look at food chains having more than two species ([HP91, KH94, Pai88]). We note that although our model formulation ignores many important ecological features, it still yields some interesting conclusions. These mathematical conclusions should be further studied as hypthesis requiring confirmation. Our results are a first step towards an analysis of stochastic food chains and food webs with an arbitrary number of species.

In the current paper we use the newly developed methods from [HN18a] to analyse the persistence and extinction of species that are part of a stochastic Lotka-Volterra food chain. We assume that species can only interact with those other species which are adjacent to them in the food chain and that there is strictly positive intraspecies competition for all the species. Our main interest was to lift the results from the deterministic setting to the stochastic one and to see whether stochasticity inhibits or enhances coexistence. By studying the invasion rates of the predators $(\mathcal{I}_{2},\dots,\mathcal{I}_{n})$ we show that one can determine which species persist and which go extinct exponentially fast. Furthermore, we provide in Section 3 an algorithm for computing the invasion rates. In this way, based on the interaction coefficients of the system, one can find sufficient and (almost) necessary conditions for persistence/extinction. We show that the introduction of a new top predator into the ecosystem makes extinction more likely. This agrees with the deterministic case studied in [GH79]. Furthermore, we also note that in our setting stochasticity makes extinction more likely. However, since the invasion rates depend continuously on the covariance matrix $\Sigma$ of the environmental noise one can see that if the random perturbations are small and the associated deterministic system is persistent then the stochastic system is also persistent. Actually, our results (see Remark 3.2) show that stochasticity acts in a bottom-up way: the variance $\sigma_{11}$ affects species $\{1,\dots,n\}$ and the variance $\sigma_{jj},j=2,\dots,n$ affects the species $\{j,j+1,\dots,n\}$. As such, in our model the environmental stochasticity of a trophic level only affects the persistence and extinction of species at higher trophic levels. We note that environmental stochasticity does not always make extinction more likely. For example in [BL16] the authors show that in certain cases, the extinction of species in a deterministic setting can be reversed into coexistence by adding randomness to the system. As such,we think that the rigorous study of the stochastic system we propose, did provide insightful information.

It is noteworthy that we do not only get robust results for extinction or persistence – we also get that the convergence to the stationary distribution in the case of persistence is exponentially fast and an exact expression for the convergence rate to [math] in the case of extinction. These rates are very helpful when one wants to run numerical methods and simulations.

We have fully analysed what happens in chains of length $n\leq 4$. Intraspecific competition is shown to change both the conditions for persistence and the strength of trophic cascades.

Humans have always tried to exterminate predators: slayers of predators were seen as heroes in most mythologies; culls have been used to control seals and sea lions in order to manage fisheries; predator control agents are often hired to kill predators (wolves, coyotes, etc). We have effectively decimated and in some cases even driven to extinction entire species of predators. The effects of these exterminations are now becoming more and more clear. There are a plethora of reasons why predators are important in food webs. Predators are usually at the top of the food chain and thus can regulate the trophic levels below them. Removing predators often destabilizes the food chain, and sets off reactions that can cascade down to the lowest trophic level. In Section 3.3 we looked in depth at a plant–herbivore–predator food chain. What our computations and figures show is the following: Aything that helps the predator (decreased death rate, higher predation rate) will be detrimental to the herbivore and favorable to the plant.

Our model also leads to the observation that food chain length should increase when we increase the stochastic growth rate $\tilde{a}_{10}$ of the plant at the first trophic level.

The only cases that cannot be treated are those for which one of the invasion rates is zero, that is $\mathcal{I}_{k}=0$ for some $k\in\{2,\dots,n\}$. This is where our methods break down. As mentioned in [GH79] even in the deterministic setting, when $\kappa(n)=0$ (which would imply one of the invasion rates is [math]) the problem becomes more complicated: one can find solutions with positive initial conditions which persist when $n=3$ while when $n=4$ there are solutions which are not persistent. In the stochastic case, when $n=1$ the prey is described by the SDE

[TABLE]

If $\mathcal{I}_{1}=\tilde{a}_{10}=0$ then one can show that $X_{1}$ is null recurrent and

[TABLE]

As a result the prey $X_{1}$ is not strongly stochastically persistent (there is no invariant probability measure on $(0,\infty)$) but it also does not go extinct. It only goes extinct in the weak sense given by (4.1). We expect similar phenomena to occur in higher dimensions if one of the invasion rates is zero. One possible approach would be to try and adapt the methods used by [Bax91] where the author studies SDE where the extinction set is $\{0\}$. [Bax91] is able to show that if the leading Lyapunov exponent is zero then the process is null-recurrent. In the setting of [Bax91] one only has to study the dirac measure at [math], something which simplifies the problem significantly.

Our results generalize the results from the deterministic setting of [GH79] to their natural stochastic analogues. We are able to find an algebraically tractable criterion (just like in the deterministic setting) for persistence and extinction.

The invasion rates are shown to be closely related to the first moments of the invariant measures living on the boundary $\partial\mathbb{R}_{+}^{n}$ of the system. This result is the analogue of looking for the different equilibrium points of the deterministic system (1.1) and then studying the stability of these points.

The main simplification of our model is the fact that the dynamics of each trophic level is governed by the adjoining trophic levels which immediately precede or succeed it. This fact makes it possible to explicitly describe the structure of the ergodic invariant probability measures of the system living on the boundary $\partial\mathbb{R}_{+}^{n}$ (Lemma A.1). The key property of an invariant probability measure $\mu$ living on $\partial\mathbb{R}_{+}^{n}$ is that if predator $X_{j}$ is not present then all predators that are above $j$ (that is, $X_{i}$ with $i>j$) are also not present. This fact is biologically clear because if species $X_{j}$ does not exist then $X_{j+1}$ must go extinct since it does not have a food source.

For more complex interactions between predators and their prey (i.e. a food web instead of a food chain), even when $n=3$, the possible outcomes become much more complicated. We refer the reader to [HN18a] for a detailed discussion of the case when one has one prey and two predators and the apex predator eats both the intermediate predator and the prey.

In ecology there has been an increased interest in the spatial synchrony that appears in population dynamics. This refers to the changes in the time-dependent characteristics (i.e. abundances etc) of structured populations. One of the mechanisms which creates synchrony is the dependence of the population dynamics on a synchronous random environmental factor such as temperature or rainfall. The synchronizing effect of environmental stochasticity, or the so-called Moran effect, has been observed in multiple population models. Usually this effect is the result of random but correlated weather effects acting on populations. For many biotic and abiotic factors, like population density, temperature or growth rate, values at close locations are usually similar. We refer the reader interested in an in-depth analysis of spatial synchrony to [KBB*+*00, LKB04]. Most stochastic differential equations models appearing in the population dynamics literature treat only the case when the noise is non-degenerate (although see [Rud03, DNDY16]). Although this approach significantly simplifies the technical proofs, from a biological point of view it is not clear that the noise should not be degenerate. For example, if one models a system with multiple populations then all populations can be influenced by the same factors (a disease, changes in temperature and sunlight etc). Environmental factors can intrinsically create spatial correlations and as such it makes sense to study how these degenerate systems compare to the non-degenerate ones. In our setting the noise of the different species could be strongly correlated. Actually, in some cases it could be more realistic to have the same one-dimensional Brownian motion $(B_{t})_{t\geq 0}$ driving the dynamics of all the interacting species. Therefore, we chose to present a full analysis of the degenerate setting.

Acknowledgments. We thank an anonymous referee for comments which helped improve this manuscript and Sebastian Schreiber for helpful discussions and suggestions.

Appendix A Proofs

The following result tells us that there is no ergodic invariant probability measure $\mu$ that has a gap in the chain of predators.

Lemma A.1**.**

Suppose $\mu\in\mathcal{M}$ such that $I_{\mu}=\{n_{1},\dots,n_{k}\}$. Then $I_{\mu}$ must be of the form $\{1,2,\dots,l\}$ for some $l\geq 1$.

Proof.

We argue by contradiction. First, suppose that $n_{1}>1$. By (2.5)

[TABLE]

which is a contradiction.

Alternatively, suppose that there exists $\mu\in\mathcal{M}$ such that $I_{\mu}=\{1,\dots,u^{*},v^{*},\dots,n_{k}\}$ with $1\leq u^{*}<v^{*}-1\leq n_{k}\leq n$. As a result one can see that $v^{*}-1\notin I^{\mu}$. Then by (2.5)

[TABLE]

which is a contradiction. ∎

For $i=1,\dots,n$, denote by $\mathcal{M}_{i}$ the set of all invariant probability measures $\mu$ of $\mathbf{X}$ satisfying $\mu\left(\mathbb{R}^{(i),\circ}_{+}\right)=1$. For $i=0$, define $\mathcal{M}_{0}=\{\boldsymbol{\delta}^{*}\}$. By Lemma A.1, we have $\operatorname{Conv}(\mathcal{M})=\operatorname{Conv}(\cup_{i=0}^{n-1}\mathcal{M}_{i})$ and $\operatorname{Conv}(\cup_{i=0}^{n}\mathcal{M}_{i})$ is the set of all invariant probability measures of $\mathbf{X}$ on $\mathbb{R}^{n}_{+}$.

Lemma A.2**.**

We have the following claims.

• •

If $\mathcal{I}_{k}\leq 0$ then $\mathcal{I}_{k+1}<0$.

• •

If $\mathcal{I}_{n}\leq 0$, there $\mathbf{X}$ has no invariant probability measure on $\mathbb{R}^{n,\circ}_{+}$.

Proof.

If $\mathcal{I}_{k+1}=-\tilde{a}_{k+1,0}+a_{k+1,j}x^{(k)}_{k}\geq 0$, then $x^{(k)}_{k}>0$. We will show in Section 4 that $x^{(k)}_{k}$ has the same sign as $\mathcal{I}_{k}$. Thus, if $\mathcal{I}_{k+1}\geq 0$ then $\mathcal{I}_{k}>0$, which proves the first claim.

If $\mathbf{X}$ has an invariant probability measure $\mu$ on $\mathbb{R}^{n,\circ}_{+}$, then we must have $\int_{\mathbb{R}^{n}_{+}}x_{n}\mu(d\mathbf{x})=x^{(n)}_{n}$. As a result $x^{(n)}_{n}>0$, which leads to $\mathcal{I}_{n}>0$ since they have the same sign. The second claim is therefore proved. ∎

Lemma A.3**.**

We have the following claims.

1. (1)

For any initial condition $\mathbf{X}(0)=\mathbf{x}\in\mathbb{R}^{n}_{+}$, the family $\left\{\widetilde{\Pi}_{t}(\cdot),t\geq 1\right\}$ is tight in $\mathbb{R}^{n}_{+}$, and its weak∗-limit set, denoted by $\mathcal{U}=\mathcal{U}(\omega)$ is a family of invariant probability measures of $\mathbf{X}$ with probability 1. 2. (2)

Suppose that there is a sequence $(T_{k})_{k\in\mathbb{N}}$ such that $\lim_{k\to\infty}T_{k}=\infty$ and $\left(\widetilde{\Pi}_{T_{k}}(\cdot)\right)_{k\in\mathbb{N}}$ converges weakly to an invariant probability measure $\pi$ of $\mathbf{X}$ when $k\to\infty$ . Then for this sample path, we have $\int_{\mathbb{R}^{n}_{+}}h(\mathbf{x})\widetilde{\Pi}_{T_{k}}(d\mathbf{x})\to\int_{\mathbb{R}^{n}_{+}}h(\mathbf{x})\pi(d\mathbf{x})$ for any continuous function $h:\mathbb{R}^{n}_{+}\to\mathbb{R}$ satisfying $|h(\mathbf{x})|<K_{h}(1+\|\mathbf{x}\|)\,,\,\mathbf{x}\in\mathbb{R}^{n}_{+}$, with $K_{h}$ a positive constant and $\delta\in[0,\delta_{1})$. 3. (3)

For any $\mathbf{x}\in\mathbb{R}^{n,\circ}_{+}$

[TABLE]

and

[TABLE]

Proof.

Let $\tilde{c}_{1}=1,\tilde{c}_{i}:=\prod_{j=2}^{i}\dfrac{a_{k-1,k}}{2a_{k,k-1}}=c_{i-1}\dfrac{a_{i-1,i-1}}{2a_{i,i-1}},i\geq 2.$ Put

[TABLE]

we can easily verify that

[TABLE]

Thus, when $\|x\|$ is sufficiently large, $|\sum_{i=1}^{n}\tilde{c}_{i}f_{i}(\mathbf{x})|\geq\tilde{\gamma}\sum_{i=1}^{n}x_{i}$, which implies

[TABLE]

As a result,

[TABLE]

In other words, Assumption 1.4 of [HN18a] is satisfied by our model. Thus, the first and second claims of this lemma follow from [HN18a, Lemma 4.6, Lemma 4.7]. By Itô’s formula and the definition of $\tilde{\Pi}_{t}$, we have

[TABLE]

By the strong law of large numbers for martingales,

[TABLE]

which leads to (A.1).

(A.2) can be derived by using equation (4.22) of [HN18a] or by mimicking the proof of [DS06, Theorem 2.4]. ∎

Proof of Theorem 1.1 (i).

Since $\mathcal{I}_{n}>0$, it follows from Lemma A.2 that $\mathcal{I}_{k}>0$ for any $k=1,\dots,n$. By Lemma A.1, for any $\mu\in\operatorname{Conv}(\mathcal{M})=\operatorname{Conv}(\cup_{i=0}^{n-1}\mathcal{M}_{i})$, we can decompose $\mu=\rho_{1}\mu_{i_{1}}+\dots+\rho_{k}\mu_{i_{k}}$ where $0\leq i_{1}<\dots<i_{k}\leq n-1$ and $\mu_{i_{j}}\in\mathcal{M}_{i_{j}}$, $\rho_{j}>0$ for $j=1,\dots,k$ and $\sum\rho_{j}=1$. Since $i_{1}<i_{j}$ for $j=2,\dots,k$, we deduce from (2.5) that $\lambda_{i_{1}+1}(\mu_{i_{j}})=0$ for $j=2,\dots,k$. On the other hand, (1.4) and (2.6) imply

[TABLE]

As a result,

[TABLE]

Thus,

[TABLE]

In other words, Assumption 2.1 is satisfied. By Theorem 3.1 of [HN18a], there exist positive $p_{1},\dots,p_{n},T$ and constants $\theta,\kappa\in(0,1)$ such that

[TABLE]

where

[TABLE]

Equation (A.4) and the Markov property of $\mathbf{X}$ lead to

[TABLE]

Thus,

[TABLE]

By [HN18a, Lemma 2.1], there exists $\widehat{K}>0$ such that

[TABLE]

which together with the Markov property implies

[TABLE]

In view of (A.5) and (A.6), we have

[TABLE]

For any fixed $\varepsilon>0$, define $K:=\left\{\mathbf{x}\in\mathbb{R}^{n,\circ}_{+}:V^{\theta}(\mathbf{x})\leq\dfrac{1}{\varepsilon}\exp(\widehat{K}T)\dfrac{K}{1-\kappa}\right\}$ then $K$ is a compact subset of $\mathbb{R}^{n,\circ}_{+}$. The definition of $K$ together with the last inequality yield

[TABLE]

The stochastic persistence in probability is therefore proved.

To prove (1.5), we need to show that for any initial value $\mathbf{x}\in\mathbb{R}^{n,\circ}_{+}$, the weak-limits points of $\tilde{\Pi}_{t}$ are a subset of $\mathcal{M}_{n}$ with probability 1.

Suppose the claim is false. Then, by part (i) of Lemma A.3, we can find $\mathbf{x}\in\mathbb{R}^{n,\circ}_{+}$ and $\tilde{\Omega}_{\mathbf{x}}\subset\Omega$ with $\mathbb{P}_{\mathbf{x}}(\tilde{\Omega}_{\mathbf{x}})>0$ and such that for $\omega\in\tilde{\Omega}_{\mathbf{x}}$, there exists $t_{k}=t_{k}(\omega)$ satisfying that $\lim_{k\to\infty}t_{k}=\infty$ and $\tilde{\Pi}_{t_{k}}(\omega)$ converges weakly to $\mu(\omega)=\rho_{1}\mu_{1}+\rho_{2}\mu_{2}$ where $\mu_{1}\in\operatorname{Conv}(\mathcal{M})$ and $\mu_{2}\in\mathcal{M}_{n}$ and $\rho_{1}>0$. By Lemma A.2, $\lambda_{n}(\mu_{1})>0$. In view of (2.5), $\lambda_{n}(\mu_{2})=0$. Thus, for almost all $\omega\in\tilde{\Omega}_{\mathbf{x}}$, we have from part (ii) of Lemma A.3 that

[TABLE]

which contradicts (A.2). Thus, with probability 1, the weak-limit points of $\tilde{\Pi}_{t}$ as $t\to\infty$ must be contained in $\mathcal{M}_{n}$. Then, (1.5) follows from (2.6).

When $\Sigma$ is positive definite, it follows from [HN18a, Theorem 3.1] that the food chain $\mathbf{X}$ is strongly stochastically persistent and its transition probability converges to its unique invariant probability measure $\pi^{(n)}$ on $\mathbb{R}_{+}^{n,\circ}$ exponentially fast in total variation.

∎

Proof of Theorem 1.1 (ii).

We suppose there exists $j^{*}<n$ such that $\mathcal{I}_{j^{*}}>0$ and $\mathcal{I}_{j^{*}+1}<0$. By Lemma A.2 part (ii), there are no invariant probability measures on $\mathbb{R}^{(j),\circ}_{+}$ for $j=j^{*}+1,\dots,n$. Using Lemma A.1, we see that the set of invariant probability measures on $\mathbb{R}^{n}_{+}$ of $\mathbf{X}$ is $\operatorname{Conv}(\cup_{i=0}^{j^{*}}\mathcal{M}_{i})$.

Note that $\lambda_{j^{*}+1}(\mu)=-\tilde{a}_{j^{*}+1}<0$ if $\mu\in\mathcal{M}_{i}$ for $i<j^{*}$ and $\lambda_{j^{*}+1}(\mu)=\mathcal{I}_{j^{*}+1}<0$ if $\mu\in\mathcal{M}_{j^{*}}$. As a result, $\lambda_{j^{*}+1}(\mu)<0$ for any $\mu\in\operatorname{Conv}(\cup_{i=0}^{j^{*}}\mathcal{M}_{i})$. Similarly, $\lambda_{j}(\mu)<0$ for any $j>j^{*}+1$ and $\mu\in\operatorname{Conv}(\cup_{i=0}^{j^{*}}\mathcal{M}_{i})$. By (A.1) we have that

[TABLE]

Since

[TABLE]

we have

[TABLE]

Using (A.2) and a contradiction argument similar to that in the proof of part (i), we can show that with probability 1, the weak-limit points of $\tilde{\Pi}_{t}$ as $t\to\infty$ must be contained in $\mathcal{M}_{j^{*}}$. Thus, for $\mathbf{x}\in\mathbb{R}^{n,\circ}_{+}$, we have from (A.8), (A.9), and Lemma A.2 that

[TABLE]

and

[TABLE]

To prove the persistence in probability of $(X_{1},\dots,X_{j^{*}})$, we define

[TABLE]

We have proved that $\operatorname{Conv}\left(\bigcup_{j=0}^{j^{*}}\mathcal{M}_{j}\right)$ is the set of invariant probability measures of $\mathbf{X}$ on $\mathbb{R}^{n}_{+}$. Note that $\operatorname{Conv}\left(\bigcup_{j=0}^{j^{*}-1}\mathcal{M}_{j}\right)$ is the set of invariant probability measures of $\mathbf{X}$ on $\partial\mathbb{R}^{(j^{*}),\diamond}$. Since $\mathcal{I}_{j^{*}}>0$, applying (A.3) with $n$ replaced by $j^{*}$ we obtain

[TABLE]

Using this condition, we can imitate the proofs in [HN18a, Section 3] to construct a Lyapunov function $U(\mathbf{x}):\mathbb{R}^{(j^{*}),\diamond}_{+}\mapsto\mathbb{R}_{+}$ of the form

[TABLE]

satisfying

[TABLE]

and

[TABLE]

where $\tilde{p}_{i}>0$ for $i=1,\dots,j^{*}$, $\sum_{i=1}^{j^{*}}\tilde{p}_{i}<1$, $\tilde{\theta},\tilde{\kappa}$ are some constants in $(0,1)$, and $\tilde{T},\tilde{K},\overline{K}$ are positive constants. Using (A.11) and (A.12), we can obtain the persistence in probability of $(X_{1},\dots,X_{j^{*}})$ in the same manner as (A.7). The proof is complete. ∎

Proof of Theorem 1.1 (iii).

Let $f:\mathbb{R}^{n}_{+}\mapsto\mathbb{R}$ be a continuous function and $\sup_{\mathbf{x}\in\mathbb{R}^{n}_{+}}|f(\mathbf{x})|\leq 1$. Fix $\mathbf{x}_{0}\in\mathbb{R}^{n,\circ}_{+}$. We have to show that

[TABLE]

In part (ii), we have proved that $(X_{1},\dots,X_{j^{*}})$ is persistent in probability. Thus, for any $\varepsilon>0$, there exist $T_{1}>0$ and $H>1$ such that

[TABLE]

For $\delta\geq 0$ define

[TABLE]

Let $\overline{f}=\int_{\mathbb{R}^{n}_{+}}f(\mathbf{x}^{\prime})\pi_{j^{*}}(d\mathbf{x}^{\prime})$. In view of (1.6), there exists $T_{2}>0$ such that

[TABLE]

Since $\mathbf{X}$ is a Markov-Feller process on $\mathbb{R}^{n}_{+}$, we can find a sufficiently small $\delta=\delta(\varepsilon)>0$ such that

[TABLE]

Thus, (A.15) and (A.16) imply

[TABLE]

Since $X_{j^{*}+1},\dots,X_{n}$ converges to [math] almost surely, there exists $T_{3}>T_{1}$ such that

[TABLE]

We deduce from (A.14) and (A.18) that

[TABLE]

For any $t\geq T_{3}+T_{2}$, we have from the Chapman-Kolmogorov equation, (A.17), (A.19) and $|f(\mathbf{x})|\leq 1$ that

[TABLE]

which leads to (A.13). The proof is complete. ∎

Proof of Theorem 1.1 (iv).

If $\Sigma_{j^{*}}$ is positive definite, then by Theorem 1.1 part (i) for $\mathbf{x}\in\mathbb{R}^{(j^{*}),\circ}_{+}$ one has that as $t\to\infty$ the transition probability $P(t,x,\cdot)$ converges in total variation to a unique invariant probability measure $\pi_{j^{*}}$. Moreover, the convergence is uniform in each compact set of $\mathbb{R}^{(j^{*}),\circ}_{+}$ (due to the property of the Lyapunov function constructed in the proof). As a result (1.6) is satisfied and the conclusion follows by part (iii) of Theorem 1.1. ∎