# Dynamical Analysis of a Class of Prey-Predator Model with Beddington-DeAngelis Functional Response, Stochastic Perturbation, and Impulsive Toxicant Input

Authors: Feifei Bian; Wencai Zhao; Yi Song; Rong Yue
Journal: Complexity
 (2017)

---

## Abstract

A stochastic prey-predator system in a polluted environment with Beddington-DeAngelis functional response is proposed and analyzed. Firstly, for the system with white noise perturbation, by analyzing the limit system, the existence of boundary periodic solutions and positive periodic solutions is proved and the sufficient conditions for the existence of boundary periodic solutions and positive periodic solutions are derived. And then for the stochastic system, by introducing Markov regime switching, the sufficient conditions for extinction or persistence of such system are obtained. Furthermore, we proved that the system is ergodic and has a stationary distribution when the concentration of toxicant is a positive constant. Finally, two examples with numerical simulations are carried out in order to illustrate the theoretical results.


---

## Body

## 1. Introduction and Model Formulation

The Lotka-Volterra model [1–3] is a classical model in the study of biological mathematics, and the continuous Lotka-Volterra model which is modeled by ordinary differential equations and delay differential equations is widely used to characterize the dynamics of biological systems [4–13]. The functional response functions are important in the population ecological models [14]. In general, functional responses fall into two categories: one depends only on the density of the prey, such as Holling I–III [15–17]; the other depends on the density of both the prey and the predator, such as Beddington-DeAngelis type [18, 19]. Compared with the Holling II functional response, the Beddington-DeAngelis type functional response, F=f12xy/b2+y+W12x, has an additional term y in the denominator modeling mutual interference among predators. In other words, this type of functional response is affected by both predator and prey. Some biologists believe that if the predators compete with each other to obtain food, functional response should depend on the density of both the prey and the predator. Arditi et al. [20] and Jost et al. [21, 22] used the actual observation data to verify this point. In particular, having collected observation data from 19 predator-prey communities, Skalski and Gilliam [23] found that predator-dependent functional responses were in agreement with the observation data, and in many instances, the Beddington-DeAngelis type looked better than the others. The Beddington-DeAngelis functional response has been widely used in the modeling of ecosystems in which there is mutual interference among predators [24, 25]. In [19], DeAngelis et al. have extensively investigated the dynamical properties of the following prey-predator system:(1)x˙t=xta1-f12ytb2+yt+W12xt-d1-g1xt,y˙t=yte12f12xtb2+yt+W12xt-d2-g2yt,where xt and yt represent the density of the prey and the predator, respectively. a1 is the intrinsic growth rate of the prey, f12, b2, and W12 are the consumption rate, the saturation constant, and the saturation constant for an alternative prey, respectively. e12 is the conversion rate of nutrients into the reproduction for the predator. The parameters di and gii=1,2 are the nonpredatory loss rate and the interspecific competition rate. We refer the reader to [19] for more details.In many ecosystems, predators tend to be omnivorous, they have wide variety of food sources. For example, the giant panda is omnivorous animal, since it can eat both meat and plant such as bamboos. In the lake ecosystem, some fishes not only prey on aquatic invertebrates, but also feed on algae and other aquatic plants. Polis and Strong in [26] and McCann and Hastings in [27] studied omnivorous nature of animals in the food chain in 1996 and 1997, respectively. Based on the above literature, we established a kind of omnivorous model as follows:(2)x˙t=xtr1-b1xt-λyta+myt+nxt,y˙t=ytr2+γxta+myt+nxt-b2yt,where r2 represents the growth rate of y due to omnivorous nature and bii=1,2 denote the density-dependent coefficient of the prey and the predator, respectively. λ, a, m, n, and γ are the consumption rate, the saturation constant, the predator interference, the saturation constant for an alternative prey, and the conversion rate, respectively. All parameters are positive in system (2).It is well known that the biological population is inevitably affected by environment perturbation while the stochastic population model is more in line with the actual situation. Recently, various models based on stochastic differential equations (SDEs) have extensively been paid the attention of the researchers (see, e.g., [28–37]). Parameter perturbation induced by white noise is an important and common form to describe the effect of stochasticity (see, e.g., [37–48]). In this paper, we consider the white noise perturbation for the intrinsic growth rates of the prey and predator; that is, r1→r1+σ1B˙1t and r2→r2+σ2B˙2t, where B1t, B2t are mutually independent Brownian motions and σ1, σ2 denote the intensities of the white noise. On the other hand, it can be seen from the recent literature that the environmental pollution has an important effect on the population systems [49–60]. In 1983, Hallam et al. [61, 62] studied the influence of environmental pollution on the population and established a relationship model between environmental toxins and population. Subsequently, Hallam et al. [63, 64] studied the persistence and extinction of population in polluted environment. The mathematical model established by Hallam et al. considered only the toxins in the organism to cause a decrease in the birth rate or an individual death, which is reasonable in the case of lower concentration of the toxicant in the environment. When pollution is serious, the emission of pollutants may directly lead to the death of the species; see [65–69]. The authors in [68] added the environmental toxic term directly to the model; this is reasonable in the heavily polluted environment. For example, in a lake ecosystem, the discharge of large amounts of industrial waste water may directly lead to the death of fish, aquatic invertebrates, and so on. Therefore, we assume that the emission of pollutants to the environment is impulsive and directly affects the survival of the species in such an environment, so we get the following system:(3)dxt=xtr1-b1xt-λyta+myt+nxt-β1cetdt+σ1txtdB1t,dyt=ytr2+γxta+myt+nxt-b2yt-β2cetdt+σ2tytdB2t,dcetdt=-hcet,t≠kτ,Δxt=0,Δyt=0,Δcet=μ,t=kτ,where σ1t, σ2t are positive, nonconstant, and continuous functions of period τ, cet stands for the concentration of the toxicant in the environment, h denotes the loss rate of toxicant at time t, τ is the impulsive input period and μ is the impulsive input amount, and β1 and β2 represent the dose-response of the prey and predator to the environmental toxicant, respectively.Furthermore, the prey-predator model may be perturbed by telegraph noise which is distinguished by factors such as rain falls and nutrition and can be represented by switching among two or more regimes of environment [40, 60, 70–80]. For example, population growth rates in different seasons are not the same. The intraspecific competition coefficient varies according to the changes in nutrition and food resources. Generally, the switching between different regimes is memoryless and the waiting time for the next switch is exponentially distributed [81, 82]. Therefore, it can be described by a continuous-time Markov chain rt taking values in a finite state space S=1,2,…,m. Taking into account the influences of white noise and telegraph noise, we propose the following stochastic differential system with impulsive toxicant input:(4)dxt=xtr1rt-λrtytart+mrtyt+nrtxt-β1rtcet-b1rtxtdt+σ1rtxtdB1t,dyt=ytr2rt+γrtxtart+mrtyt+nrtxt-β2rtcet-b2rtytdt+σ2rtytdB2t,dcetdt=-hcet,t≠kτ,Δxt=0,Δyt=0,Δcet=μ,t=kτ.For any k∈S, rik, bik, βik, σiki=1,2, λk, ak, mk, nk, and γk are all positive constants. In model (4), the population is inevitably affected by severe stochastic interference such as drought; the parameter switches one state rt=i into another state rt=j and it will switch into the next regime until the next major environmental change.The rest of this paper is organized as follows. In Section2, we provide preliminaries which are used in the following sections. In Section 3, we show that system (3) admits a nontrivial positive τ-periodic solution by constructing Lyapunov function. In Section 4, we explore the sufficient conditions for extinction and permanence in mean of system (4). Finally, some examples with numerical simulations have been given to illustrate our theoretical results.


