# Asymptotic Behavior Analysis for a Three-Species Food Chain Stochastic Model with Regime Switching

Authors: Junmei Liu; Yonggang Ma
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

This paper discusses the asymptotic behavior of a class of three-species stochastic model with regime switching. Using the Lyapunov function, we first obtain sufficient conditions for extinction and average time persistence. Then, we prove sufficient conditions for the existence of stationary distributions of populations, and they are ergodic. Numerical simulations are carried out to support our theoretical results.


---

## Body

## 1. Introduction

In recent years, the dynamic relationship between the predator and prey has become one of the research hotspots in ecology and mathematical ecology because of its universal importance. In particular, the predator-prey model is a typical inhibition model, which greatly changes the understanding of the existence and development of basic laws in the biological community. The following model is one of the Volterra models for the three species of the predator-prey system:(1)x˙1t=x1ta1−b11x1t−b12x2t,x˙2t=x2t−a2+b21x1t−b22x2t−b23x3t,x˙3t=x3t−a3+b32x2t−b33x3t,where x1t,x2t,andx3t denote the densities of prey, predator, and top-predator population at time t, respectively. The parameters a1,a2, and a3 are positive constants that stand for the intrinsic growth rate of the species x1t, the death rate of the species x2t, and the death rate of the species x3t, respectively. The coefficient b11,b22,andb33 are the intraspecific competition in the resource, b12andb23 represent the rate of consumption, and b21andb32 represent the contribution of prey to the growth of predator. As a matter of fact, there are many extensive studies in the literatures concerned with three-species predator-prey systems (see e.g., [1–6]). For example, Krikorian [1] considered the Volterra predator-prey model in a three-species and explained the global properties of its solution. Zhou [2] investigated the existence and global stability of the positive periodic solutions of delayed discrete food chains with omnivory. Hsu [5] considers a three-species Lotka–Volterra food web model with omnivores, which is defined as feeding on more than one nutritional level.In addition, the population system is always affected by environmental noise, which is very important for discovering the nature of random system from the biological point of view. Generally speaking, there are various types of environmental noise, e.g., white or color noise. First, let us consider a simple color noise, such as telegraph noise (seee.g., [7–13]). This kind of colored noise can be explained by the transformation between two or more environmental models, which can be considered to be different due to rainfall or nutrition in the population model. Therefore, we can model state switching through a finite-state Markov chain. Let rt be a right-continuous Markov chain on the probability space, taking values in a finite-state space S=1,2,…,N with the generator Γ=γuvN×N given by(2)Prt+Δt=vrt=u=γuvΔt+oΔt,ifu≠v,1+γuvΔt+oΔt,ifu=v,where Δt>0 and γuv is the transition rate from state u to state v and γuv≥0 if u≠v while γuu=−∑u≠vγuv. Then, we can incorporate the regime switch into the three-species food chain model (1) to obtain(3)x˙1t=x1ta1rt−b11rtx1t−b12rtx2t,x˙2t=x2t−a2rt+b21rtx1t−b22rtx2t−b23rtx3t,x˙3t=x3t−a3rt+b32rtx2t−b33rtx3t,with initial value xi0=xi0≥0,r0=ζ∈S. Here, we assume that the coefficients aik,bijki,j=1,2,3 are all positive for k∈S.Next, we consider other types of environmental noise, namely, the white noise (seee.g., [14–22]). In particular, Mao [14] showed that different structures of white noise may have different effects on the population systems; Mao et al. [15] revealed that the environmental noise can suppress a potential population explosion. So, we assume the intrinsic growth rate aik is disturbed with(4)aik⟶aik+σikω˙it,where ω˙iti=1,2,3 are independent white noises and σik is a positive constant representing the intensity of the white noise. We assume that the Brownian motion ωt is independent of the Markov chain rt. Therefore, system (3) becomes a three-species food chain stochastic model under regime switching:(5)dx1t=x1ta1rt−b11rtx1t−b12rtx2tdt+σ1rtdω1t,dx2t=x2t−a2rt+b21rtx1t−b22rtx2t−b23rtx3tdt−σ2rtdω2t,dx3t=x3t−a3rt+b32rtx2t−b33rtx3tdt−σ3rtdω3t.In this paper, we show that system (5) has the following properties:(i)
The solution starting from anywhere inR+3 will remain in R+3 with probability 1.(ii)
For any given initial valuex0∈R+3andr0=ς∈S, there exists a positive constant Kp such that the solution xt of system (5) has the following property: limsupt⟶∞Ex1t+x2t+x3tp≤Kp,t≥0,p>1.(iii)
We show that if the noise is sufficiently large, the solution to system (5) will become extinct with probability 1. This is, limt⟶∞xit=0,a.s.i=1,2,3.And we prove that the predator of system (5) will tend to extinction almost surly in some assumptions.(i)
The persistent in time average is investigated under certain conditions, namely, the solutionxt of system (5) with any initial value x0∈R+3,r0=ς∈S has the following property:(6)liminft⟶∞1t∫0tx3sds≥∑k∈Sπkb^21b^32r1k−bˇ11b^32r2k−bˇ12b^21+bˇ11bˇ22r3kbˇ11bˇ23b^32+bˇ12b^21+bˇ11bˇ22bˇ33.(ii)
In the case of noise being relatively small, there is a unique stationary distributionμ⋅,⋅ with ergodic property:(7)Plimt⟶∞∫0txspds=∫R3xpμxdx=1,for allp>0,where μx⋅=∑k∈Sμ⋅,k is the marginal stationary distribution of a solution xt of system (5). The key method used in this paper is to construct Lyapunov functions. This Lyapunov function analysis for stochastic differential equations has been used by many authors (see [9, 15]).The paper is organized as follows. In Section2, we give the unique existence and boundedness of the solution. In Section 3, we show the sufficient conditions for extinction and persistence in time average, respectively, which have closed relations with the stationary probability distribution of the Markov chain. Then, in Section 4, by using Lyapunov function, the sufficient conditions for the stationary distribution and ergodicity of the solution of system (8) are established. Finally, we illustrate our main results through an example in Section 5.


