Stochastic Lotka-Volterra food chains

TL;DR

This paper analyzes how environmental stochasticity affects species persistence and extinction in a simple predator-prey food chain modeled by stochastic Lotka-Volterra equations, extending deterministic results to stochastic settings.

Contribution

It provides explicit criteria for species extinction or persistence in stochastic Lotka-Volterra food chains, including cases with degenerate noise and correlations.

Findings

Environmental stochasticity increases extinction risk.

Small environmental fluctuations preserve deterministic persistence.

Adding an apex predator increases extinction likelihood.

Abstract

We study the persistence and extinction of species in a simple food chain that is modelled by a Lotka-Volterra system with environmental stochasticity. There exist sharp results for deterministic Lotka-Volterra systems in the literature but few for their stochastic counterparts. The food chain we analyze consists of one prey and $n-1$ predators. The $j$th predator eats the $j-1$th species and is eaten by the $j+1$th 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 an explicit quantity depending on the interaction coefficients of the system, which species go extinct and which converge to their unique invariant probability measure. Our work can be seen as a natural extension of the deterministic results of Gard and Hallam '79 to a stochastic setting. As one consequence we show that environmental stochasticity makes species more likely to go extinct. However, if the environmental fluctuations are small, persistence in the deterministic setting is preserved in the stochastic system. Our analysis also shows that the addition of a new apex predator makes, as expected, the different species more prone to extinction. Another novelty of our analysis is the fact that we can describe the behavior the system when the noise is degenerate. This is relevant because of the possibility of strong correlations between the effects of the environment on the different species.

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_{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}{x}_n)dt\\ d{x}_n(t)& =x_n(t)(-a_{n0}+a_{n,n-1}{x}_{n-1}(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}\kappa (n)=\kappa & =a_{10}-\frac{a_{11}}{a_{21}}\left[a_{20}+\sum _{j=2}^k\left(\prod _{i=2}^j\frac{a_{2i-2,2i-1}}{a_{2i,2i-1}}\right)a_{2j,0}\right]-\sum _{j=1}^l\left(\prod _{i=1}^j\frac{a_{2i-1,2i}}{a_{2i+1,2i}}\right)a_{2j+1,0}\end{array}$$

$$k=\left\{\begin{array}[]{ll}n/2,&\mbox{if }n~{}\text{even}\\ (n-1)/2,&\mbox{if }n~{}\text{odd}\end{array}\right.$$

$$l=\left\{\begin{array}[]{ll}n/2-1,&\mbox{if }n~{}\text{even}\\ (n-1)/2,&\mbox{if }n~{}\text{odd}.\end{array}\right.$$

$$\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_{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}{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)+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_{23}{x}_3& ={\tilde{a}}_{20}\\ & \text{makebox}[0.0pt]⋮\\ a_{j-1,j-2}{x}_{j-2}-a_{j-1,j}{x}_j& ={\tilde{a}}_{j-1,0}\\ a_{j,j-1}{x}_{j-1}& ={\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\},$$

$${\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\}.$$

$${\mathbb{E}}_{\pi }{X}_i={\int }_{{\mathbb{R}}_{+}^n}{x}_i\pi (d\mathbf{x})=x_i^{(j)}\text{ for }i\le j.$$

$$\begin{array}{rl}\tilde{\kappa }(n)=\tilde{\kappa }& ={\tilde{a}}_{10}-\frac{a_{11}}{a_{21}}\left[{\tilde{a}}_{20}+\sum _{j=2}^k\left(\prod _{i=2}^j\frac{a_{2i-2,2i-1}}{a_{2i,2i-1}}\right){\tilde{a}}_{2j,0}\right]-\sum _{j=1}^l\left(\prod _{i=1}^j\frac{a_{2i-1,2i}}{a_{2i+1,2i}}\right){\tilde{a}}_{2j+1,0},k\ge 2.\end{array}$$

$$k=\left\{\begin{array}[]{ll}n/2,&\mbox{if }n~{}\text{even}\\ (n-1)/2,&\mbox{if }n~{}\text{odd}\end{array}\right.$$

$$l=\left\{\begin{array}[]{ll}n/2-1,&\mbox{if }n~{}\text{even}\\ (n-1)/2,&\mbox{if }n~{}\text{odd}.\end{array}\right.$$

$$\begin{array}{rl}\tilde{\kappa }(n)& =\kappa (n)-\frac{1}{2}{\sigma }_{11}-\frac{a_{11}}{a_{21}}\left[\frac{1}{2}{\sigma }_{22}+\sum _{j=2}^k\left(\prod _{i=2}^j\frac{a_{2i-2,2i-1}}{a_{2i,2i-1}}\right)\frac{{\sigma }_{2j,2j}}{2}\right]-\sum _{j=1}^l\left(\prod _{i=1}^j\frac{a_{2i-1,2i}}{a_{2i+1,2i}}\right)\frac{{\sigma }_{2j+1,2j+1}}{2}\\ & <\kappa (n).\end{array}$$

$$\tilde{\kappa}(n+1)=\left\{\begin{array}[]{ll}\tilde{\kappa}(n)-\left(\prod_{i=1}^{n/2}\frac{a_{2i-1,2i}}{a_{2i+1,2i}}\right)\tilde{a}_{n+1,0},&\mbox{if }n~{}\text{even}\\ \tilde{\kappa}(n)-\frac{a_{11}}{a_{21}}\left(\prod_{i=1}^{(n+1)/2}\frac{a_{2i-2,2i-1}}{a_{2i,2i-1}}\right)\tilde{a}_{n+1,0},&\mbox{if }n~{}\text{odd}\end{array}\right.$$

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

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

$$\underset{t\to \infty }{\lim }\frac{1}{t}{\int }_0^t{\mathbb{E}}_{\mathbf{x}}\left[X_k(s)\right]ds=0.$$

$$\underset{t\to \infty }{\lim }\frac{1}{t}{\int }_0^t{\mathbb{E}}_{\mathbf{x}}\left[X_k(s)\right]ds>0.$$

$$\underset{t\to \infty }{\mathrm{liminf}}\frac{1}{t}{\int }_0^t{\mathbb{P}}_{\mathbf{x}}\left\{(X_1(s),\dots ,X_k(s))\in K_{\epsilon }\right\}ds\ge 1-\epsilon$$

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

$$\underset{t\to \infty }{\lim }\frac{1}{t}{\int }_0^t{\mathbb{E}}_{\mathbf{x}}\left[X_k(s)\right]ds=0,k>j^{\ast }.$$

$$\underset{t\to \infty }{\lim }\frac{1}{t}{\int }_0^t{\mathbb{E}}_{\mathbf{x}}\left[X_k(s)\right]ds=x_k^{(j^{\ast })}>0,k\le j^{\ast }$$

$${\mathbb{P}}_{\mathbf{x}}\left\{\underset{t\to \infty }{\lim }\frac{\ln X_1(t)}{t}={\mathcal{I}}_1={\tilde{a}}_{10},\underset{t\to \infty }{\lim }\frac{\ln X_2(t)}{t}=-{\tilde{a}}_{20}\right\}=1.$$

$${\mathbb{P}}_{\mathbf{x}}\left\{\underset{t\to \infty }{\lim }\frac{\ln X_2(t)}{t}={\mathcal{I}}_2\right\}=1$$

$$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,\circ })$$

$$\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{\}}$$

Full text

As one consequence we show that environmental stochasticity makes species more likely to go extinct. However, if the environmental fluctuations are small, persistence in the deterministic setting is preserved in the stochastic system. Our analysis also shows that the addition of a new apex predator makes, as expected, the different species more prone to extinction.

Another novelty of our analysis is the fact that we can describe the behavior the system when the noise is degenerate. This is relevant because of the possibility of strong correlations between the effects of the environment on the different species.

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

1. 1 Introduction

2. 2 Mathematical framework

3. 3 Proofs

4. 4 Invasion rates

5. 5 Discussion

6. 5.0.1 Future work

1. Introduction

A fundamental problem in ecology is to determine which species go extinct and which persist in a given ecosystem. The complex interactions of the species of the community are sometimes depicted by using food webs. The food web shows the different interaction paths that connect various animals and plants. In general, food webs can be seen as unions of food chains - diagrams that follow one path of interactions. Since species interact in complex ways, food chains and food webs are simplified models of the real world, or caricatures of nature. Nevertheless, one can explain some key properties of an ecosystem by studying these simplified models. To quote [Pim82]

“Like caricatures, though their representation of nature is distorted, there is enough truth to permit a study of some of the features they represent.”

One well-known model for the interaction of a predator and its prey, the Lotka-Volterra model, has been developed by [Lot25] and [Vol28]. Even though the model of [Lot25] and [Vol28] is for one predator and one prey one can easily extend it to food chains of any length. Many of the food chain models studied in the literature are deterministic and of Lotka-Volterra type. [GH79] give criteria for persistence and extinction for Lotka-Volterra food chains while the global stability of nonnegative equilibrium points is studied by [So79, Har79]. More general deterministic food-chains have been analyzed by [Gar80, FS85].

In this paper we analyse simple 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.

Our starting point will be the deterministic Lotka-Volterra system (which has been studied by [GH79])

[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 intra-competition 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$ 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.

Remark 1.1**.**

We note that the above system ignores the intraspecies competition between predators of the same species. The analysis for the setting with intraspecies competituon is different and will be the subject of the future paper [HN17].

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]

Species $i$ goes extinct if

[TABLE]

[GH79] are able to prove using ODE and dynamical systems techniques that the persistence or extinction of (1.1) can be determined using a single parameter which depends on the interaction coefficients $(a_{ij})$. Define

[TABLE]

where

[TABLE]

and

[TABLE]

The following theorem is one of the main results of [GH79].

Theorem 1.1**.**

The food chain modelled by (1.1) is persistent when $\kappa(n)>0$; it is not persistent (that is, some species go extinct) if $\kappa(n)<0.$

In nature, the dynamics of species is inherently stochastic due to the random fluctuations of the environmental factors. The combined effects of biotic interactions and environmental fluctuations are key 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 way of studying this interplay is modelling the populations as discrete or continuous time Markov processes and looking at the long-term behavior of these processes ([Che00, ERSS13, EHS15, LES03, SLS09, SBA11, BEM07, BS09, BHS08, CM10, CCL*+*09]).

In order to take into account environmental fluctuatuons and their effect on the persistence or extinction of species one approach is to study systems that have random environmental perturbations. This can be done by studying stochastic differential equations that arise by adding noise to ordinary differential equations. For compact state spaces [SBA11] provide results for persistence. These results have been generalized by [HN16] 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 [HN16, BL16]) and for general Markov processes (see [Ben16]).

A natural stochastic analogue of (1.1) is the system

[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. We denote by $\boldsymbol{\delta}^{*}$ the probability measure putting all of its mass at the origin $(0,\dots,0)$.

Remark 1.2**.**

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 to which one adds a random noise fluctuation term. See [Tur77, Bra02, Gar88, HNY16, EHS15, ERSS13, SBA11, HN16, 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\geq 2$. For fixed $j\in\{1,\dots,n\}$ we will see that the system

[TABLE]

is intricately related to the stationary distributions of (1.3). It is easy to show that (1.4) has a unique solution, say $\left(x^{(j)}_{1},\dots,x^{(j)}_{j}\right)$. Define

[TABLE]

We will show that when (1.4) has a strictly positive solution the invasion rate of predator $X_{j+1}$ in the habitat of $(X_{1},\dots,X_{j})$ is exactly $\mathcal{I}_{j+1}$. The invasion rate of 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 for $j=1,\dots,n$

[TABLE]

and

[TABLE]

Remark 1.3**.**

We will show in Section 3 that (1.3) has unique strong positive solutions. Furthermore, we show that if $\mathbf{X}:=(X_{1},\dots,X_{n})$ has an invariant probability measure $\pi$ on $\mathbb{R}_{+}^{(j),\circ}$ then

[TABLE]

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

We define the stochastic analogue $\tilde{\kappa}$ of $\kappa$ via

[TABLE]

where

[TABLE]

and

[TABLE]

For notational simplicity we also define $\tilde{\kappa}(1):=\tilde{a}_{10}$.

Remark 1.4**.**

By comparing equations (1.2) and (1.6) one notes that $\tilde{\kappa}(n)$ is what one gets if one does the substitutions $a_{i0}\mapsto\tilde{a}_{i0},i=1,\dots,n$ in the expression for $\kappa(n)$. Note that

[TABLE]

Furthermore, if we add one extra predator $X_{n+1}$ to the food chain then one has

[TABLE]

and in particular $\tilde{\kappa}(n+1)<\tilde{\kappa}(n).$

There are different concepts regarding the persistence and extinction of species. We review some of these below.

Definition 1.1**.**

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

[TABLE]

where $\|\cdot,\cdot\|_{\text{TV}}$ is the total variation norm and $P_{\mathbf{X}}(t,\mathbf{x},\cdot)$ is the transition probability of $(\mathbf{X}(t))_{t\geq 0}$.

Definition 1.2**.**

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

[TABLE]

Definition 1.3**.**

The species $X_{k}$ goes extinct weakly in mean if

[TABLE]

Definition 1.4**.**

The species $X_{k}$ is weakly persistent in mean if

[TABLE]

Definition 1.5**.**

The species $(X_{1},\dots,X_{j^{*}})$ are time-average 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.4) with $j=j^{*}$

We refer the reader to [Sch12] for a discussion of various forms of persistence. Having defined all the necessary concepts we can present the main result of this paper.

Theorem 1.2**.**

Assume that $a_{11}>0$ and $\mathbf{X}(0)=\mathbf{x}\in\mathbb{R}_{+}^{n,\circ}$.

• (i)

If $\tilde{\kappa}(n)>0$ the food chain $\mathbf{X}$ modelled by (1.3) is time-average 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.4) with $j=n$.

Moreover, if $\Sigma$ is positive definite, then the food chain $\mathbf{X}$ is strongly stochastically persistent and converges to its unique invariant probability measure $\pi^{(n)}$ on $\mathbb{R}_{+}^{n,\circ}$.

• (ii)

If there exists $j^{*}<n$ such that $\tilde{\kappa}(j^{*})>0$ and $\tilde{\kappa}(j^{*}+1)\leq 0$ then the predators $(X_{j^{*}+1},\dots,X_{n})$ go weakly extinct in mean, that is

[TABLE]

At the same time, the species $(X_{1},\dots,X_{j^{*}})$ are time-average persistent in probability and

[TABLE]

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

• (iii)

If $n=2$ we can strengthen the extinction results as follows:

• –

If $\mathcal{I}_{1}\leq 0$ then for any $\mathbf{X}(0)=\mathbf{x}\in\mathbb{R}_{+}^{n,\circ}$ we have that $\mathbb{P}_{\mathbf{x}}$-almost surely the randomized occupation measures $\left(\tilde{\Pi}_{t}(\cdot)\right)_{t\geq 0}$ converge weakly to $\boldsymbol{\delta}^{*}$ as $t\to\infty$ and

[TABLE]

In particular, $X_{2}$ goes extinct almost surely exponentially fast. If $\mathcal{I}_{1}<0$ then we also have that $X_{1}$ goes extinct almost surely exponentially fast .

• –

If $\mathcal{I}_{1}>0$ and $\mathcal{I}_{2}<0$ then for any $\mathbf{X}(0)=\mathbf{x}\in\mathbb{R}_{+}^{n,\circ}$ we have that $\mathbb{P}_{\mathbf{x}}$-almost surely the randomized occupation measures $\left(\tilde{\Pi}_{t}(\cdot)\right)_{t\geq 0}$ converge weakly, as $t\to\infty$, to the unique invariant probability measure $\pi^{(1)}$ on $\mathbb{R}^{(1),\circ}_{+}$ and

[TABLE]

so that $X_{2}$ goes extinct almost surely exponentialy fast.

Remark 1.5**.**

We note that by Theorem 1.2 the food chain persists when $\mathcal{I}_{j}>0,j=2,\dots,n$ and goes weakly extinct when $\mathcal{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 when $\mathcal{I}_{j}>0,j=2,\dots,n$ and $\mathcal{I}_{n+1}>0$ and goes extinct when $I_{j^{*}+1}<1$ for some $j^{*}\leq n$. This means that the introduction of an apex predator makes extinction more likely.

Remark 1.6**.**

Having $\Sigma$ positive definite guarantees that the system (1.3) is nondegenerate and that the noise is truly $n$ dimensional. Otherwise the noise is degenerate.

Theorem 1.2 extends previous results on stochastic Lotka-Volterra systems in two or three dimensions (see [LB16, HN16, Rud03]) to an $n$ dimensional setting. We also generalize the work by [Gar84] where the author gives sufficient conditions for stochastic boundedness persistence of stochastic Lotka-Volterra type food web models in bounded regions of state space. We note that the main results by [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 necessary. Moreover, based on which conditions are satisfied, we can say exactly which species persist and which go extinct.

2. Mathematical framework

We rewrite (1.3) as

[TABLE]

where $\mathbf{X}:=(X_{1}(t),\dots,X_{n}(t))$. This is a stochastic process that takes values in $\mathbb{R}_{+}^{n,\circ}:=(0,\infty)^{n}$.

The random normalized occupation measures are defined as

[TABLE]

where $\mathcal{B}(\mathbb{R}_{+}^{n,\circ})$ is the set of all Borel measurable subsets of $\mathbb{R}_{+}^{n,\circ}$. Note that $\widetilde{\Pi}_{t}(B)$ tells us the fraction of time the process $\mathbf{X}$ spends in the set $B$ during the interval $[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}$.

As a result any ergodic measure $\mu\in\mathcal{M}$ must be supported in a subspace of this form. More specifically, 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}_{+}^{\mu}\setminus\mathbb{R}_{+}^{\mu,\circ}$. For the Dirac-measure $\boldsymbol{\delta}^{*}$ concentrated at the origin [math], we have $\mathcal{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 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(\mu_{i}),i=1,\dots,n$ are the Lyapunov exponents of $\mu$.

Remark 2.3**.**

Straightforward computations show that for all $\mathbf{c}\in\mathbb{R}_{+}^{n,\circ}$ and $\gamma_{b}>0$

[TABLE]

As a result Assumption 1.1 of [HN16] is violated and we have to use different methods in this setting.

3. Proofs

Lemma 3.1**.**

Suppose that a sequence $\{\nu_{k},k=1,2,\dots\}$ of probability measures on $\mathbb{R}^{n}_{+}$ converges weakly to $\nu_{0}$. Furthermore, assume that

[TABLE]

If $h:\mathbb{R}^{n}_{+}\mapsto\mathbb{R}$ is a continuous function satisfying

[TABLE]

then

[TABLE]

Proof.

Let $\varepsilon>0$. Since $\lim_{\|\mathbf{x}\|\to\infty}\dfrac{h(x)}{\|\mathbf{x}\|^{m}}=0,$ there is $\ell_{\varepsilon}>0$ such that $|h(x)|\leq\dfrac{\varepsilon\|\mathbf{x}\|^{m}}{H}$ for all $\mathbf{x}$ satisfying $\|\mathbf{x}\|>\ell_{\varepsilon}$. This implies that for any $k$

[TABLE]

Let $\phi_{l}(\cdot):\mathbb{R}^{n}_{+}\to[0,1]$ be a continuous function with compact support satisfying $\phi_{l}(\mathbf{x})=1$ if $\|\mathbf{x}\|\leq l_{\varepsilon}$. One gets that for any $k$ the following sequence of inequalities hold

[TABLE]

Since $\nu_{k}$ converges weakly to $\nu_{0}$ we get

[TABLE]

As a consequence of (3.1) and (3.2)

[TABLE]

The desired result follows by letting $\varepsilon\to 0$.

∎

Lemma 3.2**.**

We have the following claims:

• •

For any $\mathbf{x}=(x_{1},\dots,x_{n})\in\mathbb{R}_{+}^{n}$ the system (2.1) has a unique strong solution with initial value $\mathbf{X}(0)=\mathbf{x}$. The solution satisfies

[TABLE]

• •

The process $\mathbf{X}$ is a Markov-Feller process.

• •

There exist constants $p>0$ and $M_{p}>0$ such that

[TABLE]

• •

For any invariant measure $\mu$ of $\mathbf{X}$, we have

[TABLE]

Proof.

The existence and uniqueness of strong solutions with initial values $\mathbf{x}\in\mathbb{R}_{+}^{n}$ satisfying (3.4) can be shown by standard arguments such as those from [LM09, Theorem 2.1] and [LB16, Lemma 1]. Let

[TABLE]

We have

[TABLE]

If we define $A_{1}=\min_{i=2}^{n}\{a_{i0}\}>0$ and $A_{0}=\sup_{\{x_{1}>0\}}\left\{x_{1}(a_{i0}-a_{11}x_{1})+A_{1}x_{1}\right\}<\infty$ then we can see that for all $\mathbf{x}\in\mathbb{R}_{+}^{n,\circ}$

[TABLE]

Let $p>0$ be sufficiently small such that

[TABLE]

In view of (3.7), (3.8), (3.9) and Itô’s formula, we have

[TABLE]

Sine $f_{i}(\mathbf{x})$, $i=1,\dots,n$ are loally Lipschitz funtions, it follows from (3.10) and [NYZ17, Theorem 5.1] that the process $\mathbf{X}$ is a Feller-Markov process.

Moreover, using (3.10), Dynkin’s formula and a standard argument (see e.g. [LM09, Theorem 3.2] or [HN16, Lemma 2.2]), we can easily obtain that

[TABLE]

which implies (3.5). For any invariant probability measure $\mu$ of $\mathbf{X}$ and $H>0$, it follows from Fatou’s lemma that

[TABLE]

Then letting $H\to\infty$ we obtain (3.6). ∎

Lemma 3.3**.**

Suppose $\mu\in\mathcal{M}$ such that $I_{\mu}=\{n_{1},\dots,n_{k}\}$. Then

[TABLE]

for any $i\in I_{\mu}$. As a consequence, if $\mathbf{X}:=(X_{1},\dots,X_{n})$ has an invariant probability measure $\pi$ on $\mathbb{R}_{+}^{(j),\circ}$ then

[TABLE]

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

Remark 3.1**.**

The intuition behind equation (3.11) 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.

Proof.

Let $\mathbf{X}^{\mu}(t)=(X_{1}^{\mu}(t),\dots,X_{n}^{\mu}(t))$ be the stationary solution whose distribution at any time $t$ is $\mu$. By Itô’s formula, we have

[TABLE]

Since $f_{i}$ is a linear function, it follows from (3.6) that $f_{i}$ is $\mu$-measurable. By the ergodicity of $\mathbf{X}^{\mu}(t)$,

[TABLE]

On the other hand, it is well-known that almost surely

[TABLE]

Combining these limits, we obtain

[TABLE]

almost surely. If $\lambda_{i}(\mu)$ is nonzero, we have with probability 1 that

[TABLE]

Both cases contradict the fact that $\mu(\mathbb{R}^{\mu,\circ}_{+})=1.$ The first assertion is therefore proved.

To prove the second claim, suppose that $\pi$ is an invariant probability measure on $\mathbb{R}^{(j),\circ}_{+}$. Then

[TABLE]

Solving the system $\lambda_{k}(\pi)=0,k=1,\dots,j$ we obtain the desired result. ∎

Lemma 3.4**.**

Suppose we have $\mu\in\mathcal{M}$ with $I_{\mu}=\{n_{1},\dots,n_{k}\}$. Then $I_{\mu}$ must be of the form $\{1,2,\dots,l\}$ for some $1\leq l\leq n$. In other words, for $\mu\in\mathcal{M}$, there exists $l\in\{1,\dots,n\}$ suh that $\mu(\mathbb{R}^{(l),\circ}_{+})=1.$

Proof.

We argue by contradiction. First, suppose that $n_{1}>1$. By Lemma 3.3

[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 Lemma 3.3 leads to

[TABLE]

which is a contradiction. ∎

Lemma 3.5**.**

We have the following assertions.

1. (1)

If $\tilde{\kappa}(n)>0$ then for any $\pi\in\operatorname{Conv}(\mathcal{M})$ one has $\max_{i=1}^{n}\lambda_{i}(\pi)>0$. 2. (2)

If $\tilde{\kappa}(n)\leq 0$ then for any $\mu\in\mathcal{M}$, we have $\lambda_{n}(\mu)\leq 0$. Moreover, there is no invariant probability measure in $\mathbb{R}^{n,\circ}_{+}$.

Proof.

If $\tilde{\kappa}(n)>0$ then $\tilde{a}_{10}>0$. Thus, $\lambda_{1}(\boldsymbol{\delta}^{*})=\tilde{a}_{10}>0.$ In view of Lemma 3.4, other invariant measures in $\mathcal{M}$ must be supported on sets of the form $\mathbb{R}^{(i),\circ}_{+}$ for some $i\in\{1,\dots,n\}$.

Since $\tilde{\kappa}(i)\geq\tilde{\kappa}(n)>0$ for $i=1,\dots,n,$ in order to prove claim (1) of the lemma, we need show that if $\mu\in\mathcal{M}$ is an invariant measure on $\mathbb{R}^{(i),\circ}_{+}$ then $\lambda_{i+1}(\mu)>0.$ By Lemma 3.3 and the definition of $\mathcal{I}_{i}$, we have that $\lambda_{i+1}(\mu)=\mathcal{I}_{i+1}$ which, by Lemma 4.1, has the same sign as $\tilde{\kappa}(i+1)$. Thus, $\lambda_{i+1}(\mu)>0$.

For $\pi\in\operatorname{Conv}(\mathcal{M})$, we can decompose $\pi=\rho_{1}\mu_{i_{1}}+\dots+\rho_{k}\mu_{i_{k}}$ where $0\leq i_{1}<\dots<i_{k}<n$, $\rho_{j}>0$ for $j=1,\dots,k$ and $\mu_{i_{j}}\in\mathcal{M}_{i_{j}}$. Since $i_{j}\geq i_{1}+1$ for $j=2,\dots,k$ we have from Lemma 3.3 that $\lambda_{i_{1}+1}(\mu_{i_{k}})=0$. Thus, $\lambda_{i_{1}+1}(\pi)=\rho_{1}\lambda_{i_{1}+1}(\mu_{i_{1}})=\rho_{1}\mathcal{I}_{i+1}>0$. Part (i) is therefore proved.

Now, suppose that $\tilde{\kappa}(n)\leq 0$. Clearly, if $\mu\in\mathcal{M}$ and $\mu(\mathbb{R}^{(j),\circ}_{+})=1$ for $j<n+1$, then $\lambda_{n}(\mu)=-\tilde{a}_{n0}<0$. If $\mu(\mathbb{R}^{(n-1),\circ}_{+})=1$, it follows from (1.5), (2.2) and (3.12) that $\lambda_{n}(\mu)=\mathcal{I}_{n}\leq 0$ since $\mathcal{I}_{n}$ has the same sign as $\tilde{\kappa}(i+1)$.

Finally, if there is an invariant probability measure $\pi$ in $\mathbb{R}^{n,\circ}_{+}$ then $x^{(n)}_{n}=\int_{\mathbb{R}^{n}_{+}}x_{n}\pi(d\mathbf{x})>0$ is the $n$-th component of the solution to (1.4) with $j$ replaced by $n$. By Lemma 4.1, we must have $\tilde{\kappa}(n)>0$, which completes the proof.

∎

Lemma 3.6**.**

For $i=1,\dots,n$, let $\mathcal{M}_{i}$ be the set of invariant probability measures $\mu$ of $\mathbf{X}$ satisfying $\mu\left(\mathbb{R}^{(i),\circ}_{+}\right)=1$. In particular, let $\mathcal{M}_{0}=\{\boldsymbol{\delta}^{*}\}$. Suppose that $\tilde{\kappa}_{j}>0$ for some $j\in\{1,\dots,n\}$. Then, if $\mathcal{M}_{j}\neq\emptyset$, the family $\mathcal{M}_{j}$ is tight in $\mathbb{R}^{(j),\circ}_{+}$.

Proof.

Suppose that $\mathcal{M}_{j}\neq\emptyset$. Then it is tight in $\mathbb{R}^{(j)}_{+}$ due to (3.6). Suppose that $\mathcal{M}_{j}$ is not tight in $\mathbb{R}^{(j),\circ}_{+}$. Then we can extract a sequence of probability measure $\{\mu_{k},k\in\mathbb{N}\}\subset\mathcal{M}_{j}$, such that $\mu_{k}$ converges weakly to a probability measure $\tilde{\mu}$ satisfying $\mu\left(\partial\mathbb{R}^{(j)}_{+}\right)>0$. Decompose $\tilde{\mu}=\tilde{p}\tilde{\mu}_{1}+(1-\tilde{p})\tilde{\mu}_{2}$ where $\tilde{p}\in(0,1]$ and $\tilde{\mu}_{1}\left(\partial\mathbb{R}^{(j)}_{+}\right)=1,\tilde{\mu}_{2}\left(\mathbb{R}^{(j),\circ}_{+}\right)=1$. Since the $\mu_{k}$’s are invariant probability measures of $\mathbf{X}$ and $\partial\mathbb{R}^{(j)}_{+}$ and $\mathbb{R}^{(j),\circ}_{+}$ are invariant sets, we deduce that $\tilde{\mu}_{1}$ and $\tilde{\mu}_{2}$ are also invariant probability measures of $\mathbf{X}$. Using part (1) of Lemma 3.5 with $n$ replaced by $j$, we have

[TABLE]

In view of Lemma 3.3,

[TABLE]

In particular $\lambda_{i}(\tilde{\mu}_{2})=0$ for $i=1,\dots,j$. Thus,

[TABLE]

On the other hand, it follows from (3.13), (3.6) and Lemma 3.1 that

[TABLE]

which contradicts (3.14). As a result, $\mathcal{M}_{j}$ is tight in $\mathbb{R}^{(j),\circ}_{+}$. ∎

Proof of Theorem 1.2 (i).

In view of Itô’s formula,

[TABLE]

Since the expectations with respect to $\mathbb{E}_{\mathbf{x}}$ of the right hand side terms exist, the expectations $\mathbb{E}_{\mathbf{x}}\ln X_{i}(t)$, $i=1,\dots,n$ also exist. As a result of Lemma 3.2 for all $\mathbf{x}\in\mathbb{R}^{n}_{+}$ we have

[TABLE]

Consider the empirical measure

[TABLE]

In view of Lemma 3.3, for each $\mathbf{x}\in\mathbb{R}^{n}_{+}$, the family $\{\Pi^{\mathbf{x}}_{t},t\geq 0\}$ is tight. The weak limit points of $\{\Pi^{\mathbf{x}}_{t},t\geq 0\}$ are invariant probability measures of $\mathbf{X}$ (see e.g. [EHS15, Proposition 6.4]). Let $\nu$ be any weak limit point of $\{\Pi^{\mathbf{x}}_{t},t\geq 0\}$. We decompose $\nu$ as $\nu=\rho\nu_{1}+(1-\rho)\nu_{2}$ where $\rho\in(0,1)$, $\nu_{1}\in\operatorname{Conv}(\mathcal{M})$ and $\nu_{2}\in\mathcal{M}_{n}$. At this moment, we have not proved that $\mathcal{M}_{n}$ is nonempty, but the decomposition is well-defined because we can let $\rho=1$ if $\mathcal{M}_{n}$ is empty. Suppose that we have a sequence $(t_{\ell},\ell\in\mathbb{N})$ with $t_{\ell}\to\infty$ and that $\left(\Pi^{\mathbf{x}}_{t_{\ell}}\right)_{\ell}$ converges weakly to the probability measure $\nu$. Moreover, using (3.5) and Lemma 3.1, Lemma 3.2, we can show that

[TABLE]

In light of Lemma 3.4, we can write $\nu_{1}$ in the form $\nu_{1}=q_{1}\mu_{k_{1}}+\dots+q_{m}\mu_{k_{1}}$ where $\mu_{k_{i}}\in\mathcal{M}_{i},q_{i}>0$ for $i=1,\dots,m$, and $0\leq k_{1}<\dots<k_{m}<n$.

We have proved that $\lambda_{k_{1}+1}(\mu_{k_{1}})>0$. If $m\geq 2$ then we have from Lemma 3.3 that $\lambda_{j}(\mu_{k_{i}})=0$ for $i=2,\dots,m$ and $j=1,\dots,k_{i}$. In particular, $\lambda_{k_{1}+1}(\mu_{k_{i}})=0$ for $i=2,\dots,m$. Moreover $\lambda_{i}(\nu_{2})=0$ for $i=1,\dots,n$. Thus $\lambda_{k_{1}+1}(\nu)=\rho\lambda_{k_{1}+1}(\nu_{1})=\rho q_{1}\lambda_{k_{1}+1}(\mu_{k_{1}})$. Therefore

[TABLE]

This and (3.16) imply that $\rho=0$. Thus, $\mathcal{M}_{n}$ is nonempty and all the weak limit points of $\{\Pi^{\mathbf{x}}_{t},t\geq 0\}$ belong to $\mathcal{M}_{n}$.

Moreover, since $\mathcal{M}_{n}$ is tight in $\mathbb{R}^{n,\circ}_{+}$ (see Lemma 3.6), we can find for any $\varepsilon>0$ a compact set $K_{\varepsilon}\subset\mathbb{R}^{n,\circ}_{+}$ such that

[TABLE]

If $\Gamma$ is positive definite, the diffusion process $\mathbf{X}$ is nondegenerate in $\mathbb{R}^{n,\circ}_{+}$. As a result the invariant probability measure on $\mathbb{R}^{n,\circ}_{+}$ is unique and the strong stochastic persistence of $\mathbf{X}$ follows from [Kal02, Theorem 20.20]. ∎

Proof of Theorem 1.2 (ii).

Suppose that $\tilde{\kappa}(j^{*})>0$ and $\tilde{\kappa}(j^{*}+1)\leq 0$. By part (2) of Lemma 3.5, if $\mu$ is an ergodic measure of $\mathbf{X}$ then $\lambda_{j}(\mu)<0$ for $j>j^{*}$. Applying Lemma 3.3 and Lemma 3.5, we deduce that

[TABLE]

Suppose $\mathbf{x}\in\mathbb{R}^{n,\circ}_{+}$. Let $\{t_{k}\}$ be any subsequence such that $\lim_{k\to\infty}t_{k}=\infty$ and $\Pi^{\mathbf{x}}_{t_{k}}(d\mathbf{x})$ converges weakly to some measure $\nu\in\operatorname{Conv}(\mathcal{M})$. We decompose $\nu$ as $\nu=\rho\nu_{1}+(1-\rho)\nu_{2}$ where $\rho\in(0,1)$, $\nu_{1}$ is a invariant probability measure on $\partial\mathbb{R}^{(j^{*}),\circ}_{+}$ and $\nu_{2}\in\mathcal{M}_{j^{*}}$. Then, using the arguments from the proof for part (i) of Theorem 1.2 with $n$ replaced by $j^{*}$, we can show that any weak limit points of $\left(\Pi^{\mathbf{x}}_{t_{k}}(d\mathbf{x})\right)_{k}$ for $\mathbf{x}\in\mathbb{R}^{n,\circ}_{+}$ belong to $\mathcal{M}_{j^{*}}$. The time-average persistence in probability of $(X_{1},\dots,X_{j^{*}})$ then follows from Lemma 3.6. Moreover, since

[TABLE]

for $\mathbf{x}\in\mathbb{R}^{n,\circ}_{+}$, we have from Lemma 3.1 and Lemma 3.3 that

[TABLE]

∎

Proof of Theorem 1.2 (iii).

Consider the system

[TABLE]

Define the process $\tilde{X}$ via

[TABLE]

By a comparison argument (see [HNY16, GM94] and [EHS15][Theorem 6.1]), if $X_{1}(0)=\tilde{X}(0)=x\in\mathbb{R}_{+}^{\circ}$ then

[TABLE]

The symptotic behavior of $\tilde{X}$ is wel-known (see [EHS15]). Namely,

• •

If $\mathcal{I}_{1}=\tilde{a}_{10}<0$ then $\mathbb{P}_{x}\left\{\lim_{t\to\infty}\tilde{X}(t)=0\right\}=1,x\in\mathbb{R}_{+}^{\circ}.$

• •

If $\mathcal{I}_{1}=\tilde{a}_{10}=0$ then $\mathbb{P}_{x}\left\{\lim_{t\to\infty}\dfrac{1}{t}\int_{0}^{t}\tilde{X}(s)ds=0\right\}=1,x\in\mathbb{R}_{+}^{\circ}.$

• •

If $\mathcal{I}_{1}=\tilde{a}_{10}>0$ then $\mathbb{P}_{x}\left\{\lim_{t\to\infty}\dfrac{1}{t}\int_{0}^{t}\tilde{X}(s)ds=\dfrac{\tilde{a}_{10}}{a_{11}}\right\}=1.$

By Itô’s formula, equation (3.20) and the asymptotic behavior of $\dfrac{1}{t}\int_{0}^{t}\tilde{X}(s)ds$ we have that if $(X_{1}(0),X_{2}(0))=\mathbf{x}\in\mathbb{R}_{+}^{\circ}$ then

[TABLE]

Thus, if $\mathcal{I}_{1}\leq 0$ or if $\mathcal{I}_{2}=-\tilde{a}_{20}+a_{21}\left[\dfrac{\tilde{a}_{10}}{a_{11}}\right]<0$ then $X_{2}(t)$ goes to 0 as $t\to\infty$ almost surely with respect to $\mathbb{P}_{x},x\in\mathbb{R}_{+}^{\circ}$. This and the fact that $\limsup_{t\to\infty}\frac{1}{t}\int_{0}^{t}X_{1}(s)ds\leq\left[\frac{\tilde{a}_{10}}{a_{11}}\right]\vee 0$ imply that the family of random nomalized occupation measures

[TABLE]

is tight and the family of its weak∗-limits when $t\to\infty$ is a subset of $\operatorname{Conv}(\mathcal{M})$ (see e.g.[EHS15]). Using the fact (see [DNY16] or [HNY16]) that

[TABLE]

and Lemma 3.1, we deduce that if $(t_{k})_{k\in\mathbb{N}}\subset\mathbb{R}_{+}$ is a random sequence going to $\infty$ as $k\to\infty$ and $\tilde{\Pi}_{t_{k}}(\cdot)$ converges weakly to a random measure $\pi$ almost surely, then

[TABLE]

For $(X_{1}(0),X_{2}(0))=\mathbf{x}\in\mathbb{R}_{+}^{2,\circ}$, consider two cases.

• •

If $\mathcal{I}_{1}\leq 0$ then by Theorem 1.2, part (ii), $\boldsymbol{\delta}^{*}$ is the unique invariant probability measure of $(X_{1}(t),X_{2}(t))$. We have from (3.22) that as $t\to\infty$ almost surely with respect to $\mathbb{P}_{\mathbf{x}}$ the occupation measures $\tilde{\Pi}_{t}(\cdot)$ converge weakly to $\boldsymbol{\delta}^{*}$ and

[TABLE]

• •

If $\mathcal{I}_{1}>0$ and $\mathcal{I}_{2}<0$, using (3.21), (which leads to $\limsup_{t\to\infty}\frac{\ln X_{2}(t)}{t}\leq 0$ a.s), (3.22) and a contradiction argument that is similar to the one from the proof of Theorem 1.2, part (i) and (ii), we get that as $t\to\infty$ almost surely with respect to $\mathbb{P}_{\mathbf{x}}$, the occupation measures $\tilde{\Pi}_{t}(\cdot)$ converge weakly to the unique invariant probability measure on $\mathbb{R}^{(1),\circ}_{+}$ and

[TABLE]

∎

4. Invasion rates

We want to analyze the invasion rates $\mathcal{I}_{i},i=1,\dots,n+1$. For this we note by (1.5) that we have to analyze the system (1.4). 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]

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]

and

[TABLE]

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

[TABLE]

and

[TABLE]

Let $l_{i}:=\frac{c_{i}}{-f_{i}}$. Then (4.2) can be written as

[TABLE]

which yields when $j$ is even and $4\leq j\leq n$,

[TABLE]

while if $j$ is odd and $3\leq j\leq n$,

[TABLE]

Next, we want to find a formula for $d_{n}^{\prime}$. Set $e_{i}:=\frac{d_{i+1}}{-f_{i+1}c^{\prime}_{i}}$ and $g_{i}:=\frac{1}{c^{\prime}_{i}}$. Then equation (4.3)can be written as

[TABLE]

One can see that the formula for this recursion is given by

[TABLE]

where we make the convention that $\prod_{i=l}^{m}g_{i}=1$ if $l>m$. For $n\geq 3$ is odd we get

[TABLE]

while for $n\geq 4$ even

[TABLE]

As a result, if $n\geq 3$ is odd we have

[TABLE]

while for $n\geq 4$ even

[TABLE]

Lemma 4.1**.**

The quantities $x_{n}^{(n)}$, $\mathcal{I}_{n}$ and $\tilde{\kappa}(n)$ have the same sign for any $n\geq 1$.

Proof.

It follows from (4.7), (4.8) and (1.6) that

[TABLE]

Thus, $x_{n}^{(n)}$ and $\tilde{\kappa}(n)$ have the same sign. On the other hand, we deduce from (4.3) that

[TABLE]

In view of (4.2), we have $c_{n-1}^{\prime}<0.$ Thus, $x_{n}^{(n)}$ and $\mathcal{I}_{n}$ have the same sign for $n\geq 2$. The case $n=1$ is trivial. ∎

Remark 4.1**.**

The invasion rates are functions of the variances $(\sigma_{ii})_{i=1,\dots,n}$. We note that $\mathcal{I}_{j}(\sigma_{1},\dots,\sigma_{j})$ is strictly decreasing in each variable $\sigma_{i},1\leq i\leq 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 n$) the invasion rates converge to $\widehat{\mathcal{I}}_{i}$, that is $\mathcal{I}_{i}\uparrow\widehat{\mathcal{I}}_{i}$ as $\sigma_{ii}\downarrow 0$, where

[TABLE]

while for $n\geq 4$ even

[TABLE]

Even though our methods do not work in the deterministic setting, the expressions for $\widehat{\mathcal{I}}_{1},\dots,\widehat{\mathcal{I}}_{n}$ give, correctly, the deterministic invasion rates.

5. Discussion

We have analysed the persistence and extinction of species decscribed by a stochastic Lotka-Volterra food chain. Our main result, Theorem 1.2, looks at the setting when there is no intracompetition for the predator species and any species can only interact with the species that are directly above or below it in the food chain. We show that similarly to the deterministic case (see [GH79]) one single factor, $\tilde{\kappa}(n)$, determines which species persist and which go extinct in a weak sense. It is interesting to note that one can recover $\tilde{\kappa}(n)$ from the constant $\kappa(n)$ (which determines the behavior of the deterministic food-chain), by doing the substitutions $a_{10}\mapsto\tilde{a}_{10}=a_{10}-\frac{\sigma_{11}}{2}$ and $a_{j0}\mapsto\tilde{a}_{j0}=a_{j0}+\frac{\sigma_{jj}}{2},j\geq 2$. Therefore, from a persistence/extinction point of view, the effect of the stochastic environment is that it lowers the growth rate of the prey species by one half of the variance of the noise affecting the prey and increases the death rates of all the predators by one half of the variance of the respective noise terms. This shows that in this model environmental noise inhibits the coexistence of species.

For technical reasons we cannot say anything about the speed of convergence to the invariant probability measure. When one or more species go extinct we can only show that they go extinct in a weak sense (other than in dimension $n=2$, where we prove stronger results).

In Section 4 we give explicit expressions for the invasion rates. The invasion rates are closely related to the factors $(\tilde{\kappa}(i),i=1,\dots,n)$ - something that is shown in Lemma 4.1. 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 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 3.4). 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.

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].

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 [HN16] 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.