## 2. Preliminaries

Throughout this paper, letΩ,F,Ft≥0,P be a complete probability space with a filtration Ft≥0 satisfying the usual conditions, Biti=1,2 is one-dimensional Brownian motion on this space, and rt is a right-continuous Markov chain and independent of the Brownian motion Bit. The state space of this Markov chain is S=1,2,…,m. Suppose that the generator matrix of rt is Γ=qij1≤i,j≤m, where qij stands for the transition rate from state i to j and satisfies the following conditions: (5)Prt+Δt=j∣rt=i=qijΔt+oΔt,ifi≠j,1+qijΔt+oΔt,ifi=j;here, qij≥0 if i≠j, while qii=-∑i≠jqij, i,j=1,2,…,m. As a standing hypothesis, we assume that the Markov chain rt is irreducible, which means that system (4) can switch from one regime to another. Under this assumption, the Markov chain has a unique stationary distribution π=π1,π2,…,πm which is the solution of the system of linear equations πΓ=0 subject to ∑j=1mπj=1 and πj>0 for all j∈S. Hence, for any vector ϖ=ϖ1,…,ϖmT, we have that (6)limt→∞⁡1t∫0tϖrsds=∑k∈Sπkϖk.Let us consider the following stochastic differential equation with Markov conversion.(7)dxt=fxt,rtdt+gxt,rtdBt,x0=x0,r0=r0,where xt=x1t,…,xntT∈Rn, f:Rn×S→Rn, g:Rn×S→Rn×d, and Bt is a d-dimensional Brownian motion defined on the underlying probability space. The n×n matrix (8)Gx,k=gx,kgTx,k=Gijn×nis called the diffusion matrix. Let V:Rn×S→Rn be twice continuously differentiable and LVx,k which is defined as follows be the diffusion operator about Vx,k: (9)LVx,k=∑i=1nfix,k∂Vx,k∂xi+12∑i,j=1nGij∂2Vx,k∂xi∂xj+∑i≠k∈SqkiVx,i-Vx,k.Particularly, for one-dimensional stochastic system(10)dxt=xtart-brtxtdt+αrtxtdBt,x0=x0,r0=r0;the following two lemmas can be given from referring to the articles [72, 77].Lemma 1.
System (10) has a unique continuous positive solution xt. When it exists, the solution is global and stochastically ultimately bounded.Lemma 2.
Suppose thatχ=∑k∈Sπkak-1/2α2k≠0; then
(i) system (10) is stochastic permanent if and only if χ>0;
(ii) system (10) is extinct if and only if χ<0;
(iii) whenχ>0, system (10) is ergodic and there exists a unique stationary distribution v·,·, such that (11)χ=∑k∈Sbk∫R+xvdx,k.Next, we consider the following stochastic differential equation:(12)dyt=ft,ytdt+gt,ytdBt.Lemma 3 (see [78]).
Suppose that the coefficients of (12) are τ-periodic in t and there exists a function Vt,y∈C2 which is τ-periodic in t, and Vt,y satisfies the following conditions:
(i)infy>R⁡Vt,y→∞ as R→∞.
(ii)LVt,y≤-1 outside some compact set.
Then there exists a solution for (12) which is a τ-periodic Markov process.Furthermore, we introduce some results from [80, 83] in Lemmas 4 and 5, which will be used in next section.Lemma 4 (see [80]).
LetXt∈CΩ×0,+∞,R+. Then
(i) if there are two positive constantsT and m0 such that (13)ln⁡Xt≤mt-m0∫0tXsds+∑k=1nωkBktholds for all t≥T and constants ωkk=1,2,…,n, then (14)limsupt→+∞1t∫0tXsds≤mm0,a.s.ifm>0,limt→+∞⁡Xt=0,a.s.ifm<0,
(ii) if there are three positive constantsT, m, and m0 such that (15)ln⁡Xt≥mt-m0∫0tXsds+∑k=1nωkBktholds for any t≥T, then (16)liminft→+∞1t∫0tXsds≥mm0a.s.Finally, we give some basic properties of the following subsystem of system (3),(17)dcetdt=-hcet,t≠kτ,k∈Z,Δcet=μ,t=kτ,k∈Z.Lemma 5 (see [83]).
System (17) has a unique τ-periodic solution ce∗t which is globally asymptotically stable. Here ce∗t=μe-ht-kτ/1-e-hτ, t∈kτ,k+1τ, cmax=μ/1-e-hτ, and cmin=μe-hτ/1-e-hτ.For convenience and simplicity, defineα^=mini∈S⁡αi, αˇ=maxi∈S⁡αi, and φθ=1/θ∫0θφsds, where φt is an integrable function on 0,+∞. If f is a bounded function on 0,+∞, define fu=supt∈0,+∞⁡ft.


