#  Hierarchical Model with Allee Effect, Immigration, and Holling Type II
Functional Response

Eddy Kwessi111Corresponding author: Department of Mathematics, Trinity
University, 1 Trinity Place, San Antonio, TX 78212, Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we discuss a hierarchical model, based on a Ricker competition
model. The species considered are competing for resources and may be subject
to an Allee effect due to mate limitation, anti-predator vigilance or
aggression, cooperative predation or resource defense, or social
thermoregulation. The species may be classified into a more dominant species
and less dominant or “wimpy” species or just as a predator and prey. The model
under consideration also has components taking into account immigration in
both species and a more structured Holling type II functional response. Local
and global stability analyses are discussed and simulations are provided. We
also consider demographic stochasticity on the species due to environmental
fluctuations in the form of Wiener processes and we show that there are
conditions under which a global solution exists, stationary distributions
exists, and strong persistence in mean of the species is possible. We also use
Wasserstein distance to show empirically that stochasticity can acts as
bifurcation parameter.

##  1 Introduction

Complex systems in nature often have an either open or hidden hierarchy
between its parts or subparts. Therefore dissections of its structure can help
understand its dynamics. This hierarchy may be dictated by size, individual
strength, group structure or organization. For instance, in a bee colony,
there is queen who is the largest member of the colony by size, males drones,
and under-developed female called workers. Another example can be found among
different species living in the wild, for instance, there is a hierarchy in
strength between lions, hyenas, and say antelopes. In the wild especially,
there is a constant displacement of species, either individually or by group,
or due to research for food or water, due changes in the environment, or due
to mating needs. There is a vast literature on using hierarchical models in
biological systems. We will focus in this manuscript on populations dynamics
in the ecological world. In particular, we are interested studying the
dynamics between species that are subject to immigration and emigration,
subject to an Allee effect, with a more structured functional response. For
self-containment, we recall that an Allee effects (see Allee (1949)) is a
phenomenon in population dynamics where there is a positive correlation
between a population density and its relative growth rate. It is sometime
divided into weak and strong Allee effect, see for instance Hutchings (2015).
The strong Allee effect occurs when a population has a critical density $A$
below which it declines to extinction and above which it increases towards its
carrying capacity $K$. The weak Allee effect occurs when a population lacks
such a critical density, but at lower densities, the population growth rate
arises with increasing densities.

Figure 1:  Represented is the relative capita growth rate as a function of the
population density $x$. The red curve represents no Allee effect since there
is a negative correlation between the per capita growth rate and the density.
The black dashed curve represents the weak Allee effect since at lower
densities $x$, the per capita growth rate is increasing but there is no Allee
threshold $A$. The blue dashed represents the strong Allee effect since there
is positive correlation between the per capita growth rate and the existence
of a Threshold under with the population decreases to extinction.

We also recall in that in ecology, a functional response represents the intake
rate of a species as a function of food density. Holling (1959) proposed three
types of functional responses. Let $x$ and $y$ represents the density of two
populations; then a type I functional response is of the form $h_{1}(x,y)=\mu
xy$, where $\mu$ is a constant; a type II functional response is of the form
$\displaystyle h_{1}(x,y)=\frac{\alpha xy}{1+\alpha\beta xy}$, where $\alpha$
represents the attack rate and $\beta$ is the handling time, that is, the time
spent by say species $x$ searching and processing food obtained from species
$y$. A type III functional response is of the form $\displaystyle
h_{1}(x,y)=\frac{\mu}{1+\alpha e^{-\beta xy}}$, where $\alpha$ and $\beta$ are
as above and $\mu$ is a constant representing a saturation level, that is, a
rate threshold when species densities are high.

Figure 2: Represented are the three functional response types with chosen
parameters.

In ecology, there is already a vast literature on the stability analysis of
deterministic models. For hierarchical model with one, two, or more species,
the interested reader can see for instance Assas et al. (2015b), Assas et al.
(2015a), and Elaydi et al. (2018). Considering that species live in habitats
that are often subject to demographics fluctuations or perturbations, it is
sometime more accurate to consider stochastic models. In this case,
persistence and coexistence of species despite environmental fluctuations are
of particular interest. Papers such as Chesson (2009), Benaïm and Schreiber
(2009), Hening and Li (2021), Hening et al. (2022), the references therein,
and their subsequent iterations are great introduction to the understanding of
the biological motivations and necessary theoretical underpinnings. The
remainder of the paper if organized as follows: In section 2, we discuss the
deterministic model by finding its fixed points, then we discuss its stability
both local and global , and we propose some numerical results. In section 3,
we discuss the stochastic model. In particular, we show existence of a global
solution, we show strong persistence in mean, then we proved the existence of
a stationary distribution, and we propose numerical results. These numerical
results consist of phase space diagrams, histograms, and level curves of
species densities as well a comparative analysis based on the Wasserstein
distance. In section 4, we make some concluding remarks.

##  2 Deterministic Model

In Assas et al. (2015b), Assas et al. (2015a), and Elaydi et al. (2018),
hierarchical models for two and three species were discussed. In particular,
the following model was discussed

$\begin{cases}x_{t+1}&=H_{1}x_{t}+x_{t}e^{\displaystyle
r_{1}-x_{t}-\frac{m_{1}}{1+s_{1}x_{t}}}\vspace{0.125cm}\\\
y_{t+1}&=H_{2}y_{t}+y_{t}e^{\displaystyle
r_{2}-y_{t}-bx_{t}-\frac{m_{2}}{1+s_{2}y_{t}}}\end{cases}\;,$ (2.1)

In this model, $x_{t}$ and $y_{t}$ represent the densities of the two species
$x$ and $y$ under consideration with $x$ being the “stronger” species and $y$
being the weaker one. Moreover, the $H_{i}$’s are the immigration constants,
the $r_{i}$’s are the species’ growth rates, the $m_{i}$’s are the species’
mortality rates, the $s_{i}$’s are the Allee effect constants, and $b$ is a
competition constant. For $i=1,2$, the following assumptions were proposed to
guarantee a strong Allee effect on each species:

  * (i)

$0<H_{i}<1,\leavevmode\nobreak\ \leavevmode\nobreak\
r_{i}s_{i}>1+s_{i}\ln(1-H_{i})$ , and $m_{i}>r_{i}-\ln(1-H_{i})$ .

  * (ii)

$[r_{i}s_{i}-1-s_{i}\ln(1-H_{i})]^{2}>4s_{i}(m_{i}-r_{i}+\ln(1-H_{i}))$ .

Each species is subject to a strong Allee effect induced multiplicatively with
the terms $\displaystyle e^{\frac{m_{1}}{1+s_{1}x_{t}}}$ and $\displaystyle
e^{\frac{m_{2}}{1+s_{2}y_{t}}}$. In this manuscript, we propose to discuss the
following model with immigration, and Allee effect on the waker species, and a
Hollins type II functional response. The choice of a type II rather than a
Type III stems from the fact that type III tends to occurs in experimental
data but is rare in nature.

$\begin{cases}\displaystyle x_{t+1}=x_{t}e^{\displaystyle
r_{1}-x_{t}+dh(x,y)+H_{1}}\\\ y_{t+1}=y_{t}e^{\displaystyle
r_{2}-y_{t}-\frac{m}{a+y_{t}}-h(x,y)+H_{2}}\end{cases}\;.$ (2.2)

with a Holling type II functional response $\displaystyle
h(x,y)=\frac{(b+cx_{t})y_{t}}{1+py_{t}(b+cx_{t})}$ and

$\begin{array}[]{||l||l||l||}\hline\cr\mbox{Parameters}&\mbox{Denomination}&\mbox{Dimension}\\\
\hline\cr r_{1},r_{2}&\mbox{Populations growth rates}&(\mbox{Time})^{-1}\\\
\hline\cr H_{1},H_{2}&\mbox{Immigration/emigration
rates}&(\mbox{Time})^{-1}\\\ \hline\cr a&\mbox{Allee Effect
constant}&\mbox{Biomass}\\\ \hline\cr b&\mbox{Species $x$ attack
rate}&(\mbox{Biomass})^{-1}(\mbox{Time})^{-1}\\\ \hline\cr c&\mbox{Cooperation
intensity between the two species}&(\mbox{Biomass})^{-2}(\mbox{Time})^{-1}\\\
\hline\cr d&\mbox{Conversion coefficient}&\mbox{Dimensionless}\\\ \hline\cr
m&\mbox{Mortality rate of species $y$ due to the Allee
Effect}&\mbox{Time}^{-1}\\\ \hline\cr p&\mbox{Species $x$ handling
time}&\mbox{Time}\\\ \hline\cr\end{array}$

