# Analysis of Stochastic Nicholson-Type Delay System under Markovian Switching on Patches

Authors: Xin Yi; Guirong Liu
Journal: Complexity
 (2020)

---

## Abstract

Based on the influence of random environmental perturbations and the patch structure, we propose a stochastic Nicholson-type delay system under Markovian switching on patches. Existence of a global positive solution is studied. Then, we show ultimate boundedness and estimation of the sample Lyapunov exponent of the solution. Furthermore, sufficient conditions for extinction of species are established, which is the main new ingredient of this paper. Finally, some numerical examples are presented. Our results improve and generalize previous related results.


---

## Body

## 1. Introduction

In 1980, Gurney et al. [1] established Nicholson’s blowflies equation according to experimental data of Nicholson [2]. In recent decades, there have been a large amount of results related to the dynamical behaviors for this model and its modification, see [3–13].In ecosystems, the pattern of complex population dynamics is inevitably subject to some kind of environmental noises. As a matter of fact, the phenomenon of stochasticity plays a critical role in understanding the evolutionary dynamics and ecological characteristics of species. Particularly, May [14] has revealed that due to environmental fluctuations, the parameters in a system should be stochastic. Environmental noises are classified into two categories: the first is white noise, and the second one is coloured noise. Stochastic population models [15–20] are more realistic compared to deterministic population models. Wang et al. [21] first studied a scalar stochastic Nicholson’s blowflies delayed equation(1)dxt=−αxt+pxt−τe−γxt−τdt+σxtdBt.Notice, however, that white noise is unable to depict the phenomena that the species may be invaded by the alien population [22] or suffer sudden catastrophic shocks [23]. And in recent years, some significant progress has been made in the theory of the stochastic population models with regime switching, see [24–27] and the references therein. In [28], Zhu et al. considered a stochastic equation with Markovian switching:(2)dxt=−αrtxt+prtxt−τrte−γrtxt−τrtdt+σrtxtdBt,where continuous-time Markov chain rtt≥0 is defined on a state space S=1,2,…,m.On the contrary, migration is a ubiquitous phenomenon in the nature. Both continuous reaction-diffusion models and discrete patchy systems could incorporate and explain the phenomenology of spatial dispersion [29] in the literature of mathematical ecology. Objectively speaking, patch-structured models illustrate the spatial heterogeneity of species, depending on a lot of factors, such as ecological systems in different geographic types (e.g., nature reserves and other regions), various food-rich patches of habitats, and many other circumstances. Besides, models in the patchy environment include disease systems as well, such as the two-compartment model of the cancer cell population. In order to take the dispersal phenomenon into consideration, Berezansky et al. [30] introduced the Nicholson-type delay system on patches as follows:(3)x1′t=−a1x1t−b2x1t+b1x2t+p1x1t−τe−γ1x1t−τ,x2′t=−a2x2t−b1x2t+b2x1t+p2x2t−τe−γ2x2t−τ,which includes the novel two-compartment models of leukemia dynamics and the systems of marine protected areas.In particular, considering that the parametersai of system (3) are affected by the white noise, Yi and Liu [31] formulated the stochastic diffusion system which consists of two patches:(4)dx1t=−a1x1t−b2x1t+b1x2t+p1x1t−τe−γ1x1t−τdt+σ1x1tdB1t,dx2t=−a2x2t−b1x2t+b2x1t+p2x2t−τe−γ2x2t−τdt+σ2x2tdB2t.We can further model random shift in different regimes by a continuous-time Markov chainℓtt≥0 defined on a state space ℳ=1,2,…,N. Let ℓtt≥0 be right-continuous and Γ=ρijN×N be its generator of ℓtt≥0, i.e.,(5)ℙℓt+δ=j∣ℓt=i=ρijδ+oδ,ifj≠i,1+ρiiδ+oδ,ifj=i,where δ>0, ρij≥0 for i≠j, and ∑j∈ℳρij=0,i,j∈ℳ. Suppose that ℓtt≥0 is irreducible and has the unique stationary distribution π=π1,π2,…,πN. Hence, we obtain the stochastic Nicholson-type system under Markovian switching on the patch structure as follows:(6)dx1t=−a1ℓtx1t−b2ℓtx1t+b1ℓtx2t+p1ℓtx1t−τℓte−γ1ℓtx1t−τℓtdt+σ1ℓtx1tdB1t,dx2t=−a2ℓtx2t−b1ℓtx2t+b2ℓtx1t+p2ℓtx2t−τℓte−γ2ℓtx2t−τℓtdt+σ2ℓtx2tdB2t,with initial conditions(7)xt=φt=φ1t,φ2tT,t∈−τ,0,ℓ0=ℓ0∈ℳ,where φh∈C−τ,0;0,+∞ and φh0>0 for h=1,2 and τ=maxi∈ℳτi.We focus on the meaning of parameters with respect to fish population in marine protected areaA1 and fishing area A2. x1t and x2t are the number of fish populations in A1 and A2, respectively; for h=1,2 and i∈ℳ, a1i and a2i are the mortality rate in A1 and A2, respectively; let Gxht−τi=phixht−τie−γhixht−τi be the fish growth rates; p1i and p2i represent the maximum per adult yearly birth rate in A1 and A2, respectively; γhi>0; 1/γ1i and 1/γ2i are the number at which the reproduction at their maximum birth rate in A1 and A2, respectively; τi is the maturation time; Bht is the standard Brownian motion defined on the complete probability space Ω,ℱ,ℙ; and σhi≥0, for any i∈ℳ and h=1,2. We assume ℓtt⩾0 is ℱt-adapted. Nevertheless, suppose ℓtt⩾0 and Bht are independent of each other, h=1,2.Especially, system (6) can reduce to the model in [32] if τi≡τ,i∈ℳ. By contrast, our work differs from and improves [32], which will be depicted further in detail.In the field of ecology, it is important to use mathematics to study extinction of species, see [33, 34] and the references therein. However, no work has yet been done on the problem of extinction for scalar equation (1), not to mention the scalar equation with Markovian switching (2) and system (4). In order to prove the extinction of species, the conventional method is to construct a proper Lyapunov function or functional and then estimate the upper bound of the drift term of its Ito^ differential. Taking system (6) for example, x1t and x2t are likely to appear in the denominator of the expression of LV, and coefficients in front of them are positive, for a general Lyapunov function Vx1,x2. Unfortunately, this leads to some difficulties in finding the upper bound of LV. So, based on this, we give a new method for investigating extinction of species. Especially, system (6) reduces to (1), (2), (4), or the system in [32] when parameters of system (6) assume some special values. That is to say, we have derived extinction of the above systems at the same time.In this paper, system (6) is more general than the model of [21, 28, 30–32]. In addition, our results improve and generalize the corresponding results in these literature studies.The remainder of this paper is built up as follows. In Section2, we show the global existence of almost surely positive solution. The asymptotic estimates for the solution, stochastically ultimate boundedness, and boundedness for the average in time of the θth moment of the solution are then constructed in Section 3. In Section 4, we discuss the pathwise properties of the solution. Sufficient conditions for extinction of species are obtained in Section 5. Numerical investigations are then given in Section 6. The last part is a conclusion.