## 3. Existence of Periodic Solutions of System (3)

In this section, we devote our attention to the investigation of the existence of periodic solutions of system (3). From Lemma 5, we know that system (17) has a globally asymptotically stable periodic solution ce∗t; therefore, the limit system of (3) is(18)dxt=xtr1-b1xt-λyta+myt+nxt-β1ce∗tdt+σ1txtdB1t,dyt=ytr2+γxta+myt+nxt-b2yt-β2ce∗tdt+σ2tytdB2t,x0=x0,y0=y0,where σ1t, σ2t, and ce∗t are all positive and continuous functions of period τ.Now, we discuss the existence of periodic solutions of system (18).Define(19)hi=1τ∫0τri-βice∗s-12σi2sds,i=1,2.Then, we have the following theorem about periodic solutions of system (18).Theorem 6.
Ifh1<0 and h2>0, there exists a prey extinction periodic solution 0,y∗t of system (18).Proof.
From the first equation of system (18), it is easy to see (20)dxt≤xtr1-b1xt-β1ce∗tdt+σ1txtdB1t.Applying Itô’s formula and then integrating from 0 to t, we obtain(21)ln⁡xt-ln⁡x0≤∫0tr1-β1ce∗s-12σ12sds-∫0tb1xsds+Mt,where Mt=∫0tσ1sdB1s is local martingale. From strong law of large numbers for martingales (see [84]), we have (22)limt→∞⁡Mtt=0,a.s.It then follows from (21) by dividing t on both sides and letting t→∞ that (23)limsupt→∞ln⁡xtt≤r1-β1ce∗s-12σ12sτ=r1-μβ1τh-12σ12sτ=h1<0;namely, xt tends to zero exponentially almost surely.
Sincelimt→∞⁡xt=0, a.s., from the second equation of system (18), its limit system is(24)dyt=ytr2-β2ce∗t-b2yt+σ2tytdB2t.According to Theorem  4.2 in [85], when r2-β2ce∗t>0 and (25)r2-β2ce∗s-12σ22sτ=r2-μβ2τh-12σ22sτ=h2>0,(24) has a unique positive τ-periodic solution y∗t.
Overall, whenh1<0 and h2>0, there exists a prey extinction periodic solution 0,y∗t of system (18).
The proof of this theorem is completed.In order to investigate the existence of a nontrivial positiveτ-periodic solution for system (18), first of all, we assume following conditions hold.(H1)hi>0, i=1,2.(H2)λ1=r1-β1cmin+r1-β1cmax2-4r1-β1cmaxh1>0 and ξ1>γ2/4b1h2λ1.(H3)ξ2=ab2h1-λr2+γ/n-β2cmin>0 and ξ2/ab1r1-β1cmax>b2ξ1/γ2-h2+r2+γ/n-β2cmin.Theorem 7.
Suppose thatH1, H2, and H3 hold, then there exists a positive τ-periodic solution for system (18).Proof.
Obviously, the coefficients of system (18) are continuous bounded positive periodic functions in t. Now, we show that conditions (i) and (ii) of Lemma 3 hold. Define a nonnegative C2-function (26)Vt,x,y=x-r1-β1cmaxb1ln⁡x+r1-β1cmaxb1ω1t-ξ1γ2ln⁡y+ξ1γ2ω2t+qy≕V1+V2+V3,where V1=x-r1-β1cmax/b1ln⁡x+r1-β1cmax/b1ω1t, V2=-ξ1/γ2ln⁡y+ξ1/γ2ω2t, V3=qy, q=b2r1-β1cmax/b1h1+ξ1/γ2h2/r2-β2cmin+γ/n2, and ωit is a function defined on 0,∞ satisfying ωi′t=ri-1/2σi2t-βice∗t-hi and ωi0=0i=1,2. Obviously, ωit is a τ-periodic function on 0,∞. Therefore, the function Vt,x,y is τ-periodic in t and satisfies (27)liminfk→∞,x,y∈R+2∖Uk⁡Vt,x,y=∞,where Uk=1/k,k×1/k,k. Therefore, condition (i) of Lemma 3 holds. Next, we will prove that condition (ii) of Lemma 3 also holds.
Applying Itô’s formula, one has(28)LV1=xr1-b1x-λya+my+nx-β1ce∗t-r1-β1cmaxb1r1-b1x-λya+my+nx-β1ce∗t-σ12t2+r1-β1cmaxb1ω1′t≤-b1x2+r1-β1cmin+r1-β1cmaxx+λb1ar1-β1cmaxy-r1-β1cmaxb1h1,LV2=-ξ1γ2r2+γxa+my+nx-b2y-β2ce∗t-σ22t2+ξ1γ2ω2′t≤ξ1b2γ2y-ξ1γ2h2,LV3=qyr2+γxa+my+nx-b2y-β2ce∗t≤qr2+γn-β2cminy-qb2y2.Therefore,(29)LV≤-b1x2+r1-β1cmin+r1-β1cmaxx-r1-β1cmaxb1h1-qb2y2+λr1-β1cmaxb1a+qr2+γn-β2cmin+ξ1b2γ2y-ξ1γ2h2.Define a bounded closed set (30)D=x,y∈R+2:ε≤x≤1ε,ε≤y≤1ε,where 0<ε<1 is a sufficiently small number such that(31)r1-β1cmin+r1-β1cmaxε≤12r1-β1cmaxb1h1+12ξ1γ2h2-12λ/b1ar1-β1cmax+ξ1b2/γ2+qr2-β2cmin+γ/n24qb2,(32)λr1-β1cmaxb1a+qr2+γn-β2cmin+ξ1b2γ2ε<-18b1r1-β1cmin+r1-β1cmax2-4r1-β1cmaxh1-4b1h2γ2ξ1,(33)-b12ε2+K3≤-1,(34)-b2q2ε2+K4≤-1,and K3, K4 are quantities to be determined in the rest of the proof.
Denote(35)Dε1=x,y∈R+2:0<x<ε,Dε2=x,y∈R+2:0<y<ε,Dε3=x,y∈R+2:x>1ε,Dε4=x,y∈R+2:y>1ε.Note that R+2∖D=Dε1∪Dε2∪Dε3∪Dε4. Now, we prove LVt,x,y≤-1,x,y∈R+2∖D.
Case 1. If (x,y)∈Dε1, from (29), it implies that (36)LV≤-qb2y-λ/b1ar1-β1cmax+ξ1b2/γ2+qr2-β2cmin+γ/n2qb22+λ/b1ar1-β1cmax+ξ1b2/γ2+qr2-β2cmin+γ/n24qb2-r1-β1cmaxb1h1-ξ1γ2h2+r1-β1cmin+r1-β1cmaxε≤K1,where K1=1/2{λ/b1a(r1-β1cmax)+ξ1b2/γ2+q(r2-β2cmin+γ/n)2/4qb2-r1-β1cmax/b1h1-ξ1/γ2h2}<0. In fact, from condition (H3), one can get (37)λb1ar1-β1cmax+ξ1b2γ2+qr2-β2cmin+γn-2qb2r1-β1cmaxb1h1+ξ1γ2h2=ξ1b2/γ2-h2+r2+γ/n-β2cmin-ξ2/ab1r1-β1cmaxr2+γ/n-β2cmin<0;that is to say, K1<0.
Case 2. If (x,y)∈Dε2, from (29) and (32), we can get (38)LV≤-b1x-r1-β1cmin+r1-β1cmax2b12+14b1r1-β1cmin+r1-β1cmax2-4r1-β1cmaxh1-ξ1γ2h2+λr1-β1cmaxb1a+qr2+γn-β2cmin+ξ1b2γ2ε≤14b1r1-β1cmin+r1-β1cmax2-4r1-β1cmaxh1-4b1ξ1γ2h2+λr1-β1cmaxb1a+qr2+γn-β2cmin+ξ1b2γ2ε≤K2,where K2=1/8b1(r1-β1cmin+r1-β1cmax)2-4(r1-β1cmax)h1-4b1h2/γ2ξ1. Using condition (H2), one can get K2<0.
Case 3. If (x,y)∈Dε3, then (39)LV≤-b12x2-b12x2+r1-β1cmin+r1-β1cmaxx-r1-β1cmaxb1h1-b2qy2+λr1-β1cmaxab1+ξ1b2γ2+qr2-β2cmin+γny-ξ1γ2h2≤-b12ε2+K3,where (40)K3=supx,y∈R+2⁡-b12x2+r1-β1cmin+r1-β1cmaxx-r1-β1cmaxb1h1-b2qy2+λr1-β1cmaxab1+ξ1b2γ2+qr2-β2cmin+γny-ξ1γ2h2.By (33), we have LV≤-1.
Case 4. If (x,y)∈Dε4, then (41)LV≤-b2q2y2-b1x2+r1-β1cmin+r1-β1cmaxx-r1-β1cmaxb1h1-b2q2y2+λr1-β1cmaxab1+ξ1b2γ2+qr2-β2cmin+γny-ξ1γ2h2≤-qb22ε2+K4,where (42)K4=supx,y∈R+2⁡-b1x2+r1-β1cmin+r1-β1cmaxx-r1-β1cmaxb1h1-qb22y2+λr1-β1cmaxab1+ξ1b2γ2+qr2-β2cmin+γny-ξ1γ2h2.By (34), we obtain LV≤-1.
Thus,(43)LV≤min⁡K1,K2,-1,x,y∈R+2∖D.Therefore, the proof of Theorem 7 is completed.