######  Remark 1.

One could also consider a spatial or cluster model for which immigration is
dependent on the distance between the center position $\xi$ of the cluster and
the center $\eta$ from which the species is immigrating from, using a
Laplacian spatial dispersion kernel:

$K(\xi,\eta,\gamma)=\frac{1}{\gamma}e^{-\frac{\left|\xi-\eta\right|}{\gamma}}\;.$

For a given cluster/patch centered a $\eta$, we will assume that the ecosystem
has $N$ clusters each centered at $\xi_{i}$, for $j=1,2,\cdots,N$ and
individuals from each species $x$ and $y$ from these clusters move into the
cluster centered at $\eta$ according to spatial dispersion kernels
$K(\eta,\xi,\gamma)$. We will then have the model.

$\begin{cases}\displaystyle
x_{t+1}(\eta)=x_{t}(\eta)\sum_{j=1}^{N}K(\eta,\xi_{j},\gamma_{1})e^{\displaystyle
r_{1}-x_{t}(\xi_{j})+dh(x_{t}(\xi_{j}),y_{t}(\xi_{j}))}\\\ \displaystyle
y_{t+1}(\eta)=y_{t}(\eta)\sum_{j=1}^{N}K_{2}(\eta,\xi_{j},\gamma_{2})e^{\displaystyle
r_{2}-y_{t}(\xi_{j})-\frac{m}{a+y_{t}(\xi_{j})}-h(x_{t}(\xi_{j}),y_{t}(\xi_{j}))}\end{cases}\;.$
(2.3)

In the presence of multiple preys, say $M$, each with density $y_{t}^{i}$,
($i=1,\cdots M$) for the predator with density $x_{t}$ at time $t$, we could
consider the general model

$\begin{cases}\displaystyle
x_{t+1}(\eta)=x_{t}(\eta)\sum_{j=1}^{N}K(\eta,\xi_{j},\gamma_{1})e^{\displaystyle
r_{1}-x_{t}(\xi_{j})+\sum_{i=1}^{M}d_{i}h(x_{t}(\xi_{j}),y_{t}^{i}(\xi_{j}))}&\\\
\displaystyle
y_{t+1}^{i}(\eta)=y_{t}^{j}(\eta)\sum_{j=1}^{N}K_{2}(\eta,\xi_{j},\gamma_{i})e^{\displaystyle
r_{i+1}-y_{t}^{i}(\xi_{j})-\frac{m_{i}}{a_{i}+y_{t}^{i}(\xi_{j})}-h(x_{t}(\xi_{j}),y_{t}^{i}(\xi_{j}))}&\end{cases}\;.$
(2.4)

###  2.1 Stability Analysis

###  2.2 Fixed Points

Recall that $r_{1},r_{2},a,b,c,d,p,m$ are all positive parameters.
We start by finding the fixed points of the model; they are the origin $(0,0)$
and the intersections of the isoclines of equations

$\displaystyle\displaystyle r_{1}+H_{1}-x+\frac{d(b+cx)y}{1+py(b+cx)}$
$\displaystyle=$ $\displaystyle 0\;.\vspace{0.2cm}$ (2.5)
$\displaystyle\displaystyle
r_{2}+H_{2}-y-\frac{m}{a+y}-\frac{(b+cx)y}{1+py(b+cx)}$ $\displaystyle=$
$\displaystyle 0\;.$ (2.6)

After simplifications, we will have

$\displaystyle(C_{1}):\quad x$ $\displaystyle=$ $\displaystyle
r_{1}+H_{1}+d(r_{2}+H_{2})-dy-\frac{dm}{a+y}\;.$ (2.7)
$\displaystyle(C_{2}):\quad y$ $\displaystyle=$
$\displaystyle\frac{x-(r_{1}+H_{1})}{(b+cx)(-px+d+p(r_{1}+H_{1}))}\;.$ (2.8)

The fixed points will be the intersection between the curve ($C_{2}$) in the
$xy$-plane with the $(C_{1})$ in the $yx$-plane. Things to note:

  * •

$(C_{1})$ has a vertical asymptote $y=-a$ and an asymptote
$x=r_{1}+H_{1}+d(r_{2}+H_{2})-dy$ is the $yx$-plane.

  * •

$(C_{2})$ has two vertical asymptotes $\displaystyle x=-\frac{b}{c}$ and
$x=\frac{d+p(r_{1}+H_{1})}{p}$ and a horizontal asymptote $y=0$ in the
$xy$plane.

In all, we can expect the following fixed points.

  * •

The origin $E_{00}(0,0)$.

  * •

Axial fixed points:

    1. 1.

A predator-free fixed point $E_{10}(0,y)$,

    2. 2.

A prey-free fixed $E_{01}(x,0)$.

  * •

At most two interior fixed points $E_{11}(\alpha_{1x},\alpha_{1y})$ and
$E_{22}(\alpha_{2x},\alpha_{2y})$.

(a) | (b)
---|---
|
(c) | (d)
|
Figure 3: Figure (a) represents the case of no interior fixed point, Figure
(b) represents the case of one axial fixed point, Figure (b) represents the
case of one interior fixed point, and Figure (d) represents the case of two
interior fixed points.

####  2.2.1 Local Stability

Let

$\displaystyle h(x,y)$ $\displaystyle=$
$\displaystyle\frac{(b+cx)y}{1+py(b+cx)}\;,$ $\displaystyle f_{1}(x,y)$
$\displaystyle=$ $\displaystyle e^{\displaystyle r_{1}+H_{1}-x+dh(x,y)}\;,$
$\displaystyle f_{2}(x,y)$ $\displaystyle=$ $\displaystyle e^{\displaystyle
r_{2}+H_{2}-y-\frac{m}{a+y}-h(x,y)}\;,$ $\displaystyle F(x,y)$
$\displaystyle=$
$\displaystyle(F_{1}(x,y),F_{2}(x,y))=(xf_{1}(x,y),yf_{2}(x,y))\;.$

The Jacobian matrix of the system above at any given point $M(x,y)$ is given
as

$J(x,y)=\begin{pmatrix}f_{1}(x,y)\left(1+x\left(-1+\frac{dcy}{(1+py(b+cx))^{2}}\right)\right)&xf_{1}(x,y)\left(\frac{d(b+cx)}{(1+py(b+cx))^{2}}\right)\\\
yf_{2}(x,y)\left(\frac{-cy^{2}}{(1+py(b+cx))^{2}}\right)&f_{2}(x,y)\left(1+y\left(-1+\frac{b+cx}{(1+py(b+cx))^{2}}\right)\right)\end{pmatrix}\;.$

At the origin $E(0,0)$, the Jacobian is

$J(0,0)=\begin{pmatrix}f_{1}(0,0)&0\\\
0&f_{2}(0,0)\end{pmatrix}=\begin{pmatrix}e^{r_{1}+H_{1}}&0\\\
0&e^{\displaystyle r_{2}+H_{2}-\frac{m}{a}}\end{pmatrix}\;.$

Thus the eigenvalues are

$\lambda_{1}=e^{r_{1}+H_{1}}\quad\mbox{and
$\lambda_{2}=e^{r_{2}+H_{2}-\frac{m}{a}}$}\;.$

At the predator-free axial fixed point $E(0,y)$ the Jacobian is

$\begin{array}[]{ll}J(0,y)&=\begin{pmatrix}f_{1}(0,y)&0\\\
f_{2}(0,y)\left(\displaystyle\frac{-cy^{3}}{(1+pby))^{2}}\right)&f_{2}(0,y)\left(1+y\left(-1+\frac{b+}{(1+pby))^{2}}\right)\right)\end{pmatrix}\mbox{}\\\
&\mbox{}\\\ &=\begin{pmatrix}e^{r_{1}+H_{1}}&0\\\
e^{r_{2}+H_{2}-\frac{m}{a+y}}\left(\displaystyle\frac{-cy^{3}}{(1+pby))^{2}}\right)&e^{r_{2}+H_{2}-\frac{m}{a+y}}\left(1+y\left(-1+\frac{b}{(1+pby))^{2}}\right)\right)\end{pmatrix}\;.\end{array}$

Thus the eigenvalues are

$\lambda_{1}=e^{r_{1}+H_{1}}\quad\mbox{and
$\lambda_{2}=e^{r_{2}+H_{2}-\frac{m}{a+y}}\left(1+y\left(-1+\frac{b}{(1+pby))^{2}}\right)\right)$}\;.$

At the prey-free axial fixed point $E(x,0)$ the Jacobian is

