# Coexistence for a kind of stochastic three-species competitive models

• Nantian Huang , Jiabing Huang , Yuming Wei and Yongjian Liu
From the journal Open Mathematics

## Abstract

The coexistence of species sustains the ecological balance in nature. This paper focuses on sufficient conditions for the coexistence of a three-species stochastic competitive model, where the model has non-linear diffusion parts. Three values λ3z, λ3x and λ3y are introduced and calculated from the coefficients, which can be considered as threshold values. Moreover, convergence in distribution of the positive solution of the model is also addressed. A few numerical simulations are carried out to illustrate the theoretical results.

MSC 2010: 60H10; 60J70; 92D25; 93D15

## 1 Introduction

The coexistence of species plays a vital role in protecting ecological balance, and the situation is not desirable when any species extinct. Therefore, this paper aims to discuss what the conditions for coexistence is. In term of deterministic model, Lotka-Volterra model, which was applied to study species interactions, was originally proposed by A. Lotka and V. Volterra [1, 2, 3]. Later Lotka-Volterra model is widely applied in the field of chemistry. Literature [4], based on the reversible Lotka-Volterra model, investigated the dynamic behavior of oscillatory reaction. In economy, Lotka-Volterra model is also a very useful tool to analyse enterprise competition [5]. Nowadays, Lotka-Volterra model plays a vital role in engineering field.

The classical three-dimensional competitive Lotka-Volterra model can be written as

dx(t)dt=x(t)(r1a11x(t)a12y(t)a13z(t)),dy(t)dt=y(t)(r2a21x(t)a22y(t)a23z(t)),dz(t)dt=z(t)(r3a31x(t)a32y(t)a33z(t)), (1.1)

where x(t), y(t) and z(t) denote the densities of three species at time t. For i, j = 1, 2, 3, the parameters ri, aij are all positive; ri stand for intrinsic growth rates; aii describe as intra-specific competition rates; aij(ij) represent the inter-specific competition rates. For deterministic Lotka-Volterra model, quite a few scholars contribute to investigate coexistence of species (see [6, 7, 8, 9, 10, 11, 12] and references therein).

However, as the development of stochastic analysis, deterministic model may not be suitable for reality sometimes. As a matter of fact, population systems are usually subject to environmental noise. Therefore, to describe more realistically population systems, stochastic population models have attracted considerable attention [13, 14, 15, 16, 17]. For the sake of generality, suppose that all coefficients of model(1.1) are perturbed by Brownian motions, and the model becomes

dx(t)=x(t)(r1a11x(t)a12y(t)a13z(t))dt+(Δ1x(t)+σ11x2(t))dB1(t)+σ12x(t)y(t)dB2(t)+σ13x(t)z(t)dB3(t),dy(t)=y(t)(r2a21x(t)a22y(t)a23z(t))dt+σ21x(t)y(t)dB2(t)+(Δ2y(t)+σ22y2(t))dB4(t)+σ23y(t)z(t)dB5(t),dz(t)=z(t)(r3a31x(t)a32y(t)a33z(t))dt+σ31x(t)z(t)dB3(t)+σ22y(t)z(t)dB5(t)+(Δ3z(t)+σ33z2(t))dB6(t), (1.2)

where Bi (i = 1, 2, 3, 4, 5, 6) are mutually independent Brownian motions. Although there may be an appropriate Lyapunov function to provide the condition for coexistence, it is really hard to be found in practice. Motivated this, some skillful technique should be introduced to solve this problem. To reduce unnecessary computations due to notational complexity and to make our ideas more understandable but still preserve important properties, we assume that the lowest-power terms are not affected by environment noise for simplicity, which means Δ1 = Δ2 = Δ3 = 0. Thus, throughout the rest of the paper, the following model will be considered

dx(t)=x(t)(r1a11x(t)a12y(t)a13z(t))dt+σ11x2(t)dB1(t)+σ12x(t)y(t)dB2(t)+σ13x(t)z(t)dB3(t),dy(t)=y(t)(r2a21x(t)a22y(t)a23z(t))dt+σ21x(t)y(t)dB2(t)+σ22y2(t)dB4(t)+σ23y(t)z(t)dB5(t),dz(t)=z(t)(r3a31x(t)a32y(t)a33z(t))dt+σ31x(t)z(t)dB3(t)+σ22y(t)z(t)dB5(t)+σ33z2(t)dB6(t). (1.3)

To proceed, The rest of the paper is as follows. In Section 2, basic concept and three essential constants λ3z, λ3x and λ3y are provided. Sections 3 are devoted to proving the stochastic coexistence. Then a few numerical simulations are given to manifest our results in Section 4. The last but not least, brief discussions are made to conclude the paper.

## 2 Basic concept and essential constants

Let (Ω, 𝓕, {𝓕t}t≥0, ℙ) be a complete filtered probability space with the filtration {𝓕t}t≥0 satisfying the usual condition, i.e., it is increasing and right continuous while 𝓕0 contains all ℙ-null sets. In model (1.3), Bi (i = 1, 2, 3, 4, 5, 6) are 𝓕t-adapted, mutually independent Brownian motions. Suppose that aij (i, j = 1, 2, 3) are positive constants. Set σii ≠ 0 (i = 1, 2, 3) in order to make the diffusion be non-degenerate. For convenience, denote w = (x, y, z), w0 = (x0, y0, z0), and w(t) = (x(t), y(t), z(t)). Denote ab = max{a, b}, ab = min{a, b}, R+2,o = {(x, y) : x > 0, y > 0} and R+3,o = {(x, y, z) : x > 0, y > 0, z > 0}. Let ww0(t) = (xw0(t), yw0(t), zw0(t)) be the solution to (1.3) with initial value w0. Previous study [15] proved that ww0(t) remains in R+3,o with probability one if w0 R+3,o . The definition of coexistence is given as following.

## Definition 3.1