## 4. Extinction and Persistence in Mean of System (4)

In this section, we investigate the long-term dynamic behaviors of the prey-predator system (4) with white noise and telegraph noise in a polluted environment and then discuss the extinction and average persistence of prey and predator. According to Lemma 5, the periodic solution ce∗(t) of the toxicant input is globally asymptotically stable, so the limit system of (4) is(44)dxt=xtr1rt-b1rtxt-λrtytart+mrtyt+nrtxt-β1rtce∗tdt+σ1rtxtdB1t,dyt=ytr2rt+γrtxtart+mrtyt+nrtxt-b2rtyt-β2rtce∗tdt+σ2rtytdB2t,x0=x0,y0=y0,r0=r0.In order to obtain the threshold conditions of persistence and extinction of system (44), we assume thatA1:η1=∑k=1mπkr1k-1/2σ12k-cminβ1k<0,(A2):η2=∑k=1mπkr2k+γ(k)/n(k)-1/2σ22k-cminβ2k<0,(A3):η3=∑k=1mπkr2k-1/2σ22k-cminβ2k<0,(A4):η4=∑k=1mπkr1k-λ(k)/m(k)-1/2σ12k-cmaxβ1k>0,(A5):η5=∑k=1mπkr1k-1/2σ12k-cmaxβ1k>0.Theorem 8.
Given initial value(x(0),y(0),r(0))∈R+2×S for system (44), then
(i) if(A1) is established, the prey population will be extinct,
(ii) if(A2) is established, the predator population will be extinct,
(iii) if(A1) and (A3) are established, both the prey and the predator will die out.Proof.
(i) By Itô’s formula, we get(45)dln⁡xt=r1rt-b1rtxt-λrtytart+mrtyt+nrtxt-σ12rt2-β1rtce∗tdt+σ1rtdB1t≤r1rt-σ12rt2-β1rtcmindt+σ1rtdB1t;then, (46)ln⁡xt-ln⁡x0t≤1t∫0tr1rt-σ12rt2-β1rtcmindt+1t∫0tσ1rtdB1t.By the ergodic theory of the Markov chain and the strong law of large number, we have (47)limsupt→∞⁡ln⁡xtt≤∑k=1mπkr1k-12σ12k-cminβ1k=η1,a.s.;from (A1), we know (48)limt→∞⁡xt=0.a.s.
(ii) Similarly, from the second equation of system (44), we have (49)dln⁡yt=r2rt-b2rtyt+γrtxtart+mrtyt+nrtxt-σ22rt2-β2rtce∗tdt+σ2rtdB2t;then, (50)ln⁡yt-ln⁡y0t≤1t∫0tr2rt+γrtnrt-σ22rt2-β2rtcmindt+1t∫0tσ2rtdB2t;further, (51)limsupt→∞⁡ln⁡ytt≤∑k=1mπkr2k+γknk-12σ22k-cminβ2k=η2,a.s.;from (A2), we know (52)limt→∞⁡yt=0.a.s.
(iii) By the condition(A1), one can get limt→∞x(t)=0 a.s., so that the limit system of the second equation of system (44) is (53)dyt=ytr2rt-b2rtyt-β2rtce∗tdt+σ2rtytdB2t≤ytr2rt-b2rtyt-β2rtcmindt+σ2rtytdB2t.Through Lemma 2, if (A3) holds, we obtain (54)limt→∞⁡yt=0.a.s.
This completes the proof of the theorem.Remark 9.
Ifri(k)(i=1,2) remains unchanged and σi(k)(i=1,2) or cmin increases so that η1<0 or η2<0, then condition (A1) or (A2) is established. That is to say, if the intrinsic growth rate and the predation intensity are relatively fixation, the increase of white noise intensity or pollutant concentration will lead to the extinction of the biological population.Next, we will discuss the persistence of system (44). Applying Itô’s formula to the first equation of system (44), one can get (55)dln⁡xt≥r1rt-b1rtxt-λrtmrt-σ12rt2-β1rtcmaxdt+σ1rtdB1t;then, (56)ln⁡xt-ln⁡x0t≥1t∫0tr1rt-b1rtxt-λrtmrt-σ12rt2-β1rtcmaxdt+1t∫0tσ1rtdB1t,when t is large enough, we derive (57)ln⁡xtt≥∑k=1mπkr1k-λkmk-12σ12k-cmaxβ1k-bˇ1t∫0txtdt-ε,where ε is a sufficiently small positive number. In view of (A4) and Lemma 4, we deduce (58)liminft→∞⁡1t∫0txtdt≥η4bˇ1,a.s.That is, the prey population of system (44) will be persistence in mean under condition (A4).Furthermore, the persistent property of the predator species of system (44) can be investigated as follows.From the first equation of system (44), we have (59)dxt≤xtr1rt-b1rtxt-β1rtce∗tdt+σ1rtxtdB1t.Consider the following stochastic differential equations:(60)dϕt=ϕtr1rt-b1rtϕt-β1rtce∗tdt+σ1rtϕtdB1t,ϕ0=x0,(61)dψt=ψtr1rt-b1rtψt-β1rtcmaxdt+σ1rtψtdB1t,ψ0=x0.Obviously, x(t)≤ϕ(t), ψ(t)≤ϕ(t). Using Lemma 2, if η5>0, system (61) is ergodic and there exists a unique stationary distribution μψ, such that (62)∑k=1mb1k∫R+xμψdx,k=∑k=1mπkr1k-12σ12k-cmaxβ1k.Applying Itô’s formula to (60) and then integrating from 0 to t, we get (63)ln⁡ϕt-ln⁡ϕ0t=1t∫0tr1rt-σ12rt2-β1rtce∗t-b1rtϕtdt+1t∫0tσ1rtdB1t.From the first equation of system (44), it yields that (64)ln⁡xt-ln⁡x0t=1t∫0tr1rt-σ12rt2-β1rtce∗tdt-1t∫0tλrtytart+mrtyt+nrtxtdt-1t∫0tb1rtxtdt+1t∫0tσ1rtdB1t.Due to x(t)≤ϕ(t), one can get(65)1t∫0tϕt-xtdt≤1t∫0tλˇa^b^1ytdt.a.s.From the second equation of system (44) we have (66)ln⁡ytt=ln⁡y0t+1t∫0tr2rt-σ22rt2-β2rtce∗tdt-1t∫0tb2rtytdt+1t∫0tγrtϕtart+nrtϕtdt-1t∫0tγrtϕtart+nrtϕt-γrtxtart+nrtxtdt-1t∫0tγrtxtart+nrtxt-γrtxtart+mrtyt+nrtxtdt+1t∫0tσ2rtdB2t≥ln⁡y0t+1t∫0tr2rt-σ22rt2-β2rtcmaxdt+1t∫0tγrtψtart+nrtψtdt-bˇ2+γˇλˇa^2b^1+γˇmˇa^n^1t∫0tytdt+1t∫0tσ2rtdB2t.By Lemma 4, if (A6):η6=∑k=1mπkr2k-1/2σ22k-cmaxβ2k+∑k=1mπk∫R+γkx/ak+nkxμψdx,k>0holds, then(67)liminft→∞1t∫0tytdt≥η6bˇ2+γˇλˇ/a^2b^1+γˇmˇ/a^n^a.s.In summary, one gets the following.Theorem 10.
Given initial value(x(0),y(0),r(0))∈R+2×S for system (44), then
(i) if(A4) is established, the prey population will be persistent in mean,
(ii) if conditions(A5) and (A6) are satisfied, the predator population will be persistent in mean.Remark 11.
(i) It can be seen fromη4, in the case where the intrinsic growth rate and the predation intensity are relatively constant, only by reducing the intensity of white noise or pollutant concentration, so that (A4) can be established to ensure the lasting survival of the prey population.
(ii) Obviously,η4<η5; if the prey population is persistent, the predator population is persistent as long as the white noise interference intensity or the toxin concentration is small enough, such that (A6) is established. As can be seen from condition (A6), the omnivorous nature of y contributes to its permanence.In system (4), if the concentration of the toxicant in the environment remains unchanged, that is, ce(t)=c is a positive constant, then the system can be converted into the following system:(68)dxt=xtr1rt-b1rtxt-λrtytart+mrtyt+nrtxt-cβ1rtdt+σ1rtxtdB1t,dyt=ytr2rt+γrtxtart+mrtyt+nrtxt-b2rtyt-cβ2rtdt+σ2rtytdB2t,x0=x0,y0=y0,r0=r0.Lemma 12 (see [77]).
System (7) is ergodic and positive recurrent if the following conditions are satisfied:
(i) Fori≠j, qij>0.
(ii) For eachk∈S,(69)σξ2≤ξTGx,kξ≤σ-1ξ2for all ξ∈Rn, with some constant σ∈(0,1] for all x∈Rn.
(iii) There exists a bounded open setD⊂Rn with a smooth boundary satisfying that, for each k∈S, there is a twice continuously differentiable nonnegative function V(·,·):Dc→R and that for some ζ>0, LV(x,k)≤-ζ, for any (x,k)∈Dc×S.
Moreover, the Markov process(x(t),r(t)) has a unique ergodic stationary distribution μ(·,·). Hence, for any Borel measurable function H(·,·):Rn×S→R, if ∑k∈S∫RnH(x,k)μ(x,k)dx<∞, then (70)Plimt→∞⁡1t∫0tHxs,rsds=∑k∈S∫RnHx,kμx,kdx=1.Recently, the ergodicity and stationary distribution have been explored by many authors. In the following, we give sufficient conditions for the existence of stationary distribution of system (68) and prove the following theorem by showing that system (68) satisfies the three conditions in Lemma 12.Theorem 13.
Assume that fori≠j, qij>0 and(B1):∑k∈Sπkr1k-λ(k)/m(k)-σ12k/2-cβ1k>0,(B2):∑k∈Sπkr2k-σ22k/2-cβ2k>0
hold; then the stochastic process(x(t),y(t),r(t)) of system (68) is ergodic and has a unique stationary distribution in R+2×S.Proof.
By the assumptionqij>0 for i≠j in Theorem 13, condition (i) in Lemma 12 is satisfied. Let gx,y,k=diag⁡(σ1(k)x,σ2(k)y); then (71)Gx,y,k=gx,y,kgTx,y,k=diag⁡σ12kx2,σ22ky2.Define a bounded open subset as follows: (72)U=ε,1ε×ε,1ε⊂R+2,where 0<ε<1 is a constant. Let σ=min⁡{1,M1,1/M2}; here M1=min(x,y,k)∈U¯×S{σ12(k)x2+σ22(k)y2} and M2=max(x,y,k)∈U¯×S{σ12(k)x2+σ22(k)y2}. For (x,y,k)∈U×S, we have (73)σξ2≤ξTGx,y,kξ=σ12kx2ξ12+σ22ky2ξ22≤σ-1ξ2for all ξ∈R2. Thus condition (ii) in Lemma 12 holds. Therefore, it remains for us to verify condition (iii) in Lemma 12.
Define aC2-function on R+2×S, (74)Vx,y,k=1-υηkx-υ+1-υδky-υ+x+2y+M≕V4+V5,where V4=(1-υηk)x-υ+(1-υδk)y-υ, V5=(x+2)(y+M), υ is a sufficiently small positive constant satisfying υ<min⁡1/ηk,1/δk, and M≥1+1/4b^1b^2rˇ1-cβ^1+rˇ2-cβ^2+γˇ/n^2. ηk, δk are quantities to be determined below.
An application of the operatorL to V4 yields(75)LV4x,y,k=-υ1-υηkx-υr1k-b1kx-λkyak+mky+nkx-cβ1k+12υ1+υ1-υηkσ12kx-υ-υ∑i≠kqkiηi-ηkx-υ-υ1-υδky-υr2k+γkxak+mky+nkx-b2ky-cβ2k+12υ1+υ1-υδkσ22ky-υ-υ∑i≠kqkiδi-δky-υ≤υ1-υηk-r1k-λkmk-cβ1k-σ12k2-∑i≠kqkiηi-ηk+12υσ12k-υηk1-υηk∑i≠kqkiηi-ηkx-υ+υ1-υδk-r2k-cβ2k-12σ22k-∑i≠kqkiδi-δk+12υσ22k-υδk1-υδk∑i≠kqkiδi-δky-υ+υ1-υηkb1kx1-υ+υ1-υδkb2ky1-υ.Define the vectors ξ=(ξ1,ξ2,…,ξm)T and ϑ=(ϑ1,ϑ2,…,ϑm)T with ξk=-r1(k)-λ(k)/m(k)-cβ1(k)-σ12/2(k), ϑk=-r2(k)-cβ2(k)-1/2σ22(k). As the generator matrix Γ is irreducible, for each ξk and ϑk, there exists η=(η1,η2,…,ηm)T and δ=(δ1,δ2,…,δm)T, respectively, which is a solution of the Poisson system [78] (76)Γηk-ξk=-∑j=1mπjξj,Γδk-ϑk=-∑j=1mπjϑj.Therefore we have(77)-r1k-λkmk-cβ1k-σ122k-∑i≠kqkiηi-ηk=-∑k∈Sπkr1k-λkmk-σ12k2-cβ1k,(78)-r2k-cβ2k-12σ22k-∑i≠kqkiδi-δk=-∑k∈Sπkr2k-σ22k2-cβ2k.Combining (75), (77), and (78), we obtain (79)LV4x,y,k≤υ1-υηk-∑k∈Sπkr1k-λkmk-σ12k2-cβ1k+12υσ12k-υηk1-υηk∑i≠kqkiηi-ηkx-υ+υ1-υδk-∑k∈Sπkr2k-σ22k2-cβ2k+12υσ22k-υδk1-υδk∑i≠kqkiδi-δky-υ+υ1-υηkb1kx1-υ+υ1-υδkb2ky1-υ.Similarly, for V5(x,y), we calculate (80)LV5=y+Mxr1k-b1kx-λkyak+mky+nkx-cβ1k+x+2yr2k+γkxak+mky+nkx-b2ky-cβ2k≤-b^1x2+Mrˇ1-cβ^1x-b^2y2+2rˇ2-cβ^2+γˇn^y.From conditions (B1) and (B2), we can choose υ sufficiently small such that (81)-∑k∈Sπkr1k-λkmk-σ12k2-cβ1k+12υσ12k-υηk1-υηk∑i≠kqkiηi-ηk<0,-∑k∈Sπkr2k-σ22k2-cβ2k+12υσ22k-υδk1-υδk∑i≠kqkiδi-δk<0.So LV=LV4+LV5 can be estimated as follows: (82)LV≤υ1-υηk-∑k∈Sπkr1k-λkmk-σ12k2-cβ1k+12υσ12k-υηk1-υηk∑i≠kqkiηi-ηkx-υ+υ1-υηkb1kx1-υ-b^1x2+Mrˇ1-cβ1^x+υ1-υδk-∑k∈Sπkr2k-σ22k2-cβ2k+12υσ22k-υδk1-υδk∑i≠kqkiδi-δky-υ+υ1-υδkb2ky1-υ-b^2y2+2rˇ2-cβ^2+γˇn^y=f1x+f2y.
It is easy to see that(83)LV≤f1x+f2y≤f1x+f2uy⟶-∞,ifx⟶0orx⟶+∞,f1ux+f2y⟶-∞,ify⟶0ory⟶+∞.Consequently, we derive that, for a sufficiently small ε, (84)LVx,y,k≤-1,∀x,y,k∈Uc×S.Using Lemma 12, we obtain the conclusion that (x(t),y(t),r(t)) is ergodic and positive recurrent; that is, system (68) is positive recurrent and has a unique stationary distribution.
This completes the proof of Theorem13.