## 2. Preliminary Results

To simplify, denote the solution of (6) with initial values (7):(8)xt≔xt;φ,ℓ0,where xt=x1t,x2tT. Let(9)β1i=a1i+b2i,β2i=a2i+b1i,i∈ℳ.We denoteℝ+=0,+∞, ℝ+2=x1,x2T∈ℝ2:x1>0,x2>0, and ℝ+2×2=wuv2×2∈ℝ2×2:wuv>0,u,v=1,2. For any Φ:ℳ⟶ℝ, let Φ^=mini∈ℳΦi and Φˇ=maxi∈ℳΦi. Let · denote Euclidean norm in ℝ2. Denote the trace norm A=traceATA for matrix A.Lemma 1.
Given any initial values (7), system (6) has a unique solution xt∈R+2 for all t∈−τ,∞ almost surely.Proof.
We omit the proof since it is analogous to that of [31] by making use of the generalized Ito^ formula (see, e.g., Theorem 1.45 in [35]) to ∑h=12xh−1−logxh.Remark 1.
The delay stochastic Nicholson-type model under regime switching on patches (6) is a direct extension of the models in [21, 28, 30–32]. From Lemma 1, it is worthy to point out that priori conditions α>σ2/2 in [21] are unnecessary. Therefore, Lemma 1 improves and generalizes Lemma 2.2 in [21]. In addition, this lemma shows that both white noise and telegraph noise will not destroy a great property that the solution of (3) does not explode.


## 3. Boundedness