For model (1.3), stochastic coexistence takes place in three species if for any ε > 0, there exists an M = M(ε) > 1 satisfying that

lim inftPM1x(t),y(t),z(t)M1ε.

Next, a few useful properties of solutions are presented, whose proofs can be found in [15, 19, 20].

## Proposition 2.1

The following assertions hold:

1. For all z0 R+3 \ {(0, 0, 0)}, there exists an M0 > 0 satisfying that

lim suptEV(xw0(t),yw0(t),zw0(t))M0,

where V(x, y, z) = (x + y + z) + (x + y + z)−1.

2. For any ε > 0, H > 1, T > 0, there exists an H = H(ε, H, T) > 1 satisfying that

P{H¯1xw0(t)H¯t[0,T]}1εifw0[H1,H]×[0,H]×[0,H],P{H¯1yw0(t)H¯t[0,T]}1εifw0[0,H]×[H1,H]×[0,H]

and that

P{H¯1zw0(t)H¯t[0,T]}1εifw0[0,H]×[0,H]×[H1,H].
3. For any p ∈ (0, 3), there exists an Mp > 0 satisfying that

E0t|ww0(s)|p dsMp(t+|w0|)z0R+2,t0,

where |(x,y,z)|=x2+y2+z2.

To further study, the equations on the boundary need to be considered. On the x-axis, one has

dφ(t)=φ(t)(r1a11φ(t))dt+σ11φ2(t)dB1(t). (2.1)

There exists a unique stationary distribution π1 in (0, ∞) with density

f1(u)=c1u4exp2a11σ1121ur1σ1121u2,u>0,

where c1=0x2exp2a11σ112xr1σ112x2dx1. Owing to the ergodicity, for any measurable function h( ⋅) : ℝ+ → ℝ satisfying that 0|h(u)|f1(u)du<, the following equation holds,

PlimT1T0Th(φx0(t))dt=0h(u)f1(u)du=1,x0>0, (2.2)

where φx0(t) is the solution to (2.1) starting at x0. Specially, for any p ∈ (−∞, 3),

PlimT1T0Tφx0p(t)dt=Qpx:=0upf1(u)du<=1,x0>0. (2.3)

Similarly, on the y-axis, one has

dψ(t)=ψ(t)(r2a22ψ(t))dt+σ22ψ2(t)dB4(t), (2.4)

which in (0, ∞) has a unique stationary distribution π2 with density

f2(v)=c2v4exp2a22σ2221vr2σ2221v2,v>0,c2=0x2exp2a 22σ222xr2σ222x2dx1.

By the ergodicity, for p ∈ (−∞, 3), one gets

PlimT1T0Tψy0p(t)dt=Qpy:=0vpf2(v)dv<=1,x0>0, (2.5)

where ψy0(t) is the solution to (2.4) starting at y0.

On the z-axis, one has

dρ(t)=ρ(t)(r3a33ρ(t))dt+σ33ρ2(t)dB6(t), (2.6)

which in (0, ∞) has a unique stationary distribution π3 with density

f3(ϕ)=c3ϕ4exp2a33σ3321ϕr3σ3321ϕ2,ϕ>0c3=0x2exp2a33σ332xr3σ332x2dx1.

In view of the ergodicity, for p ∈ (−∞, 3), one has

PlimT1T0Tρz0p(t)dt=Qpz:=0ϕpf3(ϕ)dϕ<=1,z0>0, (2.7)

where ρz0(t) is the solution to (2.6) starting at z0. As to the xy plane, Literature [18] has considered following model

dX(t)=X(t)(r1a11X(t)a12Y(t))dt+σ11X2(t)dB1(t)+σ12X(t)Y(t)dB2(t),dY(t)=Y(t)(r2a21X(t)a22Y(t))dt+σ21X(t)Y(t)dB2(t)+σ22Y2(t)dB4(t), (2.8)

which can be seen as a case on the x - y plane for (1.3). Its result provides that (2.8) has a unique invariant measure μ1 in R+2,o if

λ1z:=0(r2a21xσ2122x2)f1(x)dx=r2a21Q1xσ2122Q2x>0

and

λ2z:=0(r1a12yσ1222y2)f2(y)dy=r1a12Q1yσ1222Q2y>0.

Moreover, for any μ1 -integrable function F(x, y) : R+2,o → ℝ, one has

limt1t0tF(Xγ(s),Yγ(s))ds=R+2,oF(x,y)μ1(dx,dy)a.s.γR+2,o. (2.9)

Let f1(x, y) be the density of probability measure μ1 in R+2,o . Then f1(x, y) has following properties,

00f1(x,y)dxdy=1,0f1(x,y)dy=f1(x)and0f1(x,y)dx=f2(y),

where f1 (x) and f2 (y) are density of π1 and π2 respectively. Define that

λ3z:=r30(a31x+σ3122x2)f1(x)dx0(a32y+σ3222y2)f2(y)dy=r3a31Q1xσ3122Q2xa32Q1yσ3222Q2y. (2.10)

Similarly, on y-z plane, one gets

dY(t)=Y(t)(r2a22Y(t)a23Z(t))dt+σ22Y2(t)dB4(t)+σ23Y(t)Z(t)dB5(t),dZ(t)=Z(t)(r3a32Y(t)a33Z(t))dt+σ32Y(t)Z(t)dB5(t)+σ33Z2(t)dB6(t). (2.11)

Using the method of [18], if

λ1x:=0(r3a32yσ3222y2)f2(y)dy=r3a32Q1yσ3222Q2y>0

and

λ2x:=0(r2a23zσ2322z2)f3(z)dϕ=r2a23Q1zσ2322Q2z>0,

then (2.11) has a unique invariant measure μ2 in R+2,o Moreover, for any μ2 -integrable function F(y, z): R+2,o → ℝ, one has

limt1t0tF(Yγ(s),Zγ(s))ds=R+2,oF(y,z)μ2(dy,dz)a.s.γR+2,o.