Of course, one would usually want to prove some stronger results for when the species go extinct and $n>2$. We conjecture the following result holds.

Conjecture 5.1**.**

Assume that $a_{11}>0$, $\Sigma$ is positive definite and $\mathbf{X}(0)=\mathbf{x}\in\mathbb{R}_{+}^{n,\circ}$. If there exists $j^{*}<n$ such that $\tilde{\kappa}(j^{*})>0$ and $\tilde{\kappa}(j^{*}+1)<0$ then the predators $(X_{j^{*}+1},\dots,X_{n})$ go extinct, that is

[TABLE]

At the same time, the normalized occupation measure of $(X_{1},\dots,X_{j^{*}})$ converges weakly to the unique invariant probability measure $\pi^{(j^{*})}$ on $\mathbb{R}_{+}^{(j^{*}),\circ}$.

We are able to prove this conjecture when there exists strictly positive intraspecies competition among the predators. These results will appear in the follow-up paper [HN17].

[GH79, MHP14] are able to analyze the system

[TABLE]

Note that this is exactly (1.1) with $a_{11}=0$. The prey species grows exponentially in the absence of predators. The persistence and extinction of species can still be categorized in the deterministic model (5.1) (see [GH79, Theorem 5]). However, our methods do not work in the stochastic generalization of this setting. We must assume in our proofs that $a_{11}>0$. This assumption is natural from a biological point of view since resources are limited so the prey should not be able to grow without bounds on its own. For results when $a_{11}=0$, in the case when one uses telegraph noise instead of white noise, the reader is referred to [TDHS06] where the authors prove a surprising result: switching between two deterministic two-dimensional predator-prey systems does not lead to persistence or extinction. We expect the stochastic version of (5.1) with white noise to exhibit similar strange properties.

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 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 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 affecting 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.