## 2. Preliminaries

Throughout this paper, unless otherwise specified, letΩ,ℱ,ℱtt≥0,P be a complete probability space with a filtration ℱtt≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while ℱ0 contains all P-null sets). Let R+3 denote the positive cone of R3, namely, R+3:=x∈R3:xi>0,i=1,2,3, and R¯+3 denote the nonnegative cone of R3, this is, R¯+3:=x∈R3:xi≥0,i=1,2,3. For convenience and simplicity in the following discussion, denote xt=x1t,x2t,x3t. If A is a vector or matrix, its transpose is denoted by AT. Let b^ij=mink∈Sbijk,bˇij=maxk∈Sbijk,σ^i=mink∈Sσik,σˇi=maxk∈Sσik,i,j=1,2,3.C1,2R¯+×R3×S denote the family of all nonnegative real-value function Vt,x,k which are continuously twice differentiable in x and once in t.Furthermore, as a standing hypothesis, we assume that the Markov chainrt is irreducible in this paper. This is very reasonable as it means that the system will switch from any regime to any other regime. This is equivalent to the condition that, for any u,v∈S, one can find finite numbers i1,i2,…,ik∈S such that γu,i1γi1,i2…γik,v>0. Note that Γ always has an eigenvalue 0. The algebraic interpretation of irreducibility is rank Γ=N−1. Under this condition, the Markov chain has a unique stationary distribution π=π1,π2,…,πN∈R1×N, which can be determined by solving the following linear equation πΓ=0 subject to ∑k=1Nπk=1πk>0,∀k∈S.For convenience, system (5) can be rewritten into the following form:(8)dxt=fxt,rtdt+gxt,rtdωt,where(9)fxt,rt=x1ta1rt−b11rtx1t−b12rtx2tx2t−a2rt+b21rtx1t−b22rtx2t−b23rtx3tx3t−a3rt+b32rtx2t−b33rtx3t,gxt,rt=diagx1tσ1rt,−x2tσ2rt,−x3tσ3rt,ωt=ω1t,ω2t,ω3tT.For any twice continuously differentiable functionVt,x,k∈C1,2R¯+×R3×S, we define L by(10)LVt,x,k=Vtt,x,k+Vxt,x,kfxt,rt+12gTxt,rtVxxt,x,kgxt,rt+∑l=1NγklVt,x,l,where(11)Vtt,x,k=∂Vt,x,k∂t,Vxt,x,k=∂Vt,x,k∂x1,∂Vt,x,k∂x2,∂Vt,x,k∂x3,Vxxt,x,k=∂Vt,x,k∂xi∂xj3×3.From the biological sense, we are only interested in the positive solutions of equation (5). Therefore, we first prove the existence and uniqueness of the global positive solution.Theorem 1.
For any given initial valuex0∈R+3,r0=ς∈S, there is a unique positive solution xt of system (5), and the solution will remain in R+3 with probability 1.Proof.
Note that the coefficients of system (5) are local Lipschitz continuous for the given initial value x0∈R+3,r0=ς∈S. So, there is a unique maximal local solution xt on t∈0,τe, where τe is the explosion time (see [8, 15]). To show this solution is global, we need to show that τe=∞ a.s. Since the initial value is positive and bounded, there is a number m0>0 large enough such that xi0∈1/m0,m0,i=1,2,3. For each integer m≥m0, define the stopping time:(12)τm=inft∈0,∞:xit∉1m,m,for somei=1,2,3,where infϕ=∞ (as usual ϕ denotes the empty set). Clearly, τm is increasing as m⟶∞. Set τ∞=limm⟶∞τm, hence τ∞≤τe a.s. If we can show that τ∞=∞ a.s., then τe=∞ a.s. and xt∈R+3 a.s. for all t≥0. In other words, to complete the proof all we need to show that τ∞=∞ a.s. If this statement is false, there is a pair of constants T>0 and ε∈0,1 such that Pτ∞≤T>ε. Hence, there is an integer m1≥m0 such that(13)Pτm≤T≥ε,for allm≥m1.
Define theC2-function V:R+3×S⟶R+ by(14)Vx,k=∑i=13cikxi−1−logxi,where c1k=b21kb32k,c2k=b12kb32k,andc3k=b12kb23k are positive constants. The nonnegativity of this function can be seen from x−1−logx≥0 on x>0. Then, from the generalized Itô formula (10), we have(15)LVx,k=b21kb32kx1−1a1k−b11kx1−b12kx2+b12kb32kx2−1−a2k+b21kx1−b22kx2−b23kx3+b12kb23kx3−1−a3k+b32kx2−b33kx3+12b21kb32kσ12k+12b12kb32kσ22k+12b12kb23kσ32k+∑l=1NγklVx,l=b21kb32k−a1k+12σ12k+b12kb32ka2k+12σ22k+b12kb23ka3k+12σ32k+b21kb32ka1k+b21kb32kb11k−b12kb32kb21kx1+b21kb32kb12k+b12kb32kb22k−b12kb23kb32k−b12kb32ka2kx2+b12kb32kb23k+b12kb23kb33k−b12kb23ka3kx3−b21kb32kb11kx12−b12kb32kx22−b12kb23kb33kx32+∑l=1NγklVx,l=ΔHx,k+∑l=1NγklVx,l.
Letq^=maxcil/cik:1≤i≤3,1≤k,l≤N. By the definition of Vx,k, for any k,l∈S, we have(16)Vx,l=∑i=13cilxi−1−logxi≤∑i=13q^cikxi−1−logxi≤q^Vx,k.
Thus,(17)∑l=1NγklVx,l≤q^∑l=1NγklVx,k.
In addition, it is straightforward to see thatlimx⟶∞Hx,k<0 and there exists a constant Mk>0 such that maxx∈R+3Hx,k≤Mk. Then, it follows from (15) and (17) that(18)LVx,k≤Mk+q^∑l=1NγklVx,k≤M∗1+Vx,k,where M∗=maxMk,q^∑l=1Nγkl is a positive constant. Then, from Itô formula (10), we have(19)EVxτm∧T,rτm∧T≤Vx0,r0+M∗Eτm∧T+M∗E∫0τm∧TVxτm∧T,rτm∧Tdt≤Vx0,r0+M∗T+M∗∫0TEVxτm∧T,rτm∧Tdt.
The Gronwall inequality [9] implies that(20)EVxτm∧T,rtm∧T≤Vx0,r0+M∗TeM∗T.
SetΩm=τm≤T for m≥m1 and by (13), we know PΩm≥ε. Note that, for every ω∈Ωm, there is some xiτm,ω which equals either m or 1/m; hence, Vxτm,ω is no less than c^m−1−logm or c^1/m−1−logm, where c^=mink∈Scik,i=1,2,3.
Consequently,(21)Vxτm,ω,rtm,ω≥c^m−1−logm∧1m−1−logm.
It then follows from (20) that(22)Vx0,r0+M∗TeM∗T≥E1ΩmωVxτm,ω,rtm,ω≥εc^m−1−logm∧1m−1−logm,where 1Ωm is the indicator function of Ωm. Letting m⟶∞ leads to the contradiction(23)∞>Vx0,r0+M∗TeM∗T=∞.
This completes the proof.Theorem 2.
For any given initial valuex0∈R+3,r0=ς∈S, there is a constant Kp>0 such that the solution of system (5) satisfies(24)lim supt⟶∞Ex1t+x2t+x3tp≤Kp,t≥0,p>1.Proof.
By Theorem1, the solution xt will remain in R+3 for all t≥0 with probability 1.
Define a function:(25)Vx,k=c1kx1t+c2kx2t+c3kx3t,where c1k=b21kb32k,c2k=b12kb32k,andc3k=b12kb23k are positive constants. By virtue of the generalized Itô’s formula (10), we have(26)dVx,k=c1kx1ta1k−b11kx1t−b12kx2tdt+σ1kdω1t+c2kx2t−a2k+b21kx1t−b22kx2t−b23kx3tdt−σ2kdω2t+∑l=1NγklVx,ldt+c3kx3t−a3k+b32kx2t−b33kx3tdt−σ3kdω3t,
This is(27)dVx,k=c1ka1kx1t−c2ka2kx2t−c3ka3kx3t+∑l=1NγklVx,ldt+−c1kb11kx12t−c2kb22kx22t−c3kb33kx32tdt+c1kσ1kx1tdω1t−c2kσ2kx2tdω2t−c3kσ3kx3tdω3t,
Furthermore, for any given positive constantp>1, we have(28)dVpx,k=pVp−1x,kc1ka1kx1t−c2ka2kx2t−c3ka3kx3t+∑l=1NγklVx,ldt+pVp−1x,k−c1kb11kx12t−c2kb22kx22t−c3kb33kx32tdt+pVp−1x,kc1kσ1kx1tdω1t−c2kσ2kx2tdω2t−c3kσ3kx3tdω3t+12pp−1Vp−2x,kc12kσ12kx12t+c22kσ22kx22t+c32kσ32kx32tdt+∑l=1NγklVpx,ldt.
Note that(29)c1kb11kx12t+c2kb22kx22t+c3kb33kx32t≥minb^11,b^22,b^33c1k+c2k+c3kc1kx1t+c2kx2t+c3kx3t2=minb^11,b^22,b^33c1k+c2k+c3kV2x,k.
Similar to (17), we obtain(30)∑l=1NγklVpx,l≤q^∑l=1NγklVpx,k.
Combining (28)–(30), we can obtain(31)dVpx,k≤pVp−1x,ka1kVx,k−minb^11,b^22,b^33c1k+c2k+c3kV2x,kdt+12pp−1Vp−2x,kmaxσˇ12,σˇ22,σˇ32V2x,kdt+pVp−1x,kc1kσ1kx1tdω1t−c2kσ2kx2tdω2t−c3kσ3kx3tdω3t+pq^∑l=1NγklVpx,k+q^p∑l=1NγklVpx,kdt≤pa1k+12p−1maxσˇ12,σˇ22,σˇ32+q^+q^pp∑l=1NγklVpx,k−pminb^11,b^22,b^33c1k+c2k+c3kVp+1x,kdt+pVp−1x,kc1kσ1kx1tdω1t−c2kσ2kx2tdω2t−c3kσ3kx3tdω3t.
Hence,(32)EdVpx,kdt≤pa1k+12p−1maxσˇ12,σˇ22,σˇ32+q^+q^pp∑l=1NγklEVpx,k−pminb^11,b^22,b^33c1k+c2k+c3kEVp+1x,k≤pa1k+12p−1maxσˇ12,σˇ22,σˇ32+q^+q^pp∑l=1NγklEVpx,k−pminb^11,b^22,b^33c1k+c2k+c3kEVpx,kp+1/p.
By comparison theorem [9], we obtain(33)limsupt⟶∞EVpx,k≤2a1k+p−1maxσˇ12,σˇ22,σˇ32+q^+q^p/p∑l=1Nγklc1k+c2k+c3kminb^11,b^22,b^33p=ΔLp,which implies that there is a T>0 such that(34)Ex1t+x2t+x3tp≤Lpmink∈Sc1k,c2k,c3k,t>T.
In addition,Ex1t+x2t+x3tp is continuous and there exists a ιp such that(35)Ex1t+x2t+x3tp≤ιp,t∈0,T.
LetKp=maxLp/mink∈Sc1k,c2k,c3k,ιp, then(36)limsupt⟶∞Ex1t+x2t+x3tp≤Kp,t≥0,p>1.
The proof is complete.