## 5. Conclusions and Numerical Simulations

In this article, we discussed the dynamics of stochastic prey-predator models with Beddington-DeAngelis functional response in polluted environment.Firstly, for system (3), there are the following properties:(
1
) If h1<0 and h2>0, the limit system of (3) has a prey extinction periodic solution (0,y∗(t)).(
2
) If conditions (H1), (H2), and (H3) are established, the limit system of (3) has a positive periodic solution.Secondly, system (4) possesses the following properties:(
1
) If η1<0, the prey population x(t) will be extinct.(
2
) If η2<0, the predator population y(t) will be extinct.(
3
) If η1<0 and η3<0, then prey population and predator population will die out.(
4
) If η4>0, the prey population x(t) will be persistent in mean.(
5
) If η5>0 and η6>0, the predator population y(t) will be persistent in mean.To verify the correctness of the theoretical analysis, numerical simulations are employed in the following examples.Assume that the Markov chainr(t) take values in S={1,2} with the generator (85)Γ=-112-2.By the linear equations πΓ=0, we can see (π1,π2)=2/3,1/3 which is the stationary distribution of r(t). Furthermore, in the following examples, we suppose τ=0.5 and h=0.1, consistently.


### 5.1. The Existence of Periodic Solutions of System (3)

Example 14.
Assumer1=1.2, r2=2.3, β1=0.18, β2=0.22, b1=0.5, b2=0.8, λ=0.3, γ=0.5, a=1.4, m=0.8, and n=1.8.Case 1.
We choose the density of white noise as the following:σ1(t)=0.09+0.8sin⁡2π/τt, σ2(t)=0.1+0.01sin⁡2π/τt, and let μ=0.45.
Note thath1=-0.5841<0, h2=0.3150>0. The conditions of Theorem 6 hold, so there exists a boundary periodic solution (0,y∗(t)) of system (3) (see Figure 1).Figure 1
Sample paths ofxt and yt with initial conditions x0=2.5 and y0=0.8.Case 2.
We change the density of the white noise toσ1(t)=0.02+0.2sin⁡2π/τt, σ2(t)=0.07+0.1sin⁡2π/τt, and μ=0.13. This gives h1=0.7218>0, h2=1.7230>0, λ1=0.0634>0, and ξ2=0.2024>0; choose ξ1=0.1 according to ξ1>γ2/4b1h1λ1 and b2ξ1/γ2-h2+r2+γ/n-β2cmin<ξ2/ab1(r1-β1cmax).
From Theorem7, we know that there exists a positive τ-periodic solution of system (3) (see Figure 2).Figure 2
Sample paths ofxt and yt with initial conditions x0=0.8 and y0=2.5.