Because of resource constraints, asymptotic boundedness is the core of the research in ecosystems. And it is the main purpose of the present section. For simplicity, we use the following notations. For anyi∈ℳ, denote(10)K1θ,i≔a1i−θ−1θb1i−b2i−12θ−1σ12i−θ−1,K2θ,i≔a2i−θ−1θb2i−b1i−12θ−1σ22i−θ−1,Hhθ,i≔θ·supxh∈R+−Khθ,ixhθ+phiγhiexhθ−1,Hˇhθ≔maxi∈ℳHhθ,i,A1θ≔2θ/2θθ−1∑i∈ℳπiH1θ,i+H2θ,i,A2θ≔∑h=12Hˇhθθθ−1,A3θ≔2θ/2·∑h=12Hˇθθθ−1,h=1,2.Firstly, inspired by the work of Wang and Chen [32], we give this theorem.Theorem 1.
Letθ>1 such that Khθ,i>0,h=1,2,i∈ℳ. Given any initial values (7), solution x1t,x2t of (6) satisfies(11)limsupt⟶∞1t∫0tExhsθds≤A1θ,h=1,2,and(12)limsupt⟶∞Ex1θt+x2θt≤A2θ.In particular,(13)limsupt⟶∞Extθ≤A3θ.That is, system (6) is ultimately bounded.Proof.
Define(14)V1x1,x2=x1θ+x2θ.
The generalized Ito^ formula, together with the fact phiyhe−γhiyh≤phi/γhie and the elementary inequality AεB1−ε≤Aε+B1−ε for any A,B≥0 and ε∈0,1, yields(15)LV1x1,x2,y1,y2,i=−θa1ix1θ−θb2ix1θ+θb1ix1θ−1x2+12θθ−1σ12ix1θ−θa2ix2θ−θb1ix2θ+θb2ix1x2θ−1+12θθ−1σ22ix2θ+θ∑h=12phixhθ−1yhe−γhiyh≤θ−a1i+θ−1θb1i−b2i+12θ−1σ12i+θ−1x1θ+−a2i+θ−1θb2i−b1i+12θ−1σ22i+θ−1x2θ+∑h=12phiγhiexhθ−1−θθ−1x1θ+x2θ≤∑h=12θ−Khθ,ixhθ+phiγhiexhθ−1−θθ−1V1x1,x2≤∑h=12Hhθ,i−θθ−1V1x1t,x2t.
Therefore, fort>0,(16)dV1x1t,x2t≤∑h=12Hhθ,ℓt−θθ−1V1x1t,x2tdt+∑h=12θσhℓtxhθtdBht.
Then, (16) implies(17)0≤EV1x1,x2≤φ1θ0+φ2θ0+∫0t∑h=12Hhθ,ℓsds−θθ−1∫0tEV1x1s,x2sds.
Noting that the Markov chainℓt has an invariant distribution π=πi.i∈ℳ and applying the ergodic property of the Markov chain, it yields(18)limsupt⟶∞1t∫0tEV1x1s,x2sds≤limsupt⟶∞1θθ−11tV1φ10,φ20+1t∫0t∑h=12Hhθ,ℓsds≤1θθ−1∑i∈ℳπiH1θ,i+H2θ,i,a.s.
Furthermore, we have(19)deθθ−1tV1x1t,x2t≤eθθ−1t∑h=12Hhθ,ℓtdt+eθθ−1t∑h=12θσhℓtxhθtdBht≤∑h=12Hˇhθeθθ−1tdt+eθθ−1t∑h=12θσhℓtxhθtdBht.
Hence,(20)eθθ−1tEV1x1t,x2t≤V1φ10,φ20+∑h=12Hˇhθθθ−1eθθ−1t−∑h=12Hˇhθθθ−1.
Consequently, we infer immediately that (12) holds. On the contrary, according to (12), (18), and the fact that(21)xθ≤2θ/2maxx1θ,x2θ≤2θ/2V1x,it follows that (11) and (13) hold. The proof is therefore complete.Remark 2.
In Theorem1, the parameter θ is greater than 1 in the result. Although ultimate boundedness in the θth moment was derived for θ restricted to the precondition θ>1, θth moment of system (6) can be obtained when θ≤1 by Hölder’s equality.Remark 3.
Without regime switching or without migration and regime switching, Theorem1 improves the corresponding results in [21, 31]. If τi≡τ, system (6) is a direct extension of the model in [32]. Besides, no proof of ultimate boundedness in the pth moment is given in [32], which is shown in Theorem 1. Therefore, this theorem extends and improves Theorem 3.1 in [21], Theorem 2.2 in [28], Theorem 3.3 in [31], and Theorem 3.2 in [32].Theorem 2.
Given any initial values (7), solution x1t,x2t of (6) satisfies(22)limsupt⟶∞Ext≤limsupt⟶∞Ex1t+x2t≤p˜1γ^1eλ+p˜2γ^2eλ,where λ=mina^1,a^2. That is, (6) is ultimately bounded in mean.Proof.
LetV¯1t,x1,x2=eλtx1+x2. Then,(23)Ex1t+x2t≤e−λtV¯10,φ10,φ20+pˇ1γ^1e+pˇ2γ^2e∫0tes−tλds.
Finally, (22) follows by letting t⟶∞. The proof is therefore complete.Remark 4.
Compared with Theorem1, this theorem describes the case that θ=1, which does not require any conditions. If τi≡τ, we get pˇ1/γ^1eλ+pˇ2/γ^2eλ≤c/a, where c/a is defined in [32]. So, this theorem improves and extends Theorem 3.1 in [21] and Theorem 3.1 in [32].Theorem 3.
System (6) is stochastically ultimately bounded.Proof.
By (22), we derive(24)limsupt⟶∞Exht≤pˇ1γ^1eλ+pˇ2γ^2eλ,h=1,2.
By the Chebyshev inequality, it yields, for anyε∈0,1,(25)limsupt⟶∞ℙxht≥H≤H−1p⌄1γ^1eλ+p⌄2γ^2eλ=ε,where H=1/εp˜1/γ^1eλ+p˜2/γ^2eλ. The proof is therefore complete.Remark 5.
Theorem3 can be seen as the extension and improvement of [31, 32].