## 3. Extinction and Persistence in Time Average

In the previous section, we have proved that the solution of system (5) with a positive initial value remains in the positive cone R+3. In order to further study asymptotic properties of system (5), in this section, we investigate the persistence in time average and extinction of system (5) under a certain condition. We first give some assumptions and related definitions.Assumption 1.
b^12≥bˇ21,b^23≥bˇ32.Assumption 2.
∑k∈Sπka1k−1/2σ12k>0,∑k∈Sπka1k−1/2σ12k−b^11/bˇ21a2k+1/2σ22k<0.Assumption 3.
r1k=a1k−1/2σ12k,rik=aik+σ22k/2,i=2,3,(37)∑k∈Sπkr1k>0,∑k∈Sπkr1k−bˇ11b^21r2k−bˇ12b^21+bˇ11bˇ22b^21b^32r3k>0.Definition 1.
The three species are extinct iflimsupt⟶∞xit=0,a.s.i=1,2,3.
The three species are persistent in time average if.liminft⟶∞1/t∫0tx3sds>0,a.s.Theorem 3.
If Assumption1 holds, then for any given initial value x0∈R+3,r0=ς∈S, the solution xt of system (5) has the property that(38)limsupt⟶∞logx1t+x2t+x3tt≤∑k∈Sπia1k−16σ12k∧σ22k∧σ32k.Particularly, if∑k∈Sπia1k−1/6σ12k∧σ22k∧σ32k<0 holds, then(39)limt⟶∞xit=0a.s.i=1,2,3.Proof.
By the generalized Itô formula (10), we yield(40)dlogx1t+x2t+x3t=x1tx1t+x2t+x3ta1k−b11kx1t−b12x2tdt+x2tx1t+x2t+x3t−a2k+b21kx1t−b22kx2t−b23kx3tdt+x3tx1t+x2t+x3t−a3k+b32kx2t−b33kx3tdt−12σ12kx12t+σ22kx22t+σ32kx32tx1t+x2t+x3t2dt+σ1kx1tx1t+x2t+x3tdω1t−σ2kx2tx1t+x2t+x3tdω2t−σ3kx3tx1t+x2t+x3tdω3t.
By Assumption1 and basic inequality x1t+x2t+x3t2≤3x12t+x22t+x32t, integrating (40) from 0 to t, we obtain(41)logx1t+x2t+x3t≤logx10+x20+x30+∫0ta1k−16σ12k∧σ22k∧σ32kdt+M1t−M2t−M3t.where(42)M1t=∫0tσ1kx1tx1t+x2t+x3tdω1t,M2t=∫0tσ2kx2tx1t+x2t+x3tdω2t,M3t=∫0tσ3kx3tx1t+x2t+x3tdω3t.
Noticing that(43)M1,M1=∫0tσ1kx1tx1t+x2t+x3t2ds≤maxk∈Sσ1kt,by the strong law of large numbers for martingales [9, 15], we therefore have(44)limt⟶∞M1tt=0,a.s.
Similarly,limt⟶∞Mit/t=0,a.s.i=2,3.
It finally follows from (41), by dividing by t on both sides and then letting t⟶∞, and we obtain(45)limsupt⟶∞logx1t+x2t+x3tt≤limsupt⟶∞1t∫0ta1k−16σ12k∧σ22k∧σ32kdt=∑k∈Sπia1k−16σ12k∧σ22k∧σ32k.
So, we obtain the desired assertion.
Let us consider the following system:(46)dϕt=ϕta1k−b11kϕtdt+σ1kϕtdω1t,ϕ0=x10>0.
According to Luo [11], we have shown that(47)limt⟶∞ϕt=0,a.s.if∑k∈Sπka1k−12σ12k<0,limsupt⟶∞1t∫0tϕsds≤∑k∈Sπka1k−1/2σ12kb^11,a.s.if∑k∈Sπka1k−12σ12k>0.
By the comparison principle, we know thatx1t≤ϕt,t≥0.Theorem 4.
If Assumption2 holds, then for any initial value x0∈R+3,r0=ς∈S, the predator of system (5) will tend to extinction almost surly.Proof.
Note that(48)dx1t≤x1ta1k−b11kx1tdt+σ1kx1tdω1t.
By the comparison principle and (47), we have(49)limsupt⟶∞1t∫0tx1sds≤∑k∈Sπka1k−1/2σ12kb^11,a.s.
Taking Itô formula to the second equation of system (5), we obtain(50)dlogx2t=−a2k−σ22k2+b21kx1t−b22kx2t−b23kx3tdt−σ2kdω2t≤−a2k−σ22k2+b21kx1tdt−σ2kdω2t.
Integrating both sides of (50) results in(51)logx2t−logx20t≤1t∫0t−a2k−σ22k2+b21kx1sds−σ2kω2tt.
Merging (49) to (51), by Assumption 2 and limt⟶∞ω2t/t=0, we obtain(52)limsupt⟶∞logx2tt≤∑k∈Sπk−a2k−σ22k2+bˇ21ka1k−1/2σ12kb^11<0.
That is,(53)limt⟶∞x2t=0,a.s.
Similarly,(54)limsupt⟶∞logx3tt<0.
This implies(55)limt⟶∞x3t=0,a.s.
The result is confirmed.Remark 1.
It can be seen from the proof process of Theorem4 that when the following conditions are met(56)∑k∈Sπk−a2k−σ22k2+bˇ21ka1k−1/2σ12kb^11>0.
Speciesx1 and x2 coexist and species x3 is extinct.Lemma 1.
If Assumption3 holds, for any initial value x0∈R+3,r0=ς∈S, the solution xt of system (5) has following property:(57)limsupt⟶∞logxitt=0,i=1,2,3.Proof.
The proof is a modification of that of Mao [16], Zhu and Yin [12], or Luo [11]. So, we omit it.Theorem 5.
If Assumption3 holds, the solution xt of system (5) with any initial value x0∈R+3,r0=ς∈S has the following property:(58)liminft⟶∞1t∫0tx3sds≥∑k∈Sπkb^21b^32r1k−bˇ11b^32r2k−bˇ12b^21+bˇ11bˇ22r3kbˇ11bˇ23b^32+bˇ12b^21+bˇ11bˇ22bˇ33=x3∗,where x∗=x1∗,x2∗,x3∗ is the only nonnegative solution of the following equation:(59)a1k−σ12k2−b11kx1t−b12x2t=0,−a2k−σ22k2+b21kx1t−b22kx2t−b23kx3t=0,−a3k−σ32k2+b32kx2t−b33kx3t=0.Proof.
From system (5), we have(60)dc1logx1t+c2logx2t+c3logx3t=c1a1k−σ12k2−b11kx1t−b12x2tdt+σ1kdω1t+c2−a2k−σ22k2+b21kx1t−b22kx2t−b23kx3tdt−σ2kdω2t+c3−a3k−σ32k2+b32kx2t−b33kx3tdt−σ3kdω3t≥c1r1k−c2r2k−c3r3k−c2bˇ23+c3bˇ33x3tdt+c1σ1kdω1t−c2σ2kdω2t−c3σ3kdω3t,where c1=b^21,c2=bˇ11andc3=bˇ12b^21+bˇ11bˇ22/b^32.
By Assumption3, we know(61)c1r1k−c2r2k−c3r3k>0.
Integrating the both sides of (60) and dividing by t on both sides, we yield(62)c1logx1t−logx10+c2logx2t−logx20+c3logx3t−logx3tt≥1t∫0tc1r1k−c2r2k−c3r3kdt−c2bˇ23+c3bˇ33∫0tx3tdt+c1σ1kdω1t−c2σ2kdω2t−c3σ3kdω3tt.
According to Lemma1 and limt⟶∞ωit/t=0,i=1,2,3,(63)liminft⟶∞1t∫0tx3sds≥∑k∈Sπkc1r1k−c2r2k−c3r3kc2bˇ23+c3bˇ33=∑k∈Sπkb^21b^32r1k−bˇ11b^32r2k−bˇ12b^21+bˇ11bˇ22r3kbˇ11bˇ23b^32+bˇ12b^21+bˇ11bˇ22bˇ33=x3∗.
The result is confirmed.