Let f2(y, z) be the density of probability measure μ1 in R+2,o , and then f2(y, z) has following properties,

00f2(y,z)dydz=1,0f2(y,z)dz=f2(y)and0f2(y,z)dy=f3(z),

where f2 (y) and f3 (z) are density of π2 and π3 respectively. Define that

λ3x:=r10(a12y+σ1222y2)f2(y)dy0(a13z+σ1322z2)f3(z)dz=r1a12Q1yσ1222Q2ya13Q1zσ1322Q2z. (2.12)

Analogously, on x - z plane, one has

dX(t)=X(t)(r1a11X(t)a13Z(t))dt+σ11X2(t)dB1(t)+σ13X(t)Z(t)dB3(t),dZ(t)=Z(t)(r3a31X(t)a33Z(t))dt+σ31X(t)Z(t)dB3(t)+σ33Z2(t)dB6(t). (2.13)

If

λ1y:=0(r3a31xσ3122x2)f2(x)dx=r3a31Q1xσ3122Q2x>0

and

λ2y:=0(r1a13zσ1322z2)f3(z)dz=r1a13Q1zσ1322Q2z>0,

then (2.13) has a unique invariant measure μ3 in R+2,o Moreover, for any μ3 -integrable function F(x, z) : R+2,o → ℝ, one has

limt1t0tF(Xγ(s),Zγ(s))ds=R+2,oF(x,z)μ3(dx,dz)a.s.γR+2,o.

Let f3(x, z) be the density of probability measure μ3 in R+2,o , and then f3(x, z) has following properties,

00f3(x,z)dxdz=1,0f3(x,z)dz=f1(x)and0f3(x,z)dy=f3(z),

where f1 (x) and f3 (z) are density of π1 and π3 respectively. Define that

λ3y:=r20(a21x+σ2122x2)f1(x)dx0(a23z+σ2322z2)f3(z)dz=r2a21Q1xσ2122Q2xa23Q1zσ2322Q2z. (2.14)

The reason of definition of λ3z, λ3x and λ3y is as following. To determine whether zw0(t) converges to 0 or not, it is required to consider the Lyapunov exponent of zw0(t). The Itô’s formula manifests that

lnzw0(T)T=lnz0T+1T0Tr3a31xw0(t)σ3122xw02(t)a32yw0(t)σ3222yw02(t)a33zw0(t)σ3322zw02(t)dt+1T0T(σ31xw0(t)dB3(t)+σ32yw0(t)dB5(t)+σ33zw0(t)dB6(t)). (2.15)

When T is large enough, the first and the third terms on the right-side of (2.15) are small. Clearly, if zw0(t) is small in [0, T], the solution xw0(t) is close to X(x0,y0)(t) and yw0(t) is close to Y(x0,y0)(t), where X(x0,y0)(t) and Y(x0,y0)(t) are solutions to (2.8) starting at (x0, y0). Employing the ergodicity (2.9), lnzw0(T)T is close to λ3z. The definitions of λ3x and λ3y are also gotten in the similar way.

## 3 Main result

This section focuses on providing the conditions for stochastic coexistence in model (1.3).

## Theorem 3.1

If λ3z, λ3x and λ3y are all positive, the three species coexist. Moreover, in R+3,o , the solution process w(t) has a unique invariant measure m* satisfying that

1. the transition probability ℙ(t, w0, ⋅) of z(t) converges to m*, ∀ w0 R+3,o ;

2. for any m*-integrable function J(w) : R+3,o → ℝ, one has

limt1t0tJ(ww0(s))ds=R+3,oJ(w)m(dw)a.s.w0R+3,o.

To prove Theorem 3.1, some arguments need to be introduced.

The following formula is the well-known exponential martingale inequality. It asserts that for any a, b > 0,

P0tg(s)dB(s)a20tg2(s)ds>bt0eab, (3.1)

if B(s) is a 𝓕t-adapted Brownian motion while g(t) is a real-valued 𝓕t-adapted process and 0tg2(s)ds < ∞, ∀ t ≥ 0 almost surely. Ordinarily, the inequality holds in finite time, but (3.1) holds since

P0tg(s)dB(s)a20tg2(s)ds>bt0=limTP0tg(s)dB(s)a20tg2(s)ds>bt[0,T].

Let any T > 1, p ∈ (1, 1.5), and 1p+1q = 1. For A ∈ 𝓕, IA stands for the indicator function of A. Owing to part (iii) of Proposition 2.1 and Holder’s inequality, it is estimated that

EIA|lnzw0lnz0|E0TIAr3+a31xw0(t)+σ3122xw02(t)+a32yw0(t)+σ3222yw02(t)+a33zw0(t)+σ3322zw02(t)dt+EIA0T(σ31xw0(t)dB3(t)+σ32yw0(t)dB5(t))+σ33zw0(t)dB6(t))θ1P(A)T1qE0T(1+xw02p(t)+yw02p(t)+zw02p(t))dt1p+P(A)(|σ31||σ32||σ33|)E0Txw02(t)+yw02(t)+zw02(t)dtθ2T[P(A)]1q(1+|w0|)1p, (3.2)

for some constants θ1, θ2 independent of z0, T and A. In particular, when A = Ω,

Elnzw0(T)lnz0Tθ2(1+|w0|)1p, (3.3)

and consequently,

Plnyz0(T)lny0Tθ3(1+|w0|)1pεε. (3.4)

To satisfy the condition that zw0(t) is small in [0, T], it requires to introduce Lemma 3.1. Let τw0σ be the stopping time that τw0σ = inf{t ≥ 0 : zw0(t) ≥ σ}.

## Lemma 3.1

For any ε > 0, σ > 0, T > 1, there is a δ = δ(T, ε, σ) > 0 satisfying that

P{τw0σT}1ε,w0(0,)×(0,)×(0,δ].

## Proof

By (3.1), ℙ( Ω1w0 ) ≥ 1 − ε, where