$\begin{array}[]{lll}J(x,0)&=\begin{pmatrix}f_{1}(x,0)(1-x)&xdf_{1}(x,0)(b+cx)\\\
0&f_{2}(x,0)\end{pmatrix}&=\begin{pmatrix}e^{r_{1}+H_{1}-x}(1-x)&xd(b+cx)e^{r_{1}+H_{2}-x}\\\
0&e^{r_{2}+H_{2}-\frac{m}{a}}\end{pmatrix}\;.\end{array}$

Thus the eigenvalues are

$\lambda_{1}=e^{r_{1}+H_{1}-x}(1-x)\quad\mbox{and
$\lambda_{2}=e^{r_{2}+H_{2}-\frac{m}{a}}$}\;.$

At an interior fixed point $E(x,y)$ with $x,y>0$, the Jacobian is

$J(x,y)=\begin{pmatrix}1+x\left(-1+\frac{dcy}{(1+py(b+cx))^{2}}\right)&\frac{dx(b+cx)}{(1+py(b+cx))^{2}}\\\
\frac{-cy^{3}}{(1+py(b+cx))^{2}}&1+y\left(-1+\frac{b+cx}{(1+py(b+cx))^{2}}\right)\end{pmatrix}\;,$

since $f_{1}(x,y)=f_{2}(x,y)=1$. Let $D$ be the determinant and $T$ be the
trace of $J(x,y)$ respectively. Then the eigenvalues are

$\displaystyle\lambda_{1}=\frac{T-\sqrt{T^{2}-4D}}{2}\quad\mbox{and
$\displaystyle\lambda_{2}=\frac{T-\sqrt{T^{2}-4D}}{2}$}\;.$

We the have the following result:

######  Theorem 2.

Consider the deterministic system above. Then

  * $\blacktriangleright$

The origin $E_{00}(0,0)$ is locally asymptotically stable if and only if
$r_{1}+H_{1}<0$ and $r_{2}+H_{2}<\frac{m}{a}$.

  * $\blacktriangleright$

The Predator-free equilibrium $E_{01}(0,y)$ is locally asymptotically stable
if and only if $r_{1}+H_{1}<0$ and
$\left|1-y+\frac{yb}{(1+pby)^{2}}\right|<e^{\frac{m}{a}-r_{2}-H_{2}}$.

  * $\blacktriangleright$

The Prey-free equilibrium $E_{10}(0,y)$ is locally asymptotically stable if
and only if $\left|x-1\right|<e^{x-r_{1}-H_{1}}$ and
$r_{2}+H_{2}<\frac{m}{a}$.

  * $\blacktriangleright$

From the Determinant-Trace analysis, we know that

    * –

If $\left|\lambda_{1}\right|<1$ and $\left|\lambda_{2}\right|>1$ or
$\left|\lambda_{1}\right|>1$ and $\left|\lambda_{2}\right|<1$, then $E(x,y)$
is locally a saddle point.

    * –

$\left|\lambda_{1}\right|<1$ and $\left|\lambda_{2}\right|<1$, then $E(x,y)$
is locally asymptotically stable.

    * –

$\left|\lambda_{1}\right|>1$ and $\left|\lambda_{2}\right|>1$, then $E(x,y)$
is locally unstable.

###  2.3 Global Stability

In this section, we discuss global stability condition for the origin and the
interior equilibrium point.

######  Theorem 3.

Put

$\mbox{$\alpha=e^{r_{1}+H_{1}+\frac{d}{p}}$ and $\beta=e^{r_{2}+H_{2}}$}\;.$

  * $\blacktriangleright$

The origin $E_{00}(0,0)$ is globally asymptotically stable if $\alpha^{2}\leq
1$ and $\beta^{2}\leq 1$.

  * $\blacktriangleright$

The interior equilibrium $E(x_{*},y_{*})$ is globally asymptotically stable if
$\alpha^{2}=1$ and $\beta^{2}=1$.

###  2.4 Simulations

In this section, we illustrate the deterministic model for the different types
of equilibrium points. We chose eight starting points with respective $x$ and
$y$ coordinates:
$X_{0}=(0.1,0.1,0.1,10,10,10,2.5,20),Y_{0}=(10,0.5,0.1,0.1,0.5,10,0.1,3)$. The
trajectories are represented by the black, light green, light blue, light red,
light cyan, orange, purple, and magenta colors. The solid dots represent the
fixed points of the model and the solid blue and light black curves represent
the isoclines $(C_{2})$ and $C_{1}$ respectively.

(a) | (b) | (c)
---|---|---
|  |
(d) | (e) |
|  |
Figure 4: Figure (a) show convergence of trajectories to the origin $E_{0}0$,
Figure (b) shows convergence to the predator free axial point $E_{0y}$, Figure
(c) shows convergence to the prey free axial point $E_{x0}$, Figure (d) shows
convergence to the interior fixed point $E_{xy}$ Figure (e) shows two interior
fixed points, one stable and one unstable.

###  2.5 Discussion

We can make the following observation from the study above:
1) The choice of the Holling functional certainly plays a role in this model.
Our motivation for choosing type II rather than type III is that type III is
found in population dynamics if the prey density is assumed constant, which is
hardly the case here.
2) Another observation is that our choice of type II functional response is
different from the model suggested by Holling at inception. However, our
functional response accounts for an important parameter in population dynamics
such as the cooperation constant $c$.
3) Our discussion of global stability is limited to using Lyapounov functional
approach. It is important to point out that tools such monotone maps could
also be used to prove global stability, see Balreira and Luis (2014).
4) It is known that in standard Ricker model with no immigration term, the
intrinsic growth rate parameter is a bifurcation parameter (see for example
Kwessi et al. (2018)), and in some case, the standard Ricker model possesses
deterministic chaos. In the current model, the bifurcations parameters will be
$r_{1}+H_{1}+\frac{d}{p}$ and $r_{2}+H_{2}$ respectively for individual
species and it highly likely that this system also possess deterministic
chaos. Since this is beyond the investigation we are interested in pursuing,
we will not discus it further for sake of brevity.

##  3 Stochastic Model

To take into account environmental fluctuations on the species under
consideration, we let

  * •

$k\in\mathbb{N}$, we let $t_{k}=k\Delta t$ for some $\Delta t>0$.

  * •

We define $X_{k}=(x_{{}_{t_{k}}},y_{{}_{t_{k}}})$ for $k\in\mathbb{N}$.

  * •

We let the initial condition be
$X_{0}=(x_{0},y_{0})\in\mathbb{R}^{+}\times\mathbb{R}^{+}$.

  * •

We consider sequences of random variables $\displaystyle\left\\{N_{j,k}(\Delta
t),j=1,2,k\in\mathbb{N}\right\\}$, normally distributed such that for $j=1,2$
and for all $k\in\mathbb{N}$,

    1. 1.

$E[N_{j,k}(\Delta t)]=0$.

    2. 2.

$E[(N_{j,k}(\Delta t))^{2}]=\sigma_{j}^{2}\Delta t$ with for $j=1,2$.

    3. 3.

$E[(N_{j,k}(\Delta t))^{4}]=o(\Delta t)$.

  * •

We will assume that within the interval $[t_{k},t_{k+1}]$, $X_{k}$ is affected
by random perturbations $(x_{{}_{t_{k}}}N_{1,k}(\Delta
t),y_{{}_{t_{k}}}N_{2,k}(\Delta t))$.

Therefore, for $k\in\mathbb{N}$, we will have

$\begin{cases}x_{{}_{t_{k+1}}}-x_{t_{k}}=x_{{}_{t_{k}}}N_{1,k}(\Delta
t)+x_{t_{k}}\left(e^{\displaystyle
r_{1}-x_{t_{k}}+\frac{d(b+cx_{t_{k}})x_{t_{k}}}{1+py_{t_{k}}(b+cx_{t_{k}})}+H_{1}}\right)\Delta
t\\\ y_{{}_{t_{k+1}}}-y_{t_{k}}=y_{{}_{t_{k}}}N_{2,k}(\Delta
t)+y_{{}_{t_{k}}}\left(e^{\displaystyle
r_{2}-y_{{}_{t_{k}}}-\frac{m}{a+y_{{}_{t_{k}}}}-\frac{(b+cx_{{}_{t_{k}}})y_{{}_{t_{k}}}}{1+py_{{}_{t_{k}}}(b+cx_{{}_{t_{k}}})}+H_{2}}\right)\Delta
t\end{cases}\;,$

or in vector form as

$X_{k+1}-X_{k}=X_{k}F(X_{k})\Delta t+X_{k}N_{k}(\Delta t)\;.$