## 4. Stationary Distribution and Ergodicity

From Theorem3, we know that three-species tend to extinct almost surely if the noise intensity is sufficiently large. However, in the section, our main aim here is to find out what happen if the noise is relatively small. Namely, we shall discuss the sufficient conditions under which system (5) admits a unique ergodic stationary distribution. To this end, consider system (8):(64)dxt=fxt,rtdt+gxt,rtdωt,where xt is a homogeneous Markov process in 3-dimension Euclidean space R3. The diffusion matrix Ax,k=gx,kgTx,k.Next, we give the conditions for the unique existence and of stationary distribution (see [13, 23–26]).Lemma 2 (see [23]).
System (8) is positive recurrent if the following conditions are satisfied:(1)
Foru≠v,γuv>0.(2)
For eachk∈S,(65)λξ2≤ξTAx,kξ≤λ−1ξ2,for allξ∈R3,with some constantλ∈0,1 for all x∈R3,(3)
There exists a bounded open subsetD of R3 with a smooth boundary satisfying that, for each k∈S there exists a nonnegative function function V⋅,k:Dc⟶R such that V⋅,k is twice continuously differentiable and that for some κ>0(66)LVx,k≤−κ,for anyx,k∈Dc×S.Moreover, the positive recurrent processxt,rt has a unique ergodic stationary distribution μ. That is, if Zxt,rt be a function integrable with respect to the measure μ, then(67)Plimt⟶∞1t∫0tZxs,rsds=∑k∈S∫R3Zx,kdx=1.We assume that the following condition holds:Assumption 4.
(1) u≠v,γuv>0,u,v∈S.(68)2Π=∑k∈Sπkr1k−r2k−r3k−a1k+b11k−b21k24b11k−−a2k+b12k+b22k−b32k24b22k−−a3k+b23k+b32k24b33k>0.
Then, we can prove.Theorem 6.
If Assumption1 and 4 are satisfied, then system (8) with initial value in R+3×S is positive recurrent and has a unique stationary distribution μ⋅,⋅ with ergodic property:(69)Plimt⟶∞∫0txspds=∫R3xpμxdx=1,for allp>0,where μx⋅=∑k=1Nμ⋅,k is the marginal stationary distribution of a solution xt of system (8).Proof.
It is easy to see if all conditions in Lemma2 are satisfied, then system (8) is positive recurrent. Obviously, we only need to verify (3) and (5). Recall that diffusion matrix is Ax,k=diagσ12x12,σ22x22,σ32x32, Let σ^=mink∈Sσi2,i=1,2,3andσˇ=maxk∈Sσi2,i=1,2,3. Hence, we have(70)λξ2≤ξTAx,kξ≤λ−1ξ2,where λ=minσ^,σˇ−1,1, which verifies condition (3) in Lemma 2. To complete the proof, considering the following bounded open subset:(71)Eρ=1ρ,ρ×1ρ,ρ×1ρ,ρ⊂R+3,where ρ is a sufficiently large number. Then, E¯ρ⊂R+3.
Define the functionV:R+3×S⟶R+ by(72)Vx,k=∑13xi−1−logxi+ϖk+ϖ=V1x+V2k,x,k∈R+3−E¯ρ×S.Here, ϖ=ϖ1,ϖ2,…,ϖNT,|ϖ|=ϖ12+⋯+ϖN2, and ϖkk∈S are to be determined later and the reason for ϖ being here is to make ϖk+ϖ nonnegative.
ForV1, we calculate the following operator to obtain(73)LV1x=x1t−1a1k−b11kx1t−b12x2t+12σ12k+x2t−1−a2k+b21kx2t−b22kx2t−b23x3t+12σ22k+x3t−1−a3k+b32kx2t−b33x3t+12σ32k.
Using Assumption1, we have(74)LV1x≤−b11kx1t−a1k+b11k−b21k2b11k2−b22kx2t−−a2k+b12k+b22k−b32k2b22k2−b33kx3t−−a3k+b33k+b23k2b33k2+a1k+b11k−b21k24b11k+−a2k+b12k+b22k−b32k24b22k+−a3k+b33k+b23k24b33k−r1k+r2k+r3k,where r1k=a1k−σ12k/2>0,rik=aik+σ22k/2,i=2,3.
On the contrary,(75)LV2k=∑u≠v∈Sϖu−ϖv.
Next, we define a vectorΞ=Ξ1,Ξ2,…,ΞNT with(76)Ξk=r1k−r2k−r3k−a1k+b11k−b21k24b11k−−a2k+b12k+b22k−b32k24b22k−−a3k+b33k+b23k24b33k.
Since the generator matrixΓ is irreducible and Lemma 2.3 in [13], for Ξk there exists a solution,(77)ϖ=ϖ1,ϖ2,…,ϖNT,for the following system:(78)Γϖ−Ξk=−∑k=1NπkΞk.
Thus, we obtain(79)∑u≠v∈Sϖu−ϖv−r1k−r2k−r3k−a1k+b11k−b21k24b11k−−a2k+b12k+b22k−b32k24b22k−−a3k+b33k+b23k24b33k=−∑k∈Sπkr1k−r2k−r3k−a1k+b11k−b21k24b11k−−a2k+b12k+b22k−b32k24b22k−−a3k+b33k+b23k24b33k.
Combining (74)–(79), we obtain(80)LVx≤−Π−b11kx1t−a1k+b11k−b21k2b11k2−b22kx2t−−a2k+b12k+b22k−b32k2b22k2−b33kx3t−−a3k+b33k+b23k2b33k2=ΔΦx,k.
It then follows from Assumption4 that(81)limx|⟶0Φx,k=limx|⟶∞Φx,k=−∞.
Therefore, from (81), there exists a sufficiently small ρ such that(82)LVx,k≤−1,for allx,k∈R+3−E¯ρ×S.
It then follows Lemma2 that the solution of system (8) is positive recurrent with respect to the domain Eρ and has a unique stationary distribution. Moreover, the ergodic property follows from the moments estimation of Theorem 2.