Ω1w0=0t[σ31xw0(s)dB3(s)+σ32yw0(s)dB5(s)+σ33zw0(s)dB6(s)]120t[σ312xw02(s)+σ322yw02(s)+σ332zw02(s)]ds<ln1ε,t0.

In view of (2.15), when ω Ω1w0 , it is seen that

lnzw0(t)<lnz0+ln1ε+0tr3ds0t[a31xw0(s)+a32yw0(s)+a33zw0(s)]ds<lnz0+ln1ε+r3t.

Letting δ = σ ε er3T, if z0δ, then zw0(t) < σ for all t < T and ω Ω1w0 . The proof is complete.

## Lemma 3.2

For any ε, β > 0 and H, T > 1, there exists σ > 0 satisfying that for all w0 ∈ [H−1, H] × [H−1, H] × [0, σ],

P{|X(x0,y0)(t)xw0(t)|<β and|Y(x0,y0)(t)yw0(t)|<β,0tTτz0σ}1ε.

## Proof

In view of part (ii) of Proposition 2.1, there is an H sufficiently large such that

Pmax{X(x0,y0)(t),xw0(t),Y(x0,y0)(t),yw0(t)}H¯,tT>1ε,w0[H1,H]×[H1,H]×(0,σ].

Let ξw0 := τw0σ ∧ inf {t : max {X(x0,y0)(t), xw0(t), Y(x0,y0)(t), yw0(t)} ≤ H}. The Itô’s formula manifests that

|X(x0,y0)(s)xw0(s)|0s|X(x0,y0)(u)xw0(u)|r1+a11(X(x0,y0)(u)+xw0(u))du+a120sX(x0,y0)(u)Y(x0,y0)(u)xw0(s)yw0(u)du+a130s|xw0(u)zw0(u)|du+|σ11|0sX(x0,y0)(u)xw0(u)X(x0,y0)(u)+xw0(u)dB1(u)+|σ12|0sX(x0,y0)(u)Y(x0,y0)(u)xw0(u)yw0(u)dB2(u)+|σ13|0sxw0(u)zw0(u)dB3(u).

By the elementary inequality (Σi=1nai)22nΣi=1nai2 leads to

Esupst(X(x0,y0)(sξw0)xw0(sξw0))264E0tξw0|X(x0,y0)(u)xw0(u)|r1+a11(X(x0,y0)(u)+xw0(u))du2+64a122E0tξw0[X(x0,y0)(u)Y(x0,y0)(u)xw0(s)yw0(u)]2du+64a13E0rξw0xw02(u)zw02(u)du+64σ112Esupst0sξw0X(x0,y0)(u)xw0(u)X(x0,y0)(u)+xw0(u)dB1(u)2+64σ122Esupst0sξw0X(x0,y0)(u)Y(x0,y0)(u)xw0(u)yw0(u)dB2(u)2+64σ132Esupst0sξw0xw0(u)zw0(u)dB3(u)2. (3.5)

For t ∈ [0, T], a few estimates are made in the following,

E0tξw0|X(x0,y0)(u)xw0(u)|r1+a11(X(x0,y0)(u)+xw0(u))du2
(r1+2a11H¯)2TE0tξw0X(x0,y0)(u)xw0(u)2ds, (3.6)
E0tξw0[X(x0,y0)(u)Y(x0,y0)(u)xw0(s)yw0(u)]2du4H¯4T, (3.7)
E0rξw0xw02(u)zw02(u)duH¯2σ2T. (3.8)

It follows from the Burkholder-Davis-Gundy inequality that

Esupst0sξw0X(x0,y0)(u)xw0(u)X(x0,y0)(u)+xw0(u)dB1(u)2E0tξw0X(x0,y0)(u)xw0(u)2X(x0,y0)(u)+xw0(u)2du,4H¯2E0tξw0X(x0,y0)(u)xw0(u)2du, (3.9)
Esupst0sξw0X(x0,y0)(u)Y(x0,y0)(u)xw0(u)yw0(u)dB2(u)2E0tξw0X(x0,y0)(u)Y(x0,y0)(u)xw0(u)yw0(u)2du,4H¯4T, (3.10)
Esupst0sξw0xw0(u)zw0(u)dB3(u)24H¯2σ2T. (3.11)

Applying (3.6), (3.7), (3.8), (3.9), (3.10) and (3.11) to (3.5), one has

Esupst(X(x0,y0)(sξw0)xw0(sξw0))2θ¯σ2+E0tξw0(X(x0,y0)(u)xw0(u))2duθ¯σ2+0tEsupsu(X(x0,y0)(sξw0)xw0(sξw0))2dut[0,T],

where θ = θ(H, T) > 0. It then follows from Granwall’s inequality that

EsupsT(X(x0,y0)(sξw0)xw0(sξw0))2θ¯σ2eθ¯T.

Consequently, PsupsT(X(x0,y0)(sξw0)xw0(sξw0))2β2θ¯σ2eθ¯Tβ2<ε4 when σ is sufficiently small σ2<εβ24θ¯eθ¯T. The similar way yields that

PsupsT(Y(x0,y0)(sξw0)yw0(sξw0))2β2<ε4.

Since