## 4. Asymptotic Pathwise Estimation

We shall estimate a sample Lyapunov exponent in what follows.Lemma 2.
Ifa∈ℝ and b∈ℝ+, then ax2+bx/1+x2≤Ka for x∈ℝ, where Ka=a+a2+b2/2.By the properties of quadratic functions, the proof of this lemma is easy and so is omitted. In the process of findingKa, we know that the precondition is a−Ka<0. In this case, we can choose Ka which satisfies Ka=a+a2+b2/2. We have to mention that it has no relation with the sign of parameter a. If a<0, we get a+a2+b2/2<−b2/4a by simple computation. So, this lemma is an improvement of Lemma 1.2 in [28] and Lemma 2.1 in [32].Theorem 4.
Given any initial values (7), solution xt of (6) satisfies(26)limsupt⟶∞1tlogxht≤limsupt⟶∞1tlogxt≤Q2,a.s.where h=1,2, Q=maxi∈ℳminε∈ℝ+Q1i,ϵ+Q2i,ϵ with(27)Q1i,ε=2β1i−σ12i−b1i+b2iϵ2+4p1i/γ1ie2−2β1i−σ12i−b1i+b2iϵ2,Q2i,ε=2β2i−σ22i−b1i+b2i/ε2+4p2i/γ2ie2−2β2i−σ22i−b1i+b2i/ε2,for any positive constant ε.Proof.
The generalized Ito^ formula, together with Lemma 2 and the Cauchy–Schwarz inequality, yields(28)log1+x12t+x22t≤log1+x120+x220+∫0t−2β1ℓs−σ12ℓs−b1ℓs+b2ℓs∈x12s+2p1ℓs/γ1ℓsex1s1+x12s+x22sds+∫0t−2β2ℓs−σ22ℓs−b1ℓs+b2ℓs/∈x22s+2p1ℓs/γ1ℓsex2s1+x12s+x22sds−2∑h=12σh2ℓs∫0txh4s1+x12s+x22s2ds+M1t+M2t≤log1+x120+x220+Qt−2∑h=12∫0tσh2ℓsxh4s1+x12s+x22s2ds+M1t+M2t,where for any h∈1,2,(29)Mht=2∫0tσhℓsxh2s1+x12s+x22sdBhs,with the quadratic variation(30)Mht,Mht=4∫0tσh2ℓsxh4s1+x12s+x22s2ds.
According to the exponential martingale inequality (see, e.g., [36]), for any integer m>0, we have(31)ℙsup0≤t≤mMht−2∫0tσh2ℓsxh4s1+x12s+x22s2ds>2logm≤1m2,h=1,2.
Since∑m=1∞1/m2<∞ and Borel–Cantelli’s lemma (see, e.g., [36]), there exist Ω0∈ℱ with ℙΩ0=1 and an integer m0=m0ω such that(32)Mht≤2∫0tσh2ℓsxh4s1+x12s+x22s2ds+2logm,h=1,2,for all ω∈Ω0, 0≤t≤m. Substituting the above inequality into (28), for any ω∈Ω0,m≥m0,0≤t≤m, we have(33)log1+x12t+x22t≤log1+x120+x220+Qt+4logm,which yields(34)1tlog1+x12t+x22t≤1m−1log1+x120+x220+Qm+4logm,for all ω∈Ω0,0≤m−1≤t≤m,m≥m0. Letting m⟶∞ and using the inequality y≤1/21+y2 for any y∈−∞,+∞, we obtain(35)limsupt⟶∞1tlogxht≤limsupt⟶∞12m−1log1+x20+Qm+4logm=Q2,a.s.
The proof is therefore complete.Remark 6.
Without migrations, we getQ1i=2β1i−σ12i2+4p1i/γ1ie2−2β1i−σ12i/2. By comparison, we find that Q1i≤Ci, where Ci is defined in [28]. In addition, without migration and regime switching, we can get Q in Theorem 4 is less than K, where K is defined in [21]. Furthermore, the condition 2α1−σ12−b1+b2ϵ>0,2α2−σ22−b1+b2/ϵ>0 in [31] means that the parameter ϵ needs to be satisfied: b1+b2/2α2−σ22<ε<2α1−σ12/b1+b2. However, we know that this condition is unnecessary from the above theorem. Despite all this, if we let the parameter ε satisfy ε∈b1+b2/2α2−σ2,2α12−σ12/b1+b2, we compute that Q in Theorem 4 is less than Q in [31]. Therefore, the above work is a promotion of Theorem 4.1 in [21], Theorem 2.2 in [28], and Theorem 4.1 in [31].