## 5. Numerical Simulations

In order to verify the results above, we numerically simulate the solution of the system (5). By Euler–Maruyama scheme [27, 28], we use discretized Brownian paths over 0,T and write efficient Matlab code, then obtain the simulation figures.Example 1.
Let us assume that the Markov chainrt is on the state space S=1,2 with the generator Γ=−446−6. Then, we can get the stationary distribution of rt is π=0.6,0.4. Take the step size Δt=0.005 and the following setting for parameters.Case 1.
Whenrt=1,(83)a11=0.01,b111=0.2,b121=0.6,σ11=0.3,a21=0.8,b211=0.3,b221=0.2,b231=0.8,σ21=0.28,a31=0.5,b321=0.8,b331=0.1,σ31=0.3,when rt=2,(84)a11=0.02,b111=0.2,b121=0.5,σ11=0.29,a21=0.6,b211=0.4,b221=0.2,b231=0.8,σ21=0.3,a31=0.1,b321=0.6,b331=0.1,σ31=0.3,
The above parameters satisfy the condition of Theorem3, that is, b^12=0.3>bˇ21=0.2,∑k∈Nπia1k−1/6σ12k∧σ22k∧σ32k=−0.0026<0. As can be seen from Figure 1, when Assumption 1 holds, all three species are extinct.Figure 1
Case1: the sample path of the overall system with the initial data x0=1,0.8,0.6.Case 2.
Whenrt=1,(85)a11=0.5,b111=0.2,b121=0.6,σ11=0.2,a21=0.8,b211=0.3,b221=0.8,b231=0.8,σ21=0.28,a31=0.5,b321=0.8,b331=0.1,σ31=0.3,when rt=2,(86)a11=0.5,b111=0.2,b121=0.5,σ11=0.25,a21=0.9,b211=0.3,b221=0.2,b231=0.8,σ21=0.3,a31=0.1,b321=0.6,b331=0.1,σ31=0.3,
It is easy to calculate(87)∑k∈Nπia1k−12σ12k=0.4755>0,∑k∈Nπia1k−12σ12k−b^11bˇ21a2k+12σ22k=−0.11218<0.
Figure2 shows that, under the condition of Assumption 2, predators are almost inevitably extinct, and only the lowest predators survive.Figure 2
Case2:the sample path of the overall system with the initial data x0=1,0.8,0.6.Case 3.
Whenrt=1,(88)a11=0.9,b111=0.2,b121=0.6,σ11=0.15,a21=0.3,b211=0.3,b221=0.2,b231=0.8,σ21=0.1,a31=0.2,b321=0.8,b331=0.1,σ31=0.1,when rt=2,(89)a11=0.95,b111=0.2,b121=0.5,σ11=0.25,a21=0.2,b211=0.3,b221=0.2,b231=0.8,σ21=0.12,a31=0.1,b321=0.7,b331=0.1,σ31=0.1,
It is easy to see that(90)∑k∈Nπia1k−12σ12k=0.90075>0,∑k∈Sπkr1k−bˇ11b^21r2k−bˇ12b^21+bˇ11bˇ22b^21b^32r3k=0.5689608>0.
Figure3 shows that when Assumption 3 holds, the three species live for an average of time.
In order to further expound the influence of the regime switching on the stationary distribution of system (5), we give the following example.Figure 3
Case3:the sample path of the overall system with the initial data x0=1,0.8,0.6.Example 2.
Let us assume that the Markov chainrt is on the state space S=1,2 with the generator Γ=−555−5 and setting is as follows.
Whenrt=1,(91)a11=0.9,b111=0.6,b121=0.6,σ11=0.15,a21=0.1,b211=0.5,b221=0.2,b231=0.8,σ21=0.1,a31=0.1,b321=0.2,b331=0.8,σ31=0.1,when rt=2,(92)a11=0.95,b111=0.6,b121=0.5,σ11=0.25,a21=0.1,b211=0.5,b221=0.2,b231=0.8,σ21=0.12,a31=0.1,b321=0.1,b331=0.8,σ31=0.1.
Then, we can get the stationary distribution ofrt is π=0.5,0.5. It is easy to check that Assumption 4 is satisfied and π=0.1229>0. The stationary distribution of the Markov chain implies that the populations have an equal chance living in the component-wise environment 1 or 2. Figure 4 shows that the stable distribution of Markov chain shows that species x1 can survive in different environments. In Figure 4(a), the blue line, red line, and black line, respectively, represent the density track of species x1 under two switching states. In Figure 4(b), the black line, red line, and blue line, respectively, represent the probability density function curve of species x1 in two switching states. Figure 5 shows the density trajectory and probability density function curve of predator x2 under Markov chain switching. Figure 6 shows the density trajectory and probability density function curve of predator x3 under Markov chain switching.Figure 4
(a) The black line denotes the sample paths of the population densityx1, the red line denotes the corresponding component-wise rt=1 of x1, and the blue line denotes the corresponding another component-wise rt=2 of x1, respectively. (b) The PDF of x1 and its corresponding component-wise.Figure 5
(a) The black line denotes the sample paths of the population densityx2, the red line denotes the corresponding component-wise rt=1 of x2, and the blue line denotes the corresponding another component-wise rt=2 of x2, respectively. (b) The PDF of x2 and its corresponding component-wise.Figure 6
(a) The black line denotes the sample paths of the population densityx3, the red line denotes the corresponding component-wise rt=1 of x3, and the blue line denotes the corresponding another component-wise rt=2 of x3, respectively. (b) The PDF of x3 and its corresponding component-wise.


## 6. Conclusion

Three-species food chain model has been studied by many scholars recently. In particular, many asymptotic estimates on the sample average in time have been obtained (see e.g., [8, 29]). However, to the best of our knowledge, there is rare result about the stationary distribution of a three-species food chain stochastic model under regime switching. In this paper, we develop and analyse a three-species food chain stochastic model, which takes both white and colored noises into account. We first prove the existence of the global positive solution of the model. Then, using the stochastic Lyapunov functions, we investigate extinction in probability and persistence in time average. Furthermore, we obtain a stationary distribution of the solution. Moreover, some interesting questions deserve further investigation, such as incorporating intervention strategies into the system. We leave this for future consideration.



---

*Source: 1679018-2020-08-31.xml*