### 5.2. The Extinction and Persistence of System (4)

Example 15.
Choose parameters(r1(1),r2(1))=(2.3,1.5), λ(1)=0.06, γ(1)=0.3, (b1(1),b2(1))=(1.8,1.6), a(1)=0.4, m(1)=0.6, n(1)=0.3, and (β1(1),β2(1))=(0.4,0.2), if k=1, and (r1(2),r2(2))=(2.5,1.3), λ(2)=0.08, γ(2)=0.2, (b1(2),b2(2))=(1.8,1.6), a(2)=0.4, m(2)=0.8, n(2)=0.4, and (β1(2),β2(2))=(0.5,0.3), if k=2.Case 1.
Letμ=0.15, (σ1(1),σ2(1))=(1.7,2.1), and (σ1(2),σ2(2))=(1.4,0.9); we note η1=-0.1911<0, η2=-0.0210<0, and η3=-0.8543<0. The conditions of Theorem 8 are satisfied, so the prey and predator are both extinct (see Figure 3).
Next we only change the density of the white noise to(σ1(1),σ2(1))=(0.3,0.2) and (σ1(2),σ2(2))=(0.5,0.4) and keep μ=0.15. Simple calculation shows that η4=0.8622>0, η5=0.9622>0, and η6>0.6757>0. The conditions of Theorem 10 are satisfied, so the prey and predator are persistent (see Figure 4).
It is easy to see from Figures3 and 4 that the increase of the intensity of white noise can result in the extinction of prey and predator.Figure 3
Sample paths ofxt, yt with initial value x0=2.8 and y0=2.2.
(a)
k
=
1
(b)
k
=
2Figure 4
Sample paths ofxt, yt with initial value x0=0.4 and y0=0.3.
(a)
k
=
1
(b)
k
=
2Case 2.
Letμ=0.55, (σ1(1),σ2(1))=(0.6,0.5), and (σ1(2),σ2(2))=(0.45,0.55), which gives η1=-2.4356<0, η2=-0.3701<0, and η3=-1.2035<0. x(t) and y(t) are extinct (see Figure 5).
Next we only change the amount of toxicant toμ=0.12 and keep (σ1(1),σ2(1))=(0.6,0.5) and (σ1(2),σ2(2))=(0.45,0.55), We note η4=1.0467>0, η5=1.1467>0, and η6>0.7255>0. Thus x(t) and y(t) are persistent in mean (see Figure 6).
Figures5 and 6 show that the increase of the amount of toxicant can also result in the extinction of the prey and predator.Figure 5
Sample paths ofxt, yt with initial value x0=6 and y0=4.
(a)
k
=
1
(b)
k
=
2Figure 6
Sample paths ofxt, yt with initial value x0=1.6 and y0=1.4.
(a)
k
=
1
(b)
k
=
2