5.0.1. Future work

There are many possible directions for extending our results or adapting them to different settings. One natural generalization would be to work with more general food chains, not necessarily of Lotka-Volterra type, that have been studied in the deterministic setting (see [Gar80, FS85]) and add environmental fluctuations. We expect that the newly developed methods from [Ben16] and [HN16] will be key when trying to prove results about persistence and extinction of populations in stochastic environments.

A different problem would be the one where environmental fluctuations are modelled by telegraph noise instead of white noise. This would make our system a piecewise deterministic Markov process (PDMP). These processes have recently been studied in a biological context and have offered new insight regarding the competitive exclusion principle (see [BL16]) and the long term behavior of predator-prey communities (see [Cos16]).

One can expect that the functional response of some predator species changes according to the seasonally varying prey availability. A first result in this direction has appeared in [TL16]. However, in [TL16] the authors assume that the lengths of the seasons are constant whereas a more intuitive assumption would be that they are random. We expect to explore this direction in future work.

Our model does not account for population structure. Examples of structured populations can be found by looking at a population in which individuals can live in one of $n$ patches (e.g. fish swimming between basins of a lake or butterflies dispersing between meadows). Dispersion is viewed by many population biologists as an important mechanism for survival. Not only does dispersion allow individuals to escape unfavorable landscapes (due to environmental changes or lack of resources), it also facilitates populations to smooth out local spatio-temporal environmental changes. Lotka-Volterra systems with dispersion in the deterministic setting have been studied in [Has78]. It would be interesting generalize this to a stochastic setting.