From Itô’s Calculus (see for instance Durrett (1996), if we put
$W_{k+1,j}(\Delta t)-W_{k,j}(\Delta t):=N_{k,j}(\Delta t)$ and let $\Delta
t\to 0$, then the equation above converges to an autonomous stochastic
differential equation (SDE)

$dX(t)=a(X(t))dt+b(X(t))dW(t)\;,$ (3.1)

with initial condition $X_{0}=(x_{0},y_{0})$ and where

$a(X(t))=\left(\begin{array}[]{l}x(t)e^{\displaystyle
r_{1}-x(t)+\frac{d(b+cx(t))y(t)}{1+py(t)(b+cx(t))}+H_{1}}\\\
y(t)e^{\displaystyle
r_{2}-y(t)-\frac{m}{a+y(t)}-\frac{(b+cx(t))y(t)}{1+py(t)(b+cx(t))}+H_{2}}\end{array}\right)\;,$

and $b(X(t))=[x(t)\sigma_{1},y(t)\sigma_{2}]$ represents the per-capita
magnitude of environmental fluctuations, and $W(t)=[W_{1}(t),W_{2}(t)]$ is a
vector of Wiener processes.

######  Remark 4.

It would be an important question to ask why the stochastic model cannot be
introduced by adding a “stochastic” or “random” term to the deterministic
model or even to randomize the model parameters by assuming they are selected
from specific probability distributions. While it would be a worthwhile effort
theoretically, the best answer to this question can be found in Hening and Li
(2021), Remark 2.1. Because of its importance, let us restate it here for
self-containment: just adding a stochastic fluctuating term to a deterministic
model has some short comings because it does not usually give a mechanism on
how different species are influenced by the environment. Instead, following
the fundamental work by Turelli (1977), we see the SDE models as
“approximations for more realistic, but often analytically intractable,
models”. In particular, SDE’s can be seen as scaling limits, or
approximations, of difference equations.

We will be concerned the existence of global solution, strongly persistence in
mean, and the existence of a stationary distribution.

###  3.1 Existence of global solution

Let us start by recalling the following Theorem on the existence of global
solutions to a stochastic differential equation, see for instance Oksendal
(2014), Theorem 5.2.1, p. 66.

######  Theorem 5.

Let $T>0$ and $A:\mathbb{R}^{n}\times[0,T]\to\mathbb{R}^{n},\quad
B:\mathbb{R}^{n}\times[0,T]\to\mathbb{R}^{n\times m}$ be measurable functions.
Let $X_{0}\in\mathbb{R}^{n}$ be a random variable such that
$\mathbb{E}[\left\|X_{0}\right\|^{2}]<\infty$, where $\left\|\cdot\right\|$ is
a norm in $\mathbb{R}^{n}$. Suppose that for given $t\in[0,T]$ and
$X(t),Y(t)\in\mathbb{R}^{n}$, there exist $K_{1},K_{2}>0$ such that

  1. 1.

Linear growth condition:

$\left\|A(X(t),t)\right\|+\left\|B(X(t),t)\right\|\leq
K_{1}(1+\left\|X(t)\right\|)\;.$

  2. 2.

Local Lipschitz Condition:

$\left\|A(X(t),t)-A(Y(t),t)\right\|+\left\|B(X(t),t)-B(Y(t),t)\right\|\leq
K_{2}\left\|X(t)-Y(t)\right\|\;.$

Then the stochastic differential equation $dX(t)=A(X(t),t)+B(X(t),t)dW(t)$
with initial condition $X(0)=X_{0}$ has a unique solution $X(t)$ such that

$\mathbb{E}\left[\int_{0}^{T}\left\|X(t)\right\|^{2}dt\right]<\infty\;.$

######  Theorem 6.

For any initial value
$X_{0}=(x_{0},y_{0})\in\mathbb{R}^{+}\times\mathbb{R}^{+}$, the above
stochastic system has a unique positive global solution
$X(t)=(x(t),y(t))\in\mathbb{R}^{+}\times\mathbb{R}^{+}$.

###### Proof.

For $n=2$, consider the above non-autonomous stochastic differential equation

$dX(t)=a(X(t))+b(X(t))dW(t)\;.$

Put $X(0)=X_{0}=(x_{0},y_{0})$, for given $x_{0},y_{0}\in\mathbb{R}^{+}$. Let
$T>0$ be arbitrary.
Clearly,
$\displaystyle\mathbb{E}[\left\|X_{0}\right\|^{2}]=\left\|X_{0}\right\|^{2}<\infty$.
Recall that

$a(X(t))=[x(t)f_{1}(x(t),y(t)),y(t)f_{2}(x(t),y(t)];\quad
b(X(t))=[x(t)\sigma_{1},y(t)\sigma_{2}]\;.$

$\displaystyle f_{1}(x(t),y(t))=e^{\displaystyle
r_{1}-x(t)+\frac{d(b+cx(t))y(t)}{1+py(t)(b+cx(t))}+H_{1}}\leq
e^{r_{1}+H_{1}+\frac{d}{p}}$.
$\displaystyle f_{2}(x(t),y(t))=e^{\displaystyle
r_{2}-y(t)-\frac{m}{a+y(t)}-\frac{(b+cx(t))y(t)}{1+py(t)(b+cx(t))}+H_{2}}\leq
e^{r_{2}+H_{2}}$.
Using the Mean Value Theorem, we can choose

$K_{1}=K_{2}=\max\left\\{e^{r_{1}+H_{1}+\frac{d}{p}}+\sigma_{1},e^{r_{2}+H_{2}}+\sigma_{2}\right\\}\;.$

Since this is true for any $T>0$, then we conclude that the solution exists
and is global.

∎

###  3.2 Strong Persistence in Mean

Let us start by recalling the notion of strong persistence in mean and an
important lemma used to prove that a stochastic differential equation is
strongly persistent in mean.

######  Definition 7.

Let $X(t)$ be the solution to a stochastic differential equation. Suppose that
for all $t>0$, the normalized occupational measure or mean satisfies

$\left\langle X(t)\right\rangle=\frac{1}{t}\int_{0}^{t}X(u)du<\infty\;.$

$X(t)$ is said to be strongly persistent in mean if

$\lim_{t\to\infty}\inf\left\langle X(t)\right\rangle>0\;.$

To prove that the solution of a stochastic differential equation is strongly
persistent in mean, the following result is often used.

######  Lemma 8.

Let $W_{i}$ be a Wiener Process for all $1\leq i\leq n$. If
$X(t)\in(C[0,\infty],\mathbb{R}^{+})$ and there exist positive constants
$\gamma_{1},\gamma_{2}$ and $T>0$ such that

$\displaystyle\ln(X(t))\geq\gamma_{1}t-\gamma_{2}\int_{0}^{t}X(u)du+\sum_{i=1}^{n}\sigma_{i}W_{i}(t),\quad\mbox{\bf
for $t\geq T$}\;,$

then by the Strong Law of Large numbers,

$\lim_{t\to\infty}\inf\left\langle
X(t)\right\rangle\geq\frac{\gamma_{1}}{\gamma_{2}},\quad\mbox{almost
surely}\;.$

Now we can state the result on mean persistence pertaining to our system:

######  Theorem 9.

Let $X(t)$ be the solution of the stochastic differential equation (3.1). If

$\frac{1}{2}\max\left\\{\sigma_{1}^{2}e^{-r_{1}-H_{1}},\sigma_{2}^{2}e^{\frac{m}{a}+\frac{1}{p}-r_{2}-H_{2}}\right\\}<1\;,$

then $X(t)$ is strongly persistent in mean.

######  Remark 10.

The above theorem can be proved differently. Indeed, the key of our proof is
showing that the quantity
$\displaystyle\lambda_{i}(\mu):=\int_{0}^{\infty}\left(f_{i}(X(t))-\frac{\sigma_{i}^{2}}{2}\right)\mu(dX)>0$,
where $\mu$ is the Lebesgue measure. This is referred to in Hening et al.
(2022) as the external Lyapunov exponent, which determines the infinitesimal
per-capita growth of species not supported by the measure $\mu$. In that
paper, it is shown under certain conditions (see section 3.3 below) that if
$\lambda_{i}(\delta^{*})>0$ (where $\delta^{*}$ is the Dirac measure
concentrated at the origin), then $X$ is strongly stochastically persistent.

###  3.3 Existence of Stationary Distributions

We are now concerned with the existence of a stationary distribution for our
model. First, let us define the notion of stationary distribution and strongly
stochastic persistence along the lines of Hening and Li (2021).

######  Definition 11.

The probability measure $\xi$ is an invariant probability measure for a
process (or a stochastic differential equation with solution) $X(t)$ if,
whenever $X(0)=x$ has distribution $\xi$ , then for any time $t\geq 0$, the
distribution of $X(t)$ is given by $\xi$.

######  Definition 12.

Let $\Sigma$ be a $\sigma$-algebra on
$\mathbb{R}^{n,\circ}_{+}:=(0,\infty)^{n}$ and let $E\in\Sigma$. We define the
total variation norm as

$\displaystyle\left\|\mu(E)\right\|_{TV}=\sup_{\pi}\left\\{\sum_{A\in\pi}\left|\mu(A)\right|,\quad\mbox{$\pi$
is a countable disjoint partition of $E$}\right\\}\;.$

A process $X(t)$ with ${\bf x}=X(0)$ is said to be strongly stochastic
persistent if it has a unique invariant probability measure $\xi_{0}$ defined
on $\mathbb{R}^{n,\circ}_{+}$ and

$\lim_{t\to\infty}\left\|\mathbb{P}(t,X,\cdot)-\xi_{0}(\cdot)\right\|_{TV}=0,\quad\mbox{for
${\bf x}\in\mathbb{R}^{n,\circ}_{+}$}\;.$

We can now state the result on stationary measures related to our model.

######  Theorem 13.

Let $X(t)$ be the solution of the stochastic differential equation (3.1). If

$\frac{1}{2}\max\left\\{\sigma_{1}^{2}e^{-r_{1}-H_{1}},\sigma_{2}^{2}e^{\frac{m}{a}+\frac{1}{p}-r_{2}-H_{2}}\right\\}<1\;,$

then the stochastic differential equation above has a unique stationary
distribution $\xi_{0}$ with support on $\mathbb{R}^{2,\circ}_{+}$ and

$\lim_{t\to\infty}\left\|\mathbb{P}(t,X,\cdot)-\xi_{0}(\cdot)\right\|_{TV}=0,\quad\mbox{for
${\bf x}\in\mathbb{R}^{2,\circ}_{+}$}\;.$

######  Remark 14.

The proof relies of checking the key assumptions of Theorem 2.1 in Hening and
Li (2021). This theorem is quite strong in that it proves strong stochastic
persistence, which is a stronger notion than persistence in probability and
almost sure persistence in probability. Another advantage of Theorem 2.1 is
that it is easy to extend it to multiple species. Often, to show existence of
a stationary distribution, one would show almost sure persistence as in the
proposition below. However, it is not easier to check the second part of the
proposition, which technically amounts to proving Theorem 2.1 in Hening and Li
(2021) from scratch.

######  Proposition 15.

The Markov process $X(t)$ has a unique stationary distribution $\xi$ if there
is a bounded domain $V\in\mathbb{R}^{n}$ with a regular boundary $\Sigma$ for
which the following are hold true:

  1. 1.

In the domain $V$ and some of its neighborhood, the smallest eigenvalue of the
associated diffusion matrix $\Gamma(x)$ is far from zero.

  2. 2.

If ${\bf x}\in\mathbb{R}^{n}\setminus V$, the mean time $\tau({\bf x})$
required for any path emerging from ${\bf x}$ reaching the set $V$ is finite
and $\displaystyle\sup_{{\bf x}\in U}\;\mathbb{E}[\tau({\bf x})]<\infty$ for
every compact set $U\subset\mathbb{R}^{n}$, then for any integrable function
$f$ with respect to the measure $\xi$,

$\mathbb{P}\left(\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}f(X(t))dt=\int_{\mathbb{R}^{n}}f(x)\xi(dx)\right)=1\quad\mbox{for
all $x\in\mathbb{R}^{n}$}\;.$

###  3.4 Simulations

To simulate the behavior of the stochastic model on an interval $[0,T]$ for a
given $T>0$, we will use the following algorithm due to Milstein (1975).

  1. 1.

Select an initial condition $X_{0}=(x_{0},y_{0})$

  2. 2.

Select and integer $N>1$.

  3. 3.

Partition the interval $[0,T]$ as $0=t_{0}<t_{1}<t_{2}<\cdots<t_{N}=T$ with
$t_{k}=k\Delta t$ with $\displaystyle\Delta t=\frac{T}{N}$.

  4. 4.

Use Milstein Higher Order Scheme, that is, for $0\leq k\leq N$.

$\begin{cases}x_{k+1}=&x_{k}+x_{k}f_{1}(x_{k},y_{k})\Delta
t+x_{k}\sigma_{1}\Delta W_{1,k}\sqrt{\Delta
t}+x_{k}\frac{\sigma_{1}^{2}}{2}[(\Delta W_{1,k})^{2}-\Delta t]\\\
y_{k+1}=&y_{k}+y_{k}f_{2}(x_{k},y_{k})\Delta t+y_{k}\sigma_{2}\Delta
W_{2,k}\sqrt{\Delta t}+y_{k}\frac{\sigma_{2}^{2}}{2}[(\Delta
W_{2,k})^{2}-\Delta t]\end{cases}\;.$

where $\Delta W_{i,k}=W_{k+1}-W_{k}$ are independent normally distributed
random variables with zero mean and variance $\Delta t$.

In all figures below, we consider $N=100,T=1$, and the model parameters
corresponding to a two interior fixed points for the deterministic model:
$r_{1}=0.5,r_{2}=0.1,H_{1}=1.2,H_{2}=1.8,m=0.4,a=0.1,b=0.9,c=0.01,d=0.3,p=0.1$.
In the first panel on the right, we plot the phase space diagram of the
deterministic model overlayed with 100 random trajectories with five initial
points with coordinates given by
$X_{0}=(0.25,0.25,2.75,2.75,2.75),Y_{0}=(2,0.25,2,0.25,0.1)$. The
deterministic trajectories are represented by the solid thick red curves. The
second and fifth panels are histograms of stationary distributions of prey and
predator respectively. The third and fourth panels represent the levels curves
and three dimensional representation of the stationary distribution to have
more perspective. In the third plots, the black lines represent respectively
the estimated sample mean of densities of both predator and prey. What they
show is how stochasticity shifts the stable fixed points of deterministic
models.

####  3.4.1 Low stochasticity on both species

|  |
---|---|---
|  |

Figure 5: In this case, we let $\sigma_{1}=\sigma_{2}=0.1$. Unsurprisingly,
the trajectories are similar to that of the deterministic model, with
densities close to symmetric.

####  3.4.2 Low stochasticity on predator, high on prey

|  |
---|---|---
|  |

Figure 6: In this case, we let $\sigma_{1}=6$ and $\sigma_{2}=0.1$. In this
case, the predator survives and the prey quickly approaches limits of
extinction.

####  3.4.3 High stochasticity on predator, low on prey

|  |
---|---|---
|  |

Figure 7: In this case, we let $\sigma_{1}=6$ and $\sigma_{2}=0.1$. In this
case, the predator is persistent for trajectories above the deterministic
unstable fixed point.

####  3.4.4 High stochasticity on both species

|  |
---|---|---
|  |

Figure 8: In this case, we let $\sigma_{1}=6$ and $\sigma_{2}=6$. In this
case, both species will eventually goes extinct, though the species will be
extinct quicker than the predator.

To understand the effect of stochasticity on the densities of the species for
different system parameters, one could track the change in density of predator
and prey from successive increments of parameters $\sigma_{1}$ and
$\sigma_{2}$. We can then assess the Wasserstein distance (see Villani (2008))
between the these consecutive distribution to assess how similar or dissimilar
they are. This would provide a threshold for stochasticity. At the same time,
one could also track the average densities of predator and prey. Let $p>1$ be
a real number. Let us recall that the Wasserstein distance
$W(\theta_{1},\theta_{2})$ on a metric space $(\mathscr{M},D)$ between two
probability measures $\theta_{1}$ and $\theta_{2}$ with joint probability
measure $\theta$ defined on a set $\Gamma(\theta_{1},\theta_{2})$ is given as

$W_{p}(\theta_{1},\theta_{2})=\left(\inf_{\theta\in\Gamma(\theta_{1},\theta_{2})}\left\\{\int_{\mathscr{M}\times\mathscr{M}}\theta(u,v)[D(u,v)]^{p}dudv\right\\}\right)^{\frac{1}{p}}\;.$

For practical purposes, distributions $\theta_{1}$ and $\theta_{2}$ will be
taken as the empirical probability measures with samples
$U=(U_{(1)},U_{(2)},\cdots,U_{(K)})$ and $V=(V_{(1)},V_{(2)},\cdots,V_{(K)})$
from order statistics. Therefore, an estimator of
$W_{p}(\theta_{1},\theta_{2})$ is

$\widetilde{W}_{p}(\theta_{1},\theta_{2})=\left(\sum_{i=1}^{K}\left|U_{(i)}-V_{(i)}\right|^{p}\right)^{\frac{1}{p}}\;.$

In our case, we will use $p=2$. For the predator, we use the samples
$U_{1}=X(\sigma_{1,i})$ and $V_{1}=X(\sigma_{1,i+1})$ corresponding to a given
stochastic parameter $\sigma_{1,i}\in[0.1,8]$ for $i=1,2,\cdots,K$. Likewise,
for the prey we use the samples $U_{2}=Y(\sigma_{2,i})$ and
$V_{2}=Y(\sigma_{2,i+1})$, for $\sigma_{2,i}\in[0.1,8]$. In the figures below,
we chose the starting point of the trajectories to be $x_{0}=5,y_{0}=1$ and
$K=80$. For $i=1,\cdots,80$. The parameters
$\sigma_{1,\cdot}=\sigma_{2,\cdot}=\sigma$ will be referred to as
stochasticity. The left panels represent the Wasserstein distance between
consecutive predator and prey distributions by stochasticity and the right
panels stochasticity versus average population densities.

Figure 9:  Here,
$r_{1}=0.5,r_{2}=0.1,H_{1}=1.2,H_{2}=1.8,m=0.4,a=0.1,b=0.9,c=0.01;d=0.3;p=0.1$.
corresponding to two interior fixed points in the deterministic case.  Figure
10:
$r_{1}=0.1,r_{2}=0.1,H_{1}=0.9,H_{2}=2,m=0.1,a=0.05,b=0.9,c=0.01,d=2.5,p=1.5$
corresponding to a single interior fixed point in the deterministic case.
Figure 11:
$r_{1}=0.9;,r_{2}=0.1,H_{1}=-1,H_{2}=3.5,m=0.1,a=0.5,b=0.2,c=0.4,d=2.5,p=0.1$
corresponding to the predator-free axial fixed point in the deterministic
case.  Figure 12:
$r_{1}=0.1,r_{2}=0.8,H_{1}=0.5,H_{2}=-.9,m=0.8,a=1.2,b=0.4,c=0.6,d=1,p=10$
corresponding to the prey-free axial fixed point in the deterministic case.
Figure 13:
$r_{1}=0.1,r_{2}=0.8,H_{1}=-0.1,H_{2}=-1,m=0.9,a=0.4,b=0.5,c=1,d=2.2,p=0.2$.
corresponding to no interior fixed point in the deterministic case.

##  4 Discussion

  1. 1.

It is interesting to note the similarities in both Wasserstein distance plots
and average population densities plots.

  2. 2.

Indeed, in Figures 9–13, there are critical values for stochasticity from
which quantitatively, the average densities of both predator and prey change.
They are represented by the solid vertical lines, red for the predator and
blue for the prey. This is matched in Wasserstein plots with a drop in
distance.

  3. 3.

From a purely bifurcation point of view, this drop suggests that stochasticity
seems to produce a supercritical pitchfork bifurcation in both species.

  4. 4.

From a purely statistical point of view, the Wasserstein distance obtained is
clearly unimodal and right-skewed when a bifurcation occurs with the right
skew being clearly zero.

  5. 5.

Another important takeaway is that persistence of the prey is clearly seen
through the Wasserstein distance, see Figure 11. In fact, the Wasserstein is
unimodal with values in the same range except for a single peak due to the
disappearance of the predator.

  6. 6.

Figure 12 is interesting in that it shows that while the predator persists for
while under environmental stochasticity, it eventually goes extinct due to
lack of prey. On the other, the prey quickly goes extinct due to both
predation and environmental fluctuations.

  7. 7.

Is important to note that the simulations were done with the parameters of the
deterministic case corresponding to the two interior fixed points. Similar
observations as above can be drawn from the remaining cases.

  8. 8.

We observe that the per-capita magnitude of the environmental fluctuations are
just linear functions of the populations densities. More structured per-capita
magnitude functions can be considered, especially polynomial of higher degree
without fundamentally altering the conclusions we obtained here.

Our simulations are just an initial proof of concept that the Wasserstein
distance may be used successfully to study bifurcation in stochastic
environments. However, further investigations, both theoretically and
practically are needed establish a more solid understanding. Since the
Wasserstein distance is easy to calculate, it may be worthwhile to see if it
can be used to empirically assess chaotic events as an alternative to
calculating the largest Lyapunov exponent.

## References

  * Allee (1949) W. C. Allee.  _Principles of Animal Ecology_.  W.B. Saounders Co., Philadelphia, 1949.  ISBN 0-7216-1120-6.
  * Assas et al. (2015a) L. Assas, B. Dennis, S. Elaydi, E. Kwessi, and G. Livadiotis.  Hierarchical competition models with the allee effect ii: the case of immigration.  _Journal of Biological Dynamics_ , 9(1):288–316, 2015a.
  * Assas et al. (2015b) L. Assas, S. Elaydi, E. Kwessi, G. Livadiotis, and D. Ribble.  Hierarchical competition models with allee effects.  _Journal of Biological Dynamics_ , 9:32–44, 2015b.
  * Balreira and Luis (2014) S. Balreira, E. Elaydi and R. Luis.  Global stability of higher dimensional monotone maps.  _Journal of Difference Equations and Applications_ , 23(12):2037–2071, 2014.
  * Benaïm and Schreiber (2009) P. Benaïm and S. J. Schreiber.  Persistence of structured populations in random environments.  _Theor. Pop. Biol._ , 76(1):19–34, 2009.
  * Chesson (2009) P. Chesson.  General theory of competitive coexistence in spatially-varying environments.  _Theor. Popul. Biol._ , 58(3):211–237, 2009.
  * Durrett (1996) R. Durrett.  _Stochastic Calculus. A practical Introduction_.  CRC Press, 1996.  doi: https://doi.org/10.1201/9780203738283.
  * Elaydi et al. (2018) S. Elaydi, E. Kwessi, and G. Livadiotis.  Hierarchical competition models with allee effect iii: Multispecies.  _Journal of Biological Dynamics_ , 2018.  doi: 10.1080/17513758.2018.1439537.
  * Hening and Li (2021) A. Hening and Y. Li.  Stationary distributions of persistent ecological systems.  _J. of Math. Biol._ , 82(64), 2021.  doi: https://doi.org/10.1007/s00285-021-01613-2.
  * Hening et al. (2022) A. Hening, D. H. Nguyen, and S. J. Schreiber.  A classification of th dynamics of thre-diemnsional stochastic ecological systems.  _The Annals of Applied Probability_ , 32(2):893–931, 2022.
  * Holling (1959) C. S Holling.  Some characteristics of simple types of predation and parasitism.  _The Canadian Entomologist_ , 91(7):385–395, 1959\.
  * Hutchings (2015) J. A. Hutchings.  Thresholds for impaired species recovery.  _Proceedings of the Royal Society B_ , 282:20150654, 2015\.  doi: 10.1098/rspb.2015.0654.
  * Kwessi et al. (2018) E. Kwessi, S. Elaydi, B. Dennis, and G. Livadiotis.  Nearly exact discretization of single species population models.  _Natural Resource Modeling_ , 2018.  doi: 10.1111/nrm.12167.
  * Milstein (1975) G. N. Milstein.  Approximate integration of stochastic differential equations.  _Theory of Probability and its Applications_ , 19(3):557–000, 1975.
  * Oksendal (2014) B. Oksendal.  _Stochastic Differential Equations: An Introduction with Applications_.  (Universitext). Springer, 2014.
  * Turelli (1977) M. Turelli.  Random environments and stochastic calculus.  _Theor. Popul. Biol_ , 12(2):140–178, 1977.
  * Villani (2008) C. Villani.  _Optimal Transport, Old and New_.  Springer, 2008.  ISBN 978-3-540-71050-9.

##  5 Appendix

###  5.1 Proof of Theorem 3

###### Proof.

Let

$F(x,y)=\left(xe^{\displaystyle
r_{1}-x+\frac{d(b+cx)y}{1+py(b+cx)}+H_{1}},\quad ye^{\displaystyle
r_{2}-y-\frac{m}{a+y}-\frac{(b+cx)y}{1+py(b+cx)}+H_{2}}\right)\;.$

Consider the Lyapunov function $L(x,y)=x^{2}+y^{2}$.
Then

$\displaystyle L(F(x,y))-L(x,y)$ $\displaystyle=$ $\displaystyle
x^{2}\left[e^{\displaystyle
2\left(r_{1}+H_{1}-x+\frac{d(b+cx)y}{1+py(b+cx)}\right)}-1\right]$
$\displaystyle+$ $\displaystyle y^{2}\left[e^{\displaystyle
2\left(r_{2}+H_{2}-\frac{m}{a+y}-\frac{(b+cx)y}{1+py(b+cx)}\right)}-1\right]$
$\displaystyle\leq$ $\displaystyle x^{2}\left[e^{\displaystyle
2\left(r_{1}+H_{1}+\frac{db}{p}\right)}-1\right]+y^{2}\left[e^{\displaystyle
2\left(r_{2}+H_{2}\right)}-1\right]\;.$

Thus, $L(F(x,y))-L(x,y)\leq 0$ if $\displaystyle r_{1}+H_{1}+\frac{d}{p}<0$
and $r_{2}+H_{2}<0$ . This shows that the equilibrium point $(0,0)$ is
globally asymptotic stable. Recall that

$\displaystyle h(x,y)$ $\displaystyle=$
$\displaystyle\frac{(b+cx)y}{1+py(b+cx)}\leq\frac{d}{p}\;,$ $\displaystyle
F(x,y)$ $\displaystyle=$
$\displaystyle(F_{1}(x,y),F_{2}(x,y))=(xf_{1}(x,y),yf_{2}(x,y))\;.$

and

$\begin{cases}f_{1}(x,y)=e^{\displaystyle r_{1}+H_{1}-x+dh(x,y)}\leq
e^{r1+H_{1}+\frac{d}{p}}\\\ f_{2}(x,y)=e^{\displaystyle
r_{2}+H_{2}-y-\frac{m}{a+y}-h(x,y)}\leq e^{r_{2}+H_{2}}\end{cases}$ (5.1)

Consider the Lyapunov function

$L(x,y)=(x-x_{*})^{2}+(y-y_{*})^{2}\;.$

and let $\Delta L=L(F(x,y))-L(x,y)\;.$ Then

$\displaystyle\Delta L$ $\displaystyle=$
$\displaystyle(xf_{1}(x,y)-x_{*})^{2}+(yf_{2}(x,y)-y_{*})^{2}-(x-x_{*})^{2}-(y-y_{*})^{2}$
$\displaystyle\leq$ $\displaystyle(\alpha x-x_{*})^{2}+(\beta
x-y_{*})^{2}-(x-x_{*})^{2}-(y-y_{*})^{2}$ $\displaystyle=$
$\displaystyle(\alpha^{2}-1)x^{2}-2(\alpha-1)x_{*}x+(\beta^{2}-1)y^{2}-2(\beta-1)y_{*}y$
$\displaystyle=$
$\displaystyle(\alpha^{2}-1)\left(x-\frac{1}{\alpha+1}x_{*}\right)^{2}-\frac{\alpha-1}{\alpha+1}x_{*}^{2}$
$\displaystyle+$
$\displaystyle(\beta^{2}-1)\left(y-\frac{1}{\beta+1}y_{*}\right)^{2}-\frac{\beta-1}{\beta+1}y_{*}^{2}$

Let $a_{0}=\max\left\\{\alpha^{2}-1,\beta^{2}-1\right\\}$ and let
$\displaystyle
b_{0}=\frac{\alpha-1}{\alpha+1}x_{*}^{2}+\frac{\beta-1}{\beta+1}y_{*}^{2}$. If
$a_{0}\leq 0$and $x_{*}=y_{*}=0$ then $\Delta L\leq 0$. In particular, if
$a_{0}=0$, then clearly, for any $(x_{*},y_{*})$, we have $\Delta L\leq 0$.
Suppose $a_{0}>0$. Then $b_{0}\geq 0$. Put $r_{0}^{2}=\frac{b_{0}}{a_{0}}$.
Then $\Delta L\leq 0$ if
$\left(x-\frac{1}{\alpha+1}x_{*}\right)^{2}+\left(y-\frac{1}{\beta+1}y_{*}\right)^{2}<r_{0}^{2}$.
That is, $\Delta L\leq 0$ if $(x,y)\in
B\left(\left(\frac{1}{\alpha+1}x_{*},\frac{1}{\beta+1}y_{*}\right),r_{0}\right)$,
hence local stability. ∎

###  5.2 Proof of Theorem 9

###### Proof.

Let $0<\theta<1$. Then $e^{-x}\geq-x+\theta$.
Let $\displaystyle g_{1}(x(t)):=\ln(x(t))$. Then using Itô’s formula, we have

$\displaystyle dg_{1}(x(t))$ $\displaystyle=$
$\displaystyle\left[\frac{\partial g_{1}}{\partial
t}+x(t)f_{1}[x(t),y(t)]\cdot\frac{\partial g_{1}}{\partial
x}+\frac{\partial^{2}g_{1}}{\partial
x^{2}}\cdot\frac{\sigma_{1}^{2}}{2}\right]dt+x(t)\sigma_{1}\frac{\partial
g_{1}}{\partial x}dW_{1}(t)$ $\displaystyle=$
$\displaystyle\left[e^{r_{1}+H_{1}-x(t)+\frac{d(b+cx(t))}{1+py(t)(b+cx(t))}}-\frac{\sigma_{1}^{2}}{2}\right]dt+\sigma_{1}dW_{1}(t)$
$\displaystyle\geq$
$\displaystyle\left[e^{r_{1}+H_{1}-x(t)}-\frac{\sigma_{1}^{2}}{2}\right]dt+\sigma_{1}dW_{1}(t)$
$\displaystyle\geq$
$\displaystyle\left[e^{r_{1}+H_{1}}\left(-x(t)+\theta\right)-\frac{\sigma_{1}^{2}}{2}\right]dt+\sigma_{1}dW_{1}(t)\;.$

Integrating from 0 to $t$, we have

$\displaystyle\ln\left(\frac{x(t)}{x(0)}\right)$ $\displaystyle\geq$
$\displaystyle\left[\theta
e^{r_{1}+H_{1}}-\frac{\sigma_{1}^{2}}{2}\right]t-e^{r_{1}+H_{1}}\left\langle
x(t)\right\rangle+\sigma_{1}W_{1}(t)\;.$

Let $\gamma_{11}=\theta e^{r_{1}+H_{1}}-\frac{1}{2}\sigma_{1}^{2}$ and
$\gamma_{12}=e^{r_{1}+H_{1}}$. Let $\displaystyle g_{2}(y(t)):=\ln(y(t))$.
Then using Itô’s formula, we have

$\displaystyle dg_{2}(y(t))$ $\displaystyle=$
$\displaystyle\left[\frac{\partial g_{1}}{\partial
t}+y(t)f_{2}[x(t),y(t)]\cdot\frac{\partial g_{1}}{\partial
y}+\frac{\partial^{2}g_{1}}{\partial
y^{2}}\cdot\frac{\sigma_{2}^{2}}{2}\right]dt+y(t)\sigma_{1}\frac{\partial
g_{1}}{\partial y}dW_{1}(t)$ $\displaystyle=$
$\displaystyle\left[e^{r_{2}+H_{2}-y(t)-\frac{m}{a+y(t)}-\frac{b+cx(t)}{1+py(t)(b+cx(t))}}-\frac{1}{2}\sigma_{2}^{2}\right]dt+\sigma_{2}dW_{2}(t)$
$\displaystyle\geq$
$\displaystyle\left[e^{r_{2}+H_{2}-y(t)-\frac{m}{a}-\frac{1}{p}}-\frac{1}{2}\sigma_{2}^{2}\right]dt+\sigma_{2}dW_{2}(t)$
$\displaystyle\geq$
$\displaystyle\left[e^{r_{2}+H_{2}-\frac{m}{a}-\frac{1}{p}}(-y(t)+\theta)-\frac{1}{2}\sigma_{2}^{2}\right]dt+\sigma_{2}dW_{2}(t)\;.$

Integrating from 0 to $t$, we have

$\displaystyle\ln\left(\frac{y(t)}{y(0)}\right)$ $\displaystyle\geq$
$\displaystyle\left[\theta
e^{r_{2}+H_{2}-\frac{m}{a}-\frac{1}{p}}-\frac{1}{2}\sigma_{2}^{2}\right]t-e^{r_{2}+H_{2}-\frac{m}{a}-\frac{1}{p}}\left\langle
y(t)\right\rangle+\sigma_{2}W_{2}(t)\;.$

Let $\gamma_{21}=\theta
e^{r_{2}+H_{2}-\frac{m}{a}-\frac{1}{p}}-\frac{1}{2}\sigma_{2}^{2}$ and
$\gamma_{22}=e^{r_{2}+H_{2}-\frac{m}{a}-\frac{1}{p}}$ .
It follows that $\gamma_{11},\gamma_{21}>0$ if
$1>\theta>\frac{1}{2}\max\left\\{\sigma_{1}^{2}e^{-r_{1}-H_{1}},\sigma_{2}^{2}e^{\frac{m}{a}+\frac{1}{p}-r_{2}-H_{2}}\right\\}$.
For such $\theta$, put $\gamma_{1}=(\gamma_{11},\gamma_{12})$ and
$\gamma_{2}=(\gamma_{21},\gamma_{22})$. Let $\sigma=(\sigma_{1},\sigma_{2})$.
Since $X(t)=(x(t),y(t))$ and $W(t)=(W_{1}(t),W_{2}(t))$, we will have in
vector form

$\ln(X(t))\geq\gamma_{1}t-\gamma_{2}\int_{0}^{t}X(u)du+\sigma W(t)\;.$

Then by the Lemma 8 above, we conclude that

$\lim_{t\to\infty}\inf\left\langle
X(t)\right\rangle\geq\frac{\gamma_{1}}{\gamma_{2}}>0\quad\mbox{almost
surely}\;.$

Therefore, $X(t)$ is strongly persistent in mean.

∎

###  5.3 Proof of Theorem 13

Consider the the stochastic differential equation (3.1). It is of the form

$dZ_{i}(t)=Z_{i}(t)f_{i}({\bf Z}(t))dt+Z_{i}(t)g_{i}({\bf
Z}(t))dW_{i}(t),\quad i=1,2\;,$

where $Z_{1}(t)=x(t)$ and $Z_{2}(t)=y(t)$. We observe from Remark 10 that if

$\displaystyle\frac{1}{2}\max\left\\{\sigma_{1}^{2}e^{-r_{1}-H_{1}},\sigma_{2}^{2}e^{\frac{m}{a}+\frac{1}{p}-r_{2}-H_{2}}\right\\}<1\;,$

then $\displaystyle\max_{i=1,2}\lambda(\mu)>0$, for any ergodic invariant
measure $\mu$. For given a ${\bf z}\in\mathbb{R}^{2}$ and according to Theorem
2.1 in Hening and Li (2021), it remains to check the following assumptions

  1. $A_{1}:$

$\mbox{diag}(g_{1}({\bf z}),\cdots,g_{n}({\bf
z}))\Gamma^{T}\Gamma\mbox{diag}(g_{1}({\bf z}),\cdots,g_{n}({\bf z}))$ is a
positive definite matrix.

  2. $A_{2}:$

$f_{i}(\cdot),g_{i}(\cdot):\mathbb{R}^{2,\circ}_{+}\to\mathbb{R}$ are locally
Lipschitz functions, for $i=1,2$.

  3. $A_{3}:$

There exist ${\bf c}=(c_{1},c_{2})\in\mathbb{R}^{2,\circ}_{+}$ and
$\gamma_{b}>0$ such that

$\limsup_{\left\|z\right\|\to\infty}\left[\frac{\sum_{i=1}^{2}c_{i}z_{i}f_{i}({\bf
z})}{1+\sum_{i=1}^{2}c_{i}z_{i}}-\frac{1}{2}\frac{\sum_{i,j=1}^{2}\sigma_{ij}c_{i}c_{j}z_{i}z_{j}g_{i}({\bf
z})g_{j}({\bf
z})}{\left(1+\sum_{i=1}^{n}c_{i}z_{i}\right)^{2}}+\gamma_{b}\left(1+\sum_{i=1}^{2}\left(\left|f_{i}({\bf
z})+g_{i}^{2}({\bf z})\right|\right)\right)\right]<0\;.$

We note that in our case, $g_{i}({\bf Z}(t))=1$ for $i=1,2$.
For $A_{1}$, we have

$\mbox{diag}(g_{1}({\bf z}),\cdots,g_{n}({\bf
z}))\Gamma^{T}\Gamma\mbox{diag}(g_{1}({\bf z}),\cdots,g_{n}({\bf
z}))=\begin{pmatrix}\sigma_{1}&0\\\ 0&\sigma_{2}\end{pmatrix}\;,$

which is a positive definite matrix since $\sigma_{i}>0$ for $i=1,2$.
For $A_{2}$, the $f_{i}(\cdot)$’s are locally continuously differentiable
functions with bounded derivative, therefore they are locally Lipschitz
functions.
$A_{3}$ requires a little bit of work. We observe from (5.1), that given ${\bf
z}\in\mathbb{R}^{2,\circ}_{+}$, there exist $k_{1},k_{2}>0$ such $f_{i}({\bf
z})\leq k_{i}$. Let $k=\max\left\\{k_{1},k_{2}\right\\}$. Then we have for all
${\bf c}=(c_{1},c_{2})\in\mathbb{R}^{2,\circ}_{+}$

$\frac{\sum_{i=1}^{2}c_{i}z_{i}f_{i}({\bf z})}{1+\sum_{i=1}^{2}c_{i}z_{i}}\leq
k\frac{\sum_{i=1}^{2}c_{i}z_{i}}{\sum_{i=1}^{2}c_{i}z_{i}}\leq k\;.$

Similarly,

$1+\sum_{i=1}^{2}\left(\left|f_{i}({\bf z})+g_{i}^{2}({\bf
z})\right|\right)\leq 1+k+2=3+k\;.$

Also,

$\sum_{i,j=1}^{2}\sigma_{ij}c_{i}c_{j}z_{i}z_{j}g_{i}({\bf z})g_{j}({\bf
z})=\sum_{i=1}^{2}\sigma_{i}^{2}c_{i}^{2}z_{i}^{2}\;.$

Hence $A_{3}$ will be satisfied if we can find $\gamma_{b}>0$ and ${\bf
c}=(c_{1},c_{2})\in\mathbb{R}^{2,\circ}_{+}$ such that

$k+\gamma_{b}(3+k)\leq\frac{1}{2}\frac{\sum_{i=1}^{2}\sigma_{i}^{2}c_{i}^{2}z_{i}^{2}}{\left(1+\sum_{i=1}^{2}c_{i}z_{i}\right)^{2}}\;.$

We can use the Cauchy-Schwarz inequality to refine this condition further
more. Indeed, put $\sigma=\min\left\\{\sigma_{1},\sigma_{2}\right\\}$. Then by
the Cauchy Schwarz inequality,

$\left(\sum_{i=1}^{2}c_{i}z_{i}\right)^{2}\leq\sum_{i=1}^{2}1^{2}\sum_{i=1}^{2}c_{i}^{2}z_{i}^{2}=2\sum_{i=1}^{2}c_{i}^{2}z_{i}^{2}\;.$

Then
$\displaystyle\sum_{i=1}^{2}\sigma_{i}^{2}c_{i}^{2}z_{i}^{2}\geq\frac{1}{2}\sigma^{2}\left(\sum_{i=1}^{2}c_{i}z_{i}\right)^{2}$
. Hence $A_{3}$ will be satisfied if we can find $\gamma_{b}>0$ and ${\bf
c}=(c_{1},c_{2})\in\mathbb{R}^{2,\circ}_{+}$ such that

$k+\gamma_{b}(3+k)\leq\frac{\sigma}{4}\frac{\left(\sum_{i=1}^{2}c_{i}z_{i}\right)^{2}}{\left(1+\sum_{i=1}^{2}c_{i}z_{i}\right)^{2}}\;.$

Now let $z_{1},z_{2}\in\mathbb{R}_{+}$. We pick a $\gamma_{0}\in(0,1)$. There
exists $x\in\mathbb{R}^{+}$ such that $x\geq\frac{\gamma_{0}}{1-\gamma_{0}}$.
Therefore, we will have $\frac{x}{1+x}\geq\gamma_{0}$ which implies
$\frac{x^{2}}{(1+x)^{2}}\geq\gamma_{0}^{2}$. Choosing in particular $x$ in the
interval
$[\min\left\\{z_{2},z_{2}\right\\},\max\left\\{z_{1},z_{2}\right\\}]$, there
exists $\bm{c}=(c_{1},c_{2})\in\mathbb{R}^{2,\circ}_{+}$ such that
$x=\sum_{i=1}^{2}c_{i}z_{i}$ and
$\displaystyle\frac{\left(\sum_{i=1}^{2}c_{i}z_{i}\right)^{2}}{(1+\sum_{i=1}^{2}c_{i}z_{i})^{2}}\geq\gamma_{0}^{2}$.
To finish, we choose $\gamma_{b}$ such that
$\displaystyle\gamma_{0}^{2}=\frac{4(k+\gamma_{b}(3+k))}{\sigma}$ and the
proof is complete.