## 5.1. The Existence of Periodic Solutions of System (3)

Example 14.
Assumer1=1.2, r2=2.3, β1=0.18, β2=0.22, b1=0.5, b2=0.8, λ=0.3, γ=0.5, a=1.4, m=0.8, and n=1.8.Case 1.
We choose the density of white noise as the following:σ1(t)=0.09+0.8sin⁡2π/τt, σ2(t)=0.1+0.01sin⁡2π/τt, and let μ=0.45.
Note thath1=-0.5841<0, h2=0.3150>0. The conditions of Theorem 6 hold, so there exists a boundary periodic solution (0,y∗(t)) of system (3) (see Figure 1).Figure 1
Sample paths ofxt and yt with initial conditions x0=2.5 and y0=0.8.Case 2.
We change the density of the white noise toσ1(t)=0.02+0.2sin⁡2π/τt, σ2(t)=0.07+0.1sin⁡2π/τt, and μ=0.13. This gives h1=0.7218>0, h2=1.7230>0, λ1=0.0634>0, and ξ2=0.2024>0; choose ξ1=0.1 according to ξ1>γ2/4b1h1λ1 and b2ξ1/γ2-h2+r2+γ/n-β2cmin<ξ2/ab1(r1-β1cmax).
From Theorem7, we know that there exists a positive τ-periodic solution of system (3) (see Figure 2).Figure 2
Sample paths ofxt and yt with initial conditions x0=0.8 and y0=2.5.