## 5. Extinction

Sufficient conditions for extinction are the subject of this section. Unless otherwise stated, we hypothesizeτi≡τ,i∈ℳ in this section. We first rewrite (6) as follows:(36)dxt=f1xt,xt−τ,ℓtdt+f2xt,ℓtdBt,where the operator f1:ℝ+2×ℝ+2×ℳ⟶ℝ+2 is defined as(37)f1x,y,i=−β1ix1+b1ix2+p1iy1e−γ1iy1−β2ix2+b2ix1+p2iy2e−γ2iy2,the operator f2:ℝ+2×ℳ⟶ℝ+2×2 is defined as f2x,i=σ1ix100σ2ix2, and dBt=dB1tdB2t. We first note that(38)f10,0,i≡0,f20,i≡0,for i∈ℳ, whence (6) admits a trivial solution corresponding to φ0=0.Before our result, we give a lemma.Lemma 3.
For system (36), the terms f1x,y,i and f2x,i are locally bounded in x,y while uniformly bounded in i. That is, for any m>0, there is Km>0 satisfying(39)f1x,y,i∨f2x,i≤Km,for all i∈ℳ, x,y∈ℝ+2 with x∨y≤m.The proof is not particularly difficult, so we omit the proof.Theorem 5.
Assume that(40)2β^1>σˇ12+bˇ1+1+2bˇ2+2+12pˇ1+1+2pˇ2,2β^2>σˇ22+1+2bˇ1+bˇ2+1+2pˇ1+2+12pˇ2.Then, the solution of (36) satisfies limt⟶∞xt=0,a.s., for any initial values (7). That is, all populations in system (36) go to extinction with probability one.Proof.
Step 1: let(41)V2x1,x2=xT112121x=x12+2x1x2+x22.Obviously,V2 is positive-definite and radially unbounded. That is,(42)limx|⟶∞V2x1,x2=∞.The generalized Ito^ formula yields(43)V2x1t,x2t=V2φ10,φ20+∫0tLV2xs,xs−τℓs,ℓsds+∫0t∂∂xV2x1s,x2sf2xs,ℓsdBs.By computation, we know(44)LV2x,y,i≤−2β1i−σ12i−b1i−1+2b2i−p1i−12p2ix12−2β2i−σ22i−1+2b1i−b2i−12p1i−p2ix22+1+12p1iy12+1+12p2iy22≤−λ1x12−λ2x22+1+12pˇ1y12+1+12pˇ2y22,whereλ1=2β^1−σˇ12−bˇ1−1+2bˇ2−pˇ1−1/2pˇ2 and λ2=2β^2−σˇ22−1+2bˇ1−bˇ2−1/2pˇ1−pˇ2.It is straightforward to see from (40) that λh>0,h=1,2. For simplicity, we let(45)F1x=minλ1,λ2x2,F2x=1+12maxpˇ1,pˇ2x2.By condition (40) again, we obtain that(46)Fx≔F1x−F2x=minλ1,λ2−1+12maxpˇ1,pˇ2x2>0,x≠0.Applying (44) and (46), we derive(47)∫0tLV2xs,xs−τℓs,ℓsds≤∫−τ0F2xsds−∫0tFxsds.Substituting the preceding equality into (43), it yields(48)V2x1t,x2t≤V2φ10,φ20+∫−τ0F2xsds−∫0tFxsds+∫0t∂∂xV2x1s,x2sf2xs,ℓsdBs.Then, the nonnegative semimartingale convergence theorem (see, e.g., [37]) implies(49)limsupt⟶∞V2x1t,x2t<∞a.s.Moreover, we obtain from (48) that(50)E∫0tFxsds≤V2φ10,φ20+∫−τ0F2xsds.Then, lettingt⟶∞, together with the Fubini theorem, we have(51)E∫0∞Fxtdt<∞.LetAk=ω|Yω=∫0∞Fxs,ωds>2k, where k=1,2,…. Obviously, Ak⊃Ak+1. Combining Chebyshev’s inequality and (51), we see that ∑k=1∞ℙAk<∞. By Borel–Cantelli’s lemma, one can show that ℙlimk⟶∞Ak=ℙω|Yω=∞=0, that is,(52)∫0∞Fxtdt<∞a.s.Step 2: from (52), we observe(53)liminft⟶∞Fxt=0,a.s.One now needs to consider(54)limt⟶∞Fxt=0a.s.If the above conclusion would not hold, thenℙlimsupt⟶∞Fxt>0>0. So, there is ε∈0,1/3 satisfying(55)ℙΩ1≥3ε,where(56)Ω1=limsupt⟶∞Fxt>2ε.Noting that Lyapunov functionV2xt and the solution xt of (6) are all continuous, together with (49), it yields(57)sup−τ≤t<∞V2xt<∞,a.s.Define(58)νr=infx≥kV2x,fork>0.Clearly,(59)sup−τ≤t<∞νxt≤sup−τ≤t<∞V2xt<∞,a.s.In addition, by (42), we get(60)limk⟶∞νk=∞.So,(61)sup−τ≤t<∞xt<∞,a.s.Recalling (7), we know that the initial values satisfy φh∈C−τ,0;0,+∞ for h=1,2. We therefore could find an integer m>0, depending on ε, sufficiently large for φs<m for s∈−τ,0 almost surely, while(62)ℙΩ2≥1−ε,whereΩ2=sup−τ≤t<∞xt<m. By (55) and (62), one implies(63)ℙΩ1∩Ω2≥ℙΩ1−ℙΩ2c≥2ε,whereΩ2c is the complement of Ω2. Let(64)ρ1=inft≥0:Fxt≥2ε,ρ2j=inft≥ρ2j−1:Fxt≤ε,j=1,2,…,ρ2j+1=inft≥ρ2j:Fxt≥2ε,j=1,2,…,σm=inft≥0:xt≥m.From (53) and the definitions of Ω1 and Ω2, we have(65)Ω1∩Ω2⊂σm=∞∩∩j=1∞ρj<∞.Hence, we defineζt=xt∧σm for t>−τ, and its differential is(66)dζt=f1˜tdt+f2˜tdBt,where(67)f1˜t=f1xt,xt−τ,ℓtI0,σmt,f2˜t=f2xt,ℓtI0,σmt.Here,IA is the indicator function of A. Recalling Lemma 3, we know(68)f1˜t∨f2˜t≤Km,a.s,for anyt≥−τ, ℓt∈ℳ, and xt∨xt−τ≤m. Moreover, the definition of yields ζt≤m, t≥0. We also note that, for all ω∈Ω1∩Ω2 and j≥1,(69)Fζρ2j−1−Fζρ2j=ε,Fζt≥ε,t∈ρ2j−1,ρ2j.In the close ballS¯m=x∈ℝ2:x≤m, F· is uniformly continuous. Therefore, there exists ξ=ξε>0 small sufficiently such that(70)Fζ−Fζ¯<ε,ζ,ζ¯∈S¯mwithζ−ζ¯<ξ.Forω∈Ω1∩Ω2, we emphasize that if ζρ2j−1+t−ζρ2j−1<ξ for some T>0 and t∈0,T, then ρ2j−ρ2j−1≥T. Furthermore, let the number T=Tε,ξ,m>0 be small enough such that(71)2Km2TT+4≤εξ2.By (63) and (65), we can obtain that(72)ℙρ2j<∞≥2ε.In particular, ifρ2j<∞, then ζρ2j<m. Hence, the definition of ζt implies ρ2j<σm. So,(73)ζt,ω=xt,ω,for all0≤t≤ρ2j,ω∈ρ2j<∞. Then, from the Hölder inequality and the Burkholder–Davis–Gundy inequality, it follows that(74)EIρ2j−1<∞sup0≤t≤Tζρ2j−1+t−ζρ2j−12≤EIρ2j−1<∞sup0≤t≤T∫ρ2j−1ρ2j−1+tdζs2=EIρ2j−1<∞sup0≤t≤T∫ρ2j−1ρ2j−1+tf1˜sds+∫ρ2j−1ρ2j−1+tf2˜sdBs2≤2EIρ2j−1<∞sup0≤t≤T∫ρ2j−1ρ2j−1+tf1˜sds2+∫ρ2j−1ρ2j−1+tf2˜sdBs2≤2TEIρ2j−1<∞sup0≤t≤T∫ρ2j−1ρ2j−1+tf1˜s2ds+8EIρ2j−1<∞sup0≤t≤T∫ρ2j−1ρ2j−1+tf2˜s2dBs≤2Km2TT+4,which together with (71) and the Chebyshev inequality, imply easily that(75)ℙρ2j−1<∞∩sup0≤t≤Tζρ2j−1+t−ζρ2j−1≥ξ≤ε.Obviously, we observe thatρ2j−1<∞ if ρ2j<∞. Hence, by (72) and (75), we can derive that(76)ℙρ2j<∞∩sup0≤t≤Tζρ2j−1+t−ζρ2j−1<ξ=ℙρ2j<∞−ℙρ2j<∞∩sup0≤t≤Tζρ2j−1+t−ζρ2j−1≥ξ≥ε.This, together with (70), we can get that(77)ℙρ2j<∞∩ρ2j−ρ2j−1≥T≥ε.Recalling (51), (69), and (77), for all ω∈Ω1∩Ω2 and j≥1, we compute(78)∞>E∫0∞Fxtdt≥∑j=1∞EIρ2j<∞∫ρ2j−1ρ2jFζtdt≥ε∑j=1∞EIρ2j<∞ρ2j−ρ2j−1≥εT∑j=1∞ε=∞,which is a contradiction. Consequently, we infer (54).Step 3: by (46) and (54), we now derive limt⟶∞xt=0,a.s. The proof is therefore complete.Corollary 1.
Assume that are nonnegative constants,τrt≡τ≥0, and(79)2minrt∈Sαrt>maxrt∈Sσrt2+3+322maxrt∈Sprt.Then, solutionxt of (2) obeys(80)limt⟶∞xt=0,a.s.,for any initial value xt=ϕt,t∈−τ,0,ϕ0>0,ϕ∈C−τ,0;0,+∞. That is, all populations in equation (2) go to extinction with probability one.Corollary 2.
Assume thatah,bh,ph,γh,σh,τ are nonnegative constants, h=1,2, and(81)2a1>σ12+b1+2−1b2+2+12p1+1+2p2,2a2>σ22+2−1b1+b2+1+2p1+2+12p2.Then, solutionxt of (4) obeys(82)limt⟶∞xt=0,a.s.,for any initial value xht=ϕht,t∈−τ,0,ϕh0>0,ϕh∈C−τ,0;0,+∞. That is, all populations in model (4) go to extinction with probability one.Remark 7.
This theorem reveals that the solutions of (6) will all tend to the origin asymptotically with probability one when the intensities of noises and the parameters satisfy condition (40). However, [21, 28, 31, 32] do not study extinction of populations. Besides, this method can be extended to research extinction in the above literature studies. Corollaries 1 and 2 give the conditions of extinction of (2) and (4), respectively. Therefore, our work is the extension of [21, 28, 31, 32].