Psξw0=sτw0σ,s[0,T]PsupsT{maxxw0(s),yw0(s),X(x0,y0)(s),Y(x0,y0)(s)H¯1ε2,

it follows that

P|X(x0,y0)(s)xw0(s)|<βand|Y(x0,y0)(s)yw0(s)|<β,0tTτw0σPsξw0=sτw0σs[0,T]Psuprt(X(x0,y0)(rξw0)xw0(rξw0))2β2Psuprt(X(x0,y0)(rξw0)xw0(rξw0))2β21ε.

The proof is complete.

## Remark 3.1

Lemma 3.2 aims at describing that xw0(t) is close to X(x0,y0)(t) and yw0(t) is close to Y(x0,y0)(t) in certain conditions.

## Lemma 3.3

For any ε > 0, there exists an M > 0 satisfying that

P0Tσ31xw0(t)dB3(t)+σ32yw0(t)dB5(t)+σ33zw0(t)dB6(t)M¯εT|w0|1ε.

## Proof

Since

E0Tσ31xw0(t)dB3(t)+σ32yw0(t)dB5(t)+σ33zw0(t)dB6(t)2=E0Tσ312xz02(t)+σ322yw02(t)+σ332zw02(t)dt.

By (iii) of Proposition 2.1 and Chebyshev inequality, the desired result can be obtained.

Two propositions should be introduced to make proof of Theorem 3.1 clear.

## Proposition 3.1

Suppose λ1z λ2z and λ3z > 0. For any ε > 0, H > 1, there exist T = T(ε, H) and δ0 = δ0(ε, H) satisfying that ℙ(Ωw0) > 1 − 4ε, where

Ω¯w0=λ3z5lnzw0(t)lnz0,w0[H1,H]×[H1,H]×[0,δ0].

## Proof

In view of (2.10), for sufficiently small β, one has

R+2,or3a31(u+β)σ3122(u+β)2a32(v+β)σ3222(v+β)2f1(u,v)dudv4λ3z5.

Let M be as in Lemma 3.3. Using (2.9), there is T=T(ε,T)>25HM¯2ε2λ3z2 such that

P(ΩH)1ε,

where

ΩH=1T0r3a31(X(H,H)(t)+β)σ3122(X(H,H)(t)+β)2a31(Y(H,H)(t)+β)σ3122(Y(H,H)(t)+β)2dt3λ3z5.

In view of the uniqueness of solution, for all (x0, y0) ∈ [H−1, H] × [H−1, H], one has X(x0,y0)(t) ≤ X(H,H)(t) and Y(x0,y0)(t) ≤ Y(H,H)(t). Therefore, ℙ( Ω2w0 ) ≥ 1 − ε where

Ω2w0=1T0Tr3a31(X(x0,y0)(t)+β)σ3122(X(x0,y0)(t)+β)2a31(Y(x0,y0)(t)+β)σ3122(Y(x0,y0)(t)+β)2dt3λ3z5.

Owing to Lemma 3.2, there is a σ = σ(ε, H) > 0 satisfying that

a33σ+σ3322σ2<λ3z5

and ℙ( Ω3w0 ) ≥ 1 − ε, where

Ω3w0=|X(x0,y0)(t)xw0(t)|αand|Y(x0,y0)(t)yw0(t)|α,0tTτw0σ.

By virtue of Lemma 3.1, there exists δ0 = δ0(ε, H) satisfying that ℙ( Ω4w0 ) ≥ 1 − ε for all w0 ∈ [H−1, H] ×[H−1, H] × (0, δ0], where

Ω4w0={τw0σT}.

Since T>25M¯2Hε2λ3z2, it is derived from Lemma 3.3 that ℙ( Ω5w0 ) ≥ 1 − ε, where

Ω5w0=P0Tσ31xw0(t)dB3(t)+σ32yw0(t)dB5(t)+σ33zw0(t)dB6(t)λ3z5T.

Thus, for w0 ∈ [H−1, H]×[H−1, H] × (0, δ0] and ωi=25Ωiw0, it follows that

lnzw0(T)lnz00Tr3a31xw0(t)σ3122xw0(t)a32yw0(t)σ3222yw0(t)dt0Ta33zw0(t)σ3322zw0(t)dt0Tσ31xw0(t)dB3(t)+σ32yw0(t)dB5(t)+σ33zw0(t)dB6(t)0Ta2c2(uz0(t)+α)β222(uz0(t)+α)2dt3λ3z5Tλ3z5Tλ3z5Tλ3z5T

and P{λ3z5Tlnzw0(T)lnz0}Pi=25Ωiw014ε. The proof is complete.

## Proposition 3.2

If λ1z, λ2z and λ3z > 0, then for any Δ > 0, there exist T = T(Δ) and δ2 = δ2(Δ) satisfying that

lim supn1nk=0n1Pzw0(kT)δ2Δ,w0R+3,o.

## Proof

Let ε = ε(Δ) ∈ (0, 1) and H = H(Δ) > 1 be chosen later. Put Λ=θ2ε(1+3H)1p where θ2 is as in (3.2). In view of (3.4) and Proposition 3.1, there exist δ0 ∈ (0, 1) and T > 1 such that for all w0 ∈ [H−1, H] × [H−1, H] × (0, δ0], P(Ω0z0) > 1 − 5ε where

Ω0z0=λ3z5Tlnzw0(T)lnz0ΛT.

Let L1 = Λ T + ln H − ln δ0. By using | ln H − ln zw0(T)| ≤ |ln H − ln z0| + |ln zw0(T) − ln z0|, if w0 ∈ [H−1, H]× [H−1, H] × (0, δ0], one can be derived from (3.4) that

P{|lnHlnzw0(T)|L1}P{|lnzw0(T)lnz0|ΛT}5ε. (3.12)

Let δ1, δ2 satisfy L1 = ln H − ln δ1, L2 := L1 + Λ T = ln H − ln δ2. It’s clear that δ2 < δ1 < δ0. Define U(z) = (ln H − ln z) ∨ L1. Obviously,

U(z)U(z)|lnzlnz|. (3.13)

Now, 1T[EU(zw0(T)U(z0))] needs to be estimated for different w0. First, it follows form (3.2) and (3.13) that for any w0 R+3,o and ωΩ0w0,c=ΩΩ0w0,

1TEIΩ0w0,cU(zw0(T))U(z0)1TEIΩ0w0,c|lnzw0(T)lnz0|θ2(1+3H)1p[P(Ω0w0,c)]1qθ2(1+3H)1p(5ε)1q. (3.14)

If w0D1 := [H−1, H] × [H−1, H] × (0, δ2], then ln H − ln z0 ≥ ln H − ln δ2L1 + Λ T. It follows that in Ω0ω0 ,

lnHlnzw0(T)lnHlnz0ΛTL2ΛTL1,lnHlnzw0(T)lnHlnz0λ3z5T.

Consequently, one gets

1TEIΩ0w0U(zw0(T))U(z0)=1TP(Ω0w0)lnHU(zw0(T))lnH+U(z0)λ3z5(1ε). (3.15)

One can be derived from (3.14) and (3.15) that for all w0D1,

1TEU(zw0(T))U(z0)1TEIΩ0w0+IΩ0w0,cU(zw0(T))U(z0)λ3z5(1+ε)+θ2(1+3H)1p(5ε)1q. (3.16)

If w0D2 := [H−1, H] ×[H−1, H] × (δ2, δ1], then ln H − ln z0 ≥ ln H − ln z0 = L1. In Ω0w0 ,

lnHlnzw0(T)lnHlnz0λ3z5TlnHlnz0=U(z0).

Therefore, U(yw0(T)) − U(w0) ≤ 0 in Ω0w0 . As a result of (3.14),

1TEU(zw0(T))U(z0)θ2(1+3H)1p(5ε)1q,w0D2. (3.17)

If w0D3 := [H−1, H] × [H−1, H] × (δ1, δ0], then ln H – ln z0 ≤ ln H – ln δ1 = L1. And in Ω0w0 , ln H – ln zw0(T) ≤ ln H – ln z0L1, which means that U(zw0(T)) = U(z0) = L1. Therefore,

1TEU(zw0(T))U(z0)θ2(1+3H)1p(5ε)1q,w0D3. (3.18)

If w0D4 := [H−1, H] × [H−1, H] × (δ0, H], it leads to that ln H – ln z0 ≤ ln H – ln δ0 = L1ΛT. In the same way of obtaining (3.18), one gets

1TEU(zw0(T))U(z0)θ2(1+3H)1p(5ε)1q,w0D4. (3.19)

For any initial value w0 R+3,o , the Markov property of w(t) leads to

1TE[U(zw0(kT+T))U(zw0(kT))]R+2,oPzw0(kT)dw1TEU(zw(T))U(z0).

Putting D := [H−1, H] × [H−1, H] × [H−1, H] = i=14 Di, and using (3.2), (3.16), (3.17), (3.18) and (3.19), one gets

1TE[U(zw0(kT+T))U(zw0(kT))]Pww0(kT)D1λ3z5(1ε+ε1)+ε1Pww0(kT)DD1+EI{ww0(kT)D}1TEln(zww0(kT)(T))ln(zw0(kT))λ3z5(1ε)Pww0(kT)D1+ε1+θ21+E|ww0(kT)|1pPww0(kT)D1q, (3.20)

where θ2 is as in (3.2) and ε1=θ2(1+3H)1p(5ε)1q. It follows from Markov’s inequality and (i) of Proposition 2.1 that

lim supkP{ww0(kT)D}lim supkP{[xw0(kT)yw0(kT)zw0(kT)]>H}lim supkP{V(ww0(kT))>H}lim supkEV(ww0(kT))HM0H (3.21)

and

lim supk1+E|ww0(kT)|1pPww0(kT)D1qM0H1qlim supk[1+3EV(ww0(kT)]1p1+3M0H1q1+3M01p1+3M0H1q. (3.22)

Obviously,

lim infn1nTk=0n1EU(zw0(kT+T))U(zw0(kT))=lim infn1nTEU(zw0(nT)U(z0))0. (3.23)

It follows from (3.20), (3.21) and (3.23) that

0lim supn1nk=0n1Pww0(kT)D1λ3z5(1ε)+ε1+θ2(1+3M0)H1q
lim supn1nk=0n1Pww0(kT)D15λ3z(1ε)ε1+θ2(1+3M0)H1q, (3.24)

where ε1=θ2(1+3H)1p(5ε)1q. It should be noted that ℙ{zw0(kT) ≤ δ2} ≤ ℙ{ww0(kT) ∈ D1} + ℙ{ww0(kT) ∉ D}. As a result of (3.21) and (3.24), one has

lim supn1nk=0n1P{zw0(kT)δ2}5λ3z(1ε)ε1+θ2(1+3M0)H1q+M0H.

Consequently, there are sufficiently large H = H(Δ) and sufficiently small ε = ε(Δ) satisfying the desired result. The proof is complete.

Now, it is time to prove Theorem 3.1.

## Proof of Theorem 3.1

Set any ε > 0. Owing to symmetry of (1.3), similar to Proposition 3.2, it can be proved that if λ1x, λ2x and λ3x > 0, then there are T′ > 1 and δ2 > 0 satisfying that

lim supn1nk=0n1Pxw0(kT)δ2Δ,w0R+3,o.

And if λ1y, λ2y and λ3y > 0, then there are T″ > 1 and δ2 > 0 satisfying that

lim supn1nk=0n1Pyw0(kT)δ2Δ,w0R+3,o.

In fact, if λ3z, λ3x and λ3y > 0, then λ1z, λ2z, λ1x, λ2x, λ1y and λ2y will be positive, and λ3z, λ3x and λ3y have been defined in (2.10), (2.12) and (2.14) respectively. Furthermore, by choosing sufficiently large T and sufficiently small δ2, one has

lim supn1nk=0n1P|xw0(kT)||yw0(kT)||zw0(kT)|δ23Δ,w0R+3,o.

Combining this with (i) of Proposition 2.1 implies that there exists a compact set G R+3,o satisfying that

lim inf1nk=0n1P{ww0(kT)G}14Δ,w0R+3,o.

Owing to (ii) of Proposition 2.1, there is an M > 1 satisfying that

P{M1xw(t),yw(t),zw(t)M}1Δ,tT,wG.

By virtue of Markov property, one has

P{M1xw0(kT+t),yw0(kT+t),zw0(kT+t)M}P{ww0(kT)G}P{M1xw(t),yw(t),zw(t)M},wG,w0R+3,o

As a result, for any w0 R+3,o ,

lim infn1nT0nTP{M1xw0(t),yw0(t),zw0(t)M}dt(1Δ)lim infn1nk=0n1P{ww0(kT)G}(1Δ)(14Δ)15Δ.

It means that

lim inft1t0tP{M1xw0(s),yw0(s),zw0(s)M}ds15Δ,

and there exist an invariant probability measure consequently. Since the diffusion is non-degeneracy, the rest of the results are yield (see [21]). The proof is complete.

## 4 Examples and complexity discussion

In this section, we consider model (1.3),

dx(t)=x(t)(r1a11x(t)a12y(t)a13z(t))dt+σ11x2(t)dB1(t)+σ12x(t)y(t)dB2(t)+σ13x(t)z(t)dB3(t),dy(t)=y(t)(r2a21x(t)a22y(t)a23z(t))dt+σ21x(t)y(t)dB2(t)+σ22y2(t)dB4(t)+σ23y(t)z(t)dB5(t),dz(t)=z(t)(r3a31x(t)a32y(t)a33z(t))dt+σ31x(t)z(t)dB3(t)+σ22y(t)z(t)dB5(t)+σ33z2(t)dB6(t).

According to Theorems 3.1, λ3z, λ3x and λ3y govern the coexistence of three species; λ3z, λ3x and λ3y > 0 lead to coexistence of three species; Note that λ3z, λ3x and λ3y are defined in (2.10), (2.12) and (2.14) respectively. Examples 4.1 and the simulations illustrate theoretical results that three species coexist.

## Example 4.1

Consider (1.3) with parameters r1 = 6; r2 = 6; r3 = 6; a11 = 1; a12 = 1; a13 = 2; a21 = 2; a22 = 2; a23 = 1; a31 = 1; a32 = 3; a33 = 3; σ11 = 2; σ12 = 1; σ13 = 0.5; σ21 = 1; σ22 = 2; σ23 = 1; σ31 = 0.5; σ32 = 1; σ33 = 2. Direct calculation demonstrates that λ3z = 0.1258, λ3x = 1.7079, λ3y = 0.5557. As a result of Theorem 3.1, the three species coexist. Sample paths of xw0(t), yw0(t) and zw0(t) with w0 = (3, 3, 3) are illustrated in Fig. 1.

Figure 1

Sample paths of xw0(t), yw0(t) and zw0(t) in Example 5.1 with w0 = (3, 3, 3).

## Remark 4.1

When λ3z, λ3x and λ3y are not all positive, the properties of the solution will become very complicated. Three species may die out or coexist; see Example 4.2.

## Example 4.2

Consider (1.3) with parameters r1 = 6; r2 = 5; r3 = 4; a11 = 2; a12 = 1; a13 = 1; a21 = 1; a22 = 2; a23 = 1; a31 = 1; a32 = 1; a33 = 2; σ11 = 1; σ12 = 0.5; σ13 = 0.5; σ21 = 0.5; σ22 = 1; σ23 = 0.8; σ31 = 0.5; σ32 = 0.8; σ33 = 1. Direct calculation demonstrates that λ3z = –1.1692, λ3x = 2.2482, λ3y = 0.4026. The species z(t) may die out or three species coexist. Sample paths of xw0(t), yw0(t) and zw0(t) with w0 = (3, 3, 3) are illustrated in Fig. 2.

Figure 2

Sample paths of xw0(t), yw0(t) and zw0(t) of Example 5.2 with w0 = (3, 3, 3) in two trials.

## 5 Conclusions

This paper focuses on providing the conditions for a three-species stochastic competitive model. The obtaining results show that three values λ3z, λ3x and λ3y calculated by the coefficients can be served as threshold values. For model (1.3), three species will coexist and all positive solutions converge to a unique invariant probability measure when λ3z, λ3x and λ3y > 0. Besides, compared with the system in [18], the results in this paper is in agreement with Theorem 2.1 in [18] by chosen that r3 = a31 = a32 = a33 = σ31 = σ32 = σ33 = 0.

It is noteworthy that in Section 3, Δi = 0, i = 1, 2, 3. However, for model (1.2), that is Δi ≠ 0, i = 1, 2, 3, equation on the x-axis, y-axis and z-axis are

dφ(t)=φ(t)(r1a11φ(t))dt+(Δ1φ(t)+σ11φ2(t))dB1(t), (5.1)
dψ(t)=ψ(t)(r2a22ψ(t))dt+(Δ2ψ(t)+σ22ψ2(t))dB4(t) (5.2)

and

dρ(t)=ρ(t)(r3a33ρ(t))dt+(Δ3ρ(t)+σ33ρ2(t))dB6(t), (5.3)

respectively. Employing Theorem 3.1 of [22], one can be obtained that if riΔi22<0 (i = 1, 2, 3), then Plimtφ(t)=0orlimtψ(t)=0orlimtρ(t)=0=1 for all positive solutions φ(t), ψ(t) and ρ(t). Then, by using arguments similar to the proof of Proposition 4.1 in [18], it shows that Plimtx(t)=0orlimty(t)=0orlimtz(t)=0=1. In case riΔi22>0 (i = 1, 2, 3), (5.1), (5.2) and (5.3) respectively have unique invariant probability measure whose density f~1,f~2 and f~3 can be solved from the Fokker-Planck equation. Define

λ~3z=r30(a31x+σ3122x2)f~1(x)dx0(a32z+σ3222z2)f~2(y)dy,λ~3x=r10(a12x+σ1222x2)f~2(y)dy0(a13z+σ1322z2)f~3(z)dz

and

λ~3y=r20(a21x+σ2122x2)f~1(x)dx0(a23z+σ2322z2)f~3(z)dz.

If riΔi22>0 (i = 1, 2, 3). Using the argument in Section 3 with modifications, one can be also gotten that if λ͠3z, λ͠3z and λ͠3z > 0, then model (1.2) has an invariant probability measure in R+3,o . It implies that the results stated in Theorem 3.1 hold for (1.2) with λ3z, λ3z, λ3z replaced by λ͠3z, λ͠3x, λ͠3y. Furthermore, the results recover the main findings that Theorem 3.5 and 4.3 in [17].

Though this paper has verified that three species coexist when λ3z, λ3z and λ3z > 0 for model (1.3), three species may also coexist in the case that λ3z, λ3z and λ3z are not all positive. It therefore still required novel techniques to deal with this problem.

## Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 11561069), the Guangxi Natural Science Foundation of China (Grant No. 2018GXNSFDA281028, 2017GXNSFAA198234) and the Guangxi University High Level Innovation Team and Distinguished Scholars Program of China (Document No. [2018] 35).

## References

[1] Lotka A. J., Contribution to the theory of periodic reactions, J. Trauma. Stress, 1910, 21, 271–274.10.1021/j150111a004Search in Google Scholar

[2] Lotka A. J., Analytical note on certain rhythmic relations in organic systems, Proceedings of the National Academy of Sciences of the United States of America, 6, 1920, 410–415.10.1073/pnas.6.7.410Search in Google Scholar PubMed PubMed Central

[3] Volterra V., Variations and fluctuations of the number of individuals in animal species living together, Animal Ecology Mcgrawhill, 1928, 3, 3–51.10.1093/icesjms/3.1.3Search in Google Scholar

[4] Hering R. H., Oscillations in Lotka-Volterra systems of chemical reactions, J. Math. Chem., 1990, 5, 197–202.10.1007/BF01166429Search in Google Scholar

[5] Lakka S., Michalakelis C., Varoutas D., Martakos D., Competitive dynamics in the operating systems market: Modeling and policy implications, Technol. Forecast. Soc. Change, 2013, 80, 88–105.10.1016/j.techfore.2012.06.011Search in Google Scholar

[6] Xiao D., Li W., Limit cycles for the competitive three dimensional Lotka-Volterra system, J. Differential Equations, 2000, 164, 1–15.10.1006/jdeq.1999.3729Search in Google Scholar

[7] Faria T., Global dynamics for Lotka-Volterra systems with infinite delay and patch structure, Appl. Math. Comput., 2014, 245, 575–590.10.1016/j.amc.2014.08.009Search in Google Scholar

[8] Zeeman M. L., Hopf bifurcations in competitive three dimensional Lotka-Volterra systems, Dynam. Stabil. Syst., 1993, 8, 189–217.10.1080/02681119308806158Search in Google Scholar

[9] Zeeman E. C., Zeeman M. L., An n-dimensional competitive Lotka–Volterra system is generically determined by the edges of its carrying simplex, Nonlinearity, 2002, 15, 2019–2032.10.1088/0951-7715/15/6/312Search in Google Scholar

[10] Baigent S., Geometry of carrying simplices of 3-species competitive Lotka-Volterra systems, Nonlinearity, 2013, 26, 1001–1029.10.1088/0951-7715/26/4/1001Search in Google Scholar

[11] Soliman A. A., Al-Jarallah E. S., Asymptotic stability of solutions of Lotka-Volterra predator-prey model for four species, Appl. Math., 2015, 6, 684–693.10.4236/am.2015.64063Search in Google Scholar

[12] Baigent S., Convexity of the carrying simplex for discrete-time planar competitive Kolmogorov systems, J. Difference Equ. Appl., 2016, 22, 1–14.10.1080/10236198.2015.1125895Search in Google Scholar

[13] Liu M., Fan M., Permanence of stochastic Lotka-Volterra systems, J. Nonlinear Sci., 2017, 27, 425–452.10.1007/s00332-016-9337-2Search in Google Scholar

[14] Jiang D., Zhang Q., Hayat T., Alsaedi A., Periodic solution for a stochastic non-autonomous competitive Lotka-Volterra model in a polluted environment, Phys. A, 2017, 471, 276–287.10.1016/j.physa.2016.12.008Search in Google Scholar

[15] Mao X., Marion G., Renshaw E., Environmental Brownian noise suppresses explosions in population dynamics, Stochastic Processes and their Applications, 2002, 97, 95–110.10.1016/S0304-4149(01)00126-0Search in Google Scholar

[16] Zhang Q., Jiang D., The coexistence of a stochastic Lotka-Volterra model with two predators competing for one prey, Appl. Math. Comput., 2015, 269, 288–300.10.1016/j.amc.2015.07.054Search in Google Scholar

[17] Zhao Y., Yuan S., Ma J., Survival and stationary distribution analysis of a stochastic competitive model of three species in a polluted environment, Bulletin of Mathematical Biology, 2015, 77, 1285–1326.10.1007/s11538-015-0086-4Search in Google Scholar

[18] Nguyen D. H., Yin G., Coexistence and exclusion of stochastic competitive Lotka-Volterra models, J. Differential Equations, 2017, 262, 1192–1225.10.1016/j.jde.2016.10.005Search in Google Scholar

[19] Mao X., Sabais S., Renshaw E., Asymptotic behavior of stochastic Lotka-Volterra model, J. Math. Anal. Appl., 2003, 287, 141–156.10.1016/S0022-247X(03)00539-0Search in Google Scholar

[20] Li X., Jiang D., Mao X., Population dynamical behavior of Lotka-Volterra system under regime switching, J. Comput. Appl. Math., 2009, 232, 427–448.10.1016/j.cam.2009.06.021Search in Google Scholar

[21] Bellet L. R., Ergodic properties of Markov processes, in: Open Quantum Systems II, Springer, Berlin, Heidelberg, 2006, 1–39.10.1007/3-540-33966-3_1Search in Google Scholar

[22] Ikeda N., Watanabe S., Stochastic Differential Equations and Diffusion Processes, second edition, North-Holland Publishing Co., Amsterdam, 1989.Search in Google Scholar