## 5.2. The Extinction and Persistence of System (4)

Example 15.
Choose parameters(r1(1),r2(1))=(2.3,1.5), λ(1)=0.06, γ(1)=0.3, (b1(1),b2(1))=(1.8,1.6), a(1)=0.4, m(1)=0.6, n(1)=0.3, and (β1(1),β2(1))=(0.4,0.2), if k=1, and (r1(2),r2(2))=(2.5,1.3), λ(2)=0.08, γ(2)=0.2, (b1(2),b2(2))=(1.8,1.6), a(2)=0.4, m(2)=0.8, n(2)=0.4, and (β1(2),β2(2))=(0.5,0.3), if k=2.Case 1.
Letμ=0.15, (σ1(1),σ2(1))=(1.7,2.1), and (σ1(2),σ2(2))=(1.4,0.9); we note η1=-0.1911<0, η2=-0.0210<0, and η3=-0.8543<0. The conditions of Theorem 8 are satisfied, so the prey and predator are both extinct (see Figure 3).
Next we only change the density of the white noise to(σ1(1),σ2(1))=(0.3,0.2) and (σ1(2),σ2(2))=(0.5,0.4) and keep μ=0.15. Simple calculation shows that η4=0.8622>0, η5=0.9622>0, and η6>0.6757>0. The conditions of Theorem 10 are satisfied, so the prey and predator are persistent (see Figure 4).
It is easy to see from Figures3 and 4 that the increase of the intensity of white noise can result in the extinction of prey and predator.Figure 3
Sample paths ofxt, yt with initial value x0=2.8 and y0=2.2.
(a)
k
=
1
(b)
k
=
2Figure 4
Sample paths ofxt, yt with initial value x0=0.4 and y0=0.3.
(a)
k
=
1
(b)
k
=
2Case 2.
Letμ=0.55, (σ1(1),σ2(1))=(0.6,0.5), and (σ1(2),σ2(2))=(0.45,0.55), which gives η1=-2.4356<0, η2=-0.3701<0, and η3=-1.2035<0. x(t) and y(t) are extinct (see Figure 5).
Next we only change the amount of toxicant toμ=0.12 and keep (σ1(1),σ2(1))=(0.6,0.5) and (σ1(2),σ2(2))=(0.45,0.55), We note η4=1.0467>0, η5=1.1467>0, and η6>0.7255>0. Thus x(t) and y(t) are persistent in mean (see Figure 6).
Figures5 and 6 show that the increase of the amount of toxicant can also result in the extinction of the prey and predator.Figure 5
Sample paths ofxt, yt with initial value x0=6 and y0=4.
(a)
k
=
1
(b)
k
=
2Figure 6
Sample paths ofxt, yt with initial value x0=1.6 and y0=1.4.
(a)
k
=
1
(b)
k
=
2



---

*Source: 3742197-2017-12-05.xml*