## 6. Numerical Simulations

Based on [38], we show numerical simulations in the present section.Here, we consider model (6) with the same initial data φ10=1,φ20=0.5 and the same ℓtt≥0 on ℳ=1,2,3 with(83)Γ=−10462−3135−8.Then, we know the Markov chainℓtt≥0 is irreducible and has a unique stationary distribution π=0.1845,0.6019,0.2136.In Figure1, we give a simulation of the sample path of ℓtt≥0 with ℓ0=3.Figure 1
Sample path ofℓtt≥0 with ℓ0=3.In Figure2, we can choose a1=1,2,3; b1=0.12,0.23,0.22;p1=4,5,6; γ1=0.4,0.5,0.6;σ1=1.2,1.1,2; a2=3,1,2;b2=0.23,0.16,0.12; p2=7,4,3;γ2=0.3,0.4,0.5; σ2=0.45,1.5,0.25;τ=1,2,3. It is easy to see that 1=a22<σ222/2=1.125 when ξt=2. There is a good agreement between Lemma 1 and Figure 2. By Theorem 2, we know p1ˇ/γ1^eλ+p2ˇ/γ2^eλ=14.1020. Furthermore, we get Q=max9.5553,6.3244,3.8081 by calculation. Therefore, conditions of Theorem 4 have been checked. So, limsupt⟶∞1/tlnxht≤4.7777,a.s.Figure 2
Numerical solutions of (6) with a1=1,2,3; b1=0.12,0.23,0.22; p1=4,5,6; γ1=0.4,0.5,0.6; σ1=1.2,1.1,2; a2=3,1,2; b2=0.23,0.16,0.12; p2=7,4,3; γ2=0.3,0.4,0.5; σ2=0.45,1.5,0.25; and τ=1,2,3.In Figure3, we can choose a1=2.1,2.2,2;b1=0.11,0.13,0.12;p1=3.2,3,3.5;γ1=0.83,0.85,0.86;σ1=0.62,0.61,0.63;a2=2,2.1,2.2;b2=0.11,0.13,0.12;p2=3.1,3.3,3.2;γ2=0.81,0.8,0.83;σ2=0.71,0.75,0.71. Let θ=2. Then, conditions of Theorem 1 could be checked. By calculation, we get limsupt⟶∞1/t∫0tExhsθds≤31.9797,limsupt⟶∞Ex1θt+x2θt≤19.3208,limsupt⟶∞Extθ≤38.6416. Figure 3 clearly supports this result.Figure 3
Numerical solutions of (6) with a1=2.1,2.2,2; b1=0.11,0.13,0.12; p1=3.2,3,3.5; γ1=0.83,0.85,0.86; σ1=0.62,0.61,0.63; a2=2,2.1,2.2;b2=0.11,0.13,0.12; p2=3.1,3.3,3.2;γ2=0.81,0.8,0.83; σ2=0.71,0.75,0.71;τ=1,2,3.In Figure4, we can choose a1=17,18,19;b1=0.14,0.35,0.15; p1=4.55,6,5.43;γ1=0.8,0.5,0.6; σ1=0.5,0.45,0.12;a2=16.5,16.9,17.5; b2=0.5,0.13,0.22;p2=3.23,4.67,6; γ2=0.3,0.4,0.8;σ2=0.49,1,0.65;τ=3,3,3. By calculation, we get τ1=τ2=τ3=3, 35=2β^1>σˇ12+bˇ1+1+2bˇ2+2+1/2=pˇ1+1+2pˇ2≈32.5350, and 33.28=2β^2>σˇ22+1+2bˇ1+bˇ2+1+2pˇ1+2+1/2pˇ2≈33.0729. Therefore, conditions of Theorem 5 have been checked. Thus, from Theorem 5, all species become extinct. Figure 4 clearly supports this result.Figure 4
Numerical solutions of (6). The parameters of model (5) are a1=17,18,19;b1=0.14,0.35,0.15; p1=4.55,6,5.43;γ1=0.8,0.5,0.6; σ1=0.5,0.45,0.12;a2=16.5,16.9,17.5; b2=0.5,0.13,0.22;p2=3.23,4.67,6; γ2=0.3,0.4,0.8;σ2=0.49,1,0.65; and τ=3,3,3.


## 7. Conclusions

By the conclusion of Lemma1, it is worthy to point out that the Brownian noise and colored noise will not destroy a great property that the solution of (6) may not explode. Especially, system (6) reduces to (1)–(4) or the model in [32] when parameters of system (6) take some special values. From Lemma 1, the condition α>σ2/2 in [21] is too strict and unnecessary. In Theorem 1, we comprehensively analyze ultimate boundedness in the θth moment and boundedness for the average in time of the θth moment of solution, which is the improvement of Theorem 3.1 in [21], Theorem 2.2 in [28], Theorem 3.3 in [31], and Theorem 3.2 in [32]. In Theorem 4, we find an upper bound Q/2 of the sample Lyapunov exponent. When parameters of system (6) take some special values, we compute that the upper bound Q/2 is less than the corresponding upper bound in [21, 28]. Furthermore, we find that the condition 2α1−σ12−b1+b2ϵ>0,2α2−σ22−b1+b2/ϵ>0 in [31] is not necessary. Despite all this, if we let parameter ϵ satisfy the above conditions, we compute that Q/2 is less than the upper bound in [31]. One point should be stressed is that the method for extinction in Theorem 5 can be used successfully for the models in [21, 28, 31, 32]. And then, Corollaries 1 and 2 give the conditions of extinction of (2) and (4), respectively. From Remarks 1–7, our work is a generalization and promotion of the corresponding work in [21, 28, 30–32]. To some extent, our proposed approaches are both more robust and more efficient than the existing methods.



---

*Source: 9078471-2020-05-23.xml*