# A Stochastic Holling-Type II Predator-Prey Model with Stage Structure and Refuge for Prey

Authors: Wanying Shi; Youlin Huang; Chunjin Wei; Shuwen Zhang
Journal: Advances in Mathematical Physics
 (2021)

---

## Abstract

In this paper, we study a stochastic Holling-type II predator-prey model with stage structure and refuge for prey. Firstly, the existence and uniqueness of the global positive solution of the system are proved. Secondly, the stochastically ultimate boundedness of the solution is discussed. Next, sufficient conditions for the existence and uniqueness of ergodic stationary distribution of the positive solution are established by constructing a suitable stochastic Lyapunov function. Then, sufficient conditions for the extinction of predator population in two cases and that of prey population in one case are obtained. Finally, some numerical simulations are presented to verify our results.


---

## Body

## 1. Introduction

Population ecology is one of the most important research fields in biomathematics. There are three basic relationships among different species living in the same natural environment: (i) competitive relationship [1, 2], for example, crops and weeds compete for fertilizer, sunlight, and other resources for their own growth; (ii) reciprocal relationship [3, 4], such as rhinoceros birds feed on the parasites from rhinoceroses. For rhinoceroses, rhinoceros birds can get rid of the parasites and warn them when they are in danger; (iii) predation relationship [5, 6], as the saying goes, “big fish eat small fish, small fish eat shrimps.” Among the relationships of population interaction, the predator-prey system has been extensively studied. Lotka and Volterra put forward the predator-prey model in 1925 and 1926, respectively. Since then, many scholars have devoted themselves to the study of various predator-prey models and obtained abundant research results [7–10].In all the above predator-prey systems, it is generally assumed that predators in the same population have identical predation ability, and the prey in the same population has identical viability and fertility. However, in the real world, the life of many creatures should be divided into two stages: immature and mature. There are recognizable morphological and behavioral differences that may exist between these stages. For example, the immature creatures have no fertility and predation ability, and their survival ability and defense ability are relatively weak. On the contrary, the mature creatures not only have reproductive ability and predator ability but also have strong survival ability and defensive ability. These will have a great impact on the dynamic behavior of the population. Thus, it is more practical to study the predator-prey model with stage structure. In view of this, many scholars have studied assorted predator-prey systems with stage structure and given corresponding dynamic analysis [11–17]. In particular, the following predator-prey model with stage structure for prey has been considered in [16].
(1)ẋ1=αx2−r1x1−βx1−ηx12−β1x1x3,ẋ2=βx1−r2x2,ẋ3=x3−r+kβ1x1−η1x3,where α, r1, β, η, m, b, r2, r, k, and η1 are positive constants and ẋi=dxi/dt,i=1,2,3.x1, x2 and x3 represent the size of the immature prey population, mature prey population, and predator population, respectively. And r1, r2 and r indicate the mortality rates of immature prey population, mature prey population, and predator population, respectively. The parameter α denotes the birth rate of the immature prey population; β is the transformation rate of the mature prey population. η and η1 denote the intraspecific competition coefficient of the immature prey population and predator population, respectively. k is the digestion constant. Predators only prey on immature prey.Functional response plays an important role in describing the predator-prey model, which refers to the response of predators’ predation rate to the density of prey, that is, the predation effect of predators on prey. Model (1) uses bilinear functional response to describe the interaction between predator and prey. However, many researchers [18–20] have pointed out that nonlinear functional response can more accurately describe the trend of population density. Specifically, the growth rate of predators described by Holling-type II functional response is increasing with the increase of the number of prey [21]. When the survival space is far from saturated and the resources are sufficient, the growth rate of predators is relatively fast, but when the survival space and resources are limited, the growth rate changes less significantly and finally approaches a fixed value. Holling-type II functional response is according to some practical phenomena and takes into account the influence of density constraint. Therefore, it has been discussed by many scholars and extensively applied [13, 22–26]. In addition, it is well known that in many cases, prey can avoid predators through refuges. For example, crabs on the beach hide under sand or stones to avoid seabirds, thus increasing their survival chance. In recent years, different predator-prey systems with refuges have attracted much attention [6, 27–30]. In [27], Qi and Meng studied the threshold behavior of a stochastic predator-prey system with prey refuge and fear effect. They concluded that the survival rate of prey can be improved by increasing the strength of refuge. In [29], Ghosh et al. have studied a prey-predator model incorporating prey refuge and additional food for predators and found that slightly higher refuge is beneficial to the coexistence of species, but strong refuge will lead to the extinction of the predator population.Inspired by the above motivations, we consider the following Holling-type II predator-prey model with stage structure and refuge for prey:(2)ẋ1=αx2−r1x1−βx1−ηx12−m1−ba+x11−bx1x3,ẋ2=βx1−r2x2,ẋ3=x3−r+km1−ba+x11−bx1−η1x3,where m denotes the influence of predators on prey, a represents the half-saturation constant, b∈0,1 denotes the refuge rate of immature prey, 1−bx1 is the number of immature prey that predators can capture, and m1−b is the capture rate of predators.On the other hand, in nature, the biological population is inevitably affected by environmental noise. May [31] pointed out that due to the continuous fluctuation of the environment, the birth rate, mortality, environmental capacity, and other parameters in the model will show random fluctuations in varying degrees. However, the deterministic model does not consider the impact of environmental disturbance, so to some extent, it cannot accurately predict the dynamic behavior of the population. In order to describe the environmental noise better, many scholars consider the environmental noise as white noise, for instance, [22–27, 32, 33]. Considering the sensitivity of mortality to environmental noise, in this paper, we assume that the mortality rates r1, r2 and r are disturbed by environmental noise, i.e.,
(3)−r1⟶−r1+σ1Ḃ1t,−r2⟶−r2+σ2Ḃ2t,−r⟶−r+σ3Ḃ3t,where σii=1,2,3 are the intensities of the white noise and Biti=1,2,3 represent independent standard Brownian motions which are defined on the complete probability space Ω,F,Ft≥0,ℙ with a filtration Ft≥0 satisfying the usual conditions (that means it is increasing and right continuous while F0 contains all ℙ-null sets). Then, the stochastic model corresponding to the deterministic model (2) is derived in the following form:
(4)dx1=αx2−r1x1−βx1−ηx12−m1−ba+x11−bx1x3dt+σ1x1dB1t,dx2=βx1−r2x2dt+σ2x2dB2t,dx3=x3−r+km1−ba+x11−bx1−η1x3dt+σ3x3dB3t.Based on the analysis above, the paper focuses on the analysis of some dynamic behaviors of the stochastic Holling-type II predator-prey model with stage structure and refuge for prey. The major contributions are presented as follows. First, for the first time, the Holling-type II functional response is applied to establish a predator-prey model with stage structure and refuge for prey. Second, the existence and uniqueness as well as the stochastically ultimate boundedness of positive solution of the stochastic model are analyzed. Third, some conditions for the existence and uniqueness of ergodic stationary distribution of the stochastic model are established. Finally, some extinction conditions of the stochastic model are studied.Moreover, in [26], Xu et al. studied a stochastic two-predator one-prey model with modified Leslie-Gower and Holling-type II schemes. The authors proved the existence and uniqueness of global positive solution of their stochastic model by using comparison theorem. In this paper, we consider it by constructing a Lyapunov function. Besides, in [4], Liu et al. discussed a mutualism system in random environments. They obtained the existence and uniqueness of a stable stationary distribution by means of Markov semigroup theory and Fokker-Planck equation. In this paper, we consider the existence and uniqueness of an ergodic stationary distribution of the stochastic model (4) by the stochastic Lyapunov function method.The rest of this paper is arranged as follows. Some useful definitions and lemmas are given in Section2. The existence and uniqueness of the global positive solution are discussed in Section 3. In Section 4, we obtain sufficient conditions for stochastically ultimate bounded of the prey and predator. In Section 5, the sufficient conditions for the existence of a unique ergodic stationary distribution of the positive solution are established. In Section 6, we establish sufficient conditions for the extinction of the predator and prey in two cases. The first case is the extinction of the predator, and another case is that both the prey and the predator are extinct. To verify our results, some numerical simulations are presented in Section 7. Finally, the conclusion is given in Section 8.


## 2. Preliminaries

In this section, some definitions and lemmas are given to prepare for further work.Lemma 1 (Ito^ formula) (see [34]).
LetXt be a d-dimensional Ito^ process on t≥0 with the stochastic differential
(5)dXt=ftdt+gtdBt,where f∈L1ℝ+;ℝd and g∈L2ℝ+;ℝd×m. Let V∈C2,1ℝd×ℝ+;ℝ. Then, VXt,t is again an Ito^ process with the stochastic differential given by
(6)dVXtz,t=VtXt,t+VXXt,tft+12tracegTtVXXXt,tgtdt+VXXt,tgtdBt=LVXt,tdt+VXXt,tgtdBt,where Vt=∂V/∂t,VX=∂V/∂X1,∂V/∂X2,⋯,∂V/∂Xd,VXX=∂2V/∂Xi∂Xjd×d.Definition 2 (see [35]).
With respect to system (4), the solution is said to be stochastically ultimate bounded, if for ε∈0,1, there is a positive constant H=Hε such that for any initial data x10,x20,x30∈ℝ+3, the solution x1t,x2t,x3t has the property that
(7)limt⟶∞supℙxt≥H≤ε,where xt=x12+x22+x321/2.Definition 3 (see [36]).
The transition probability functionPs,x,t,A is said to be time-homogeneous (and the corresponding Markov process is called time-homogeneous) if the function Ps,x,t+s,A is independent of s, where 0≤s≤t,x∈ℝd and A∈B, B denotes the σ-algebra of Borel sets in ℝd.
LetXt be a regular time-homogeneous Markov process in ℝd described by the stochastic differential equation
(8)dXt=fXtdt+gXtdBt.
The diffusion matrix of the processXt is defined as follows:
(9)Ax=aijx,aij=gixgjx.Lemma 4 (see [36]).
The Markov processXt has a unique ergodic stationary distribution π· if there exists a bounded open domain U⊂ℝd with regular boundary Γ, having the following properties:
H1: the diffusion matrix Ax is strictly positive definite for all x∈U.
H2: there exists a nonnegative C2-function V such that LV is negative for any ℝd\U.


## 3. Existence and Uniqueness of the Global Positive Solution

The existence and uniqueness of global positive solution are the premises of studying the properties of population dynamics. From [34], for any given initial data x10,x20,x30∈ℝ+3, if the coefficients of SDE model (4) satisfy the linear growth condition and the local Lipschitz condition, then there exists a unique global positive solution x1t,x2t,x3t∈ℝ+3. But the coefficients of model (4) only satisfy the local Lipschitz condition; then, the solution of the system will explode in a finite time. In this section, we will prove that SDE model (4) has a unique global positive solution.Theorem 5.
For any given initial datax10,x20,x30∈ℝ+3, there is a unique solution x1t,x2t,x3t to system (4) on t≥0 and the solution will remain in ℝ+3 with probability 1.Proof.
Since the coefficients of system (4) satisfy the local Lipschitz condition, then for any given initial value x10,x20,x30∈ℝ+3, there is a unique local solution x1t,x2t,x3t∈ℝ+3 on t∈0,τe (see [34]), where τe denotes the explosion time. To show this solution is global, we only need to prove τe=∞a.s. Let n0>0 be sufficiently large such that xi0∈1/n0,n0i=1,2,3. For any integer n≥n0, we define the following stopping time:
(10)τn=inft∈0,τe:minx1t,x2t,x3t≤1normaxx1t,x2t,x3t≥n.We let inf∅=∞ (∅ denotes the empty set). From the definition of stopping time, it is easy to know that τn is increasing as n⟶∞. Set τ∞=limn⟶∞τn; hence, τ∞≤τea.s. If we can show that τ∞=∞a.s., then τe=∞a.s. and x10,x20,x30∈ℝ+3a.s. for all t≥0. Thus, to complete the proof, we only need to prove τ∞=∞a.s. If this statement is false, for any ε∈0,1, there exists a constant T>0 such that
(11)ℙτ∞≤T>ε.
Therefore, there is an integern1≥n0, for all n≥n1, we have
(12)ℙτn≤T≥ε.
Now, we define aC2-function V:
(13)Vx1,x2,x3=x1−c−clnx1c+αr2x2−1−lnx2+1kx3−1−lnx3,where c is a positive constant that will be determined later. Since
(14)u−c−clnuc,u−1−lnu≥0,∀u>0,thus, V is a nonnegative function. Applying Ito^ formula to V, we have
(15)dVx1,x2,x3=LVx1,x2,x3dt+σ1x1−cdB1t+σ2αr2x2−1dB2t+σ31kx3−1dB3t,where LV:ℝ+3⟶ℝ is defined by
(16)LVx1,x2,x3=1−cx1αx2−r1x1−βx1−ηx12−m1−ba+x11−bx1x3+αr21−1x2×βx1−r2x2+1kx3−1−r+km1−ba+x11−bx1−η1x3+c2σ12+α2r2σ22+12kσ32=αx2−r1+βx1−ηx12−m1−ba+x11−bx1x3−cαx2x1+cr1+β+cηx1+cm1−ba+x11−bx3+αβr2x1−αx2−αβr2x1x2+α−rkx3+m1−ba+x11−bx1x3−η1kx32+rk−m1−ba+x11−bx1+η1kx3+c2σ12+α2r2σ22+12kσ32=−ηx12+αβr2+cηx1−r1+βx1−η1kx32+η1kx3−rkx3+cm1−ba+x11−bx3−cαx2x1−αβr2x1x2+cr1+β+α−m1−ba+x11−bx1+rk+c2σ12+α2r2σ22+12kσ32≤−ηx12+αβr2+cηx1−η1kx32+η1kx3+cm1−ba−rkx3+cr1+β+α+rk+c2σ12+α2r2σ22+12kσ32.
Choosec=ar/km1−b, then
(17)LV≤−ηx12−αβr2+cηx1+12ηαβr2+cη2+14ηαβr2+cη2−η1kx32−x3+14+η14k+cr1+β+α+rk+c2σ12+α2r2σ22+12kσ32≤14ηαβr2+arkm1−b2+η14k+arr1+βkm1−b+α+rk+ar2km1−bσ12+α2r2σ22+12kσ32≕K.where K is a positive constant. The rest of the proof is similar to [37] and hence is omitted here. The proof is completed.


## 4. Stochastically Ultimate Boundedness

Theorem5 shows that the solution of system (4) remains in the positive cone ℝ+3. However, this nonexplosion property in a population dynamical system is often not good enough. Hence, we need to consider another important asymptotic property: stochastically ultimate boundedness, which means that the solution will be ultimately bounded with large probability.Theorem 6.
If the conditionsσ12<2r1+2β−1, σ22<2r2−1, and σ32<2r−2km−1 hold, then, the solution of system (4) is stochastically ultimately bounded for any initial data x10,x20,x30∈ℝ+3.Proof.
Forx1t,x2t,x3t∈ℝ+3, define
(18)Vx1,x2,x3=x12+x22+x32.
By Itô formula, we have(19)dVx1,x2,x3=LVx1,x2,x3dt+2σ1x12dB1t+2σ2x22dB2t+2σ3x32dB3t,where
(20)LV=2x1αx2−r1x1−βx1−ηx12−m1−ba+x11−bx1x3+2x2βx1−r2x2+2x32−r+km1−ba+x11−bx1−η1x3+σ12x12+σ22x22+σ32x32=−2ηx13+2αx1x2+σ12−2r1−2βx12−2m1−ba+x11−bx12x3+2βx1x2+σ22−2r2x22−2η1x33+σ32−2rx32+2km1−ba+x11−bx1x32≤σ12−2r1−2β+1x12+σ22−2r2+1x22+σ32−2r+2km+1x32+2α+βx1x2−x12−x22−x32.
Let(21)fx1,x2,x3=σ12−2r1−2β+1x12+σ22−2r2+1x22+σ32−2r+2km+1x32+2α+βx1x2.
According to the conditions, we can find that functionfx1,x2,x3 has an upper bound. We assume that its upper bound is as follows:
(22)M1=supx1,x2,x3∈ℝ+3fx1,x2,x3.
LetM=M1+1 and noticing f0,0,0=0, thus, M>0. According to (19), we can obtain
(23)dVx1,x2,x3≤N−x12+x22+x32dt+2σ1x12dB1t+2σ2x22dB2t+2σ3x32dB3t.
ByIto^ formula, we have
(24)detVx1,x2,x3=etVx1,x2,x3dt+etdVx1,x2,x3≤etx12+x22+x32dt+etM−etx12+x22+x32dt+2σ1etx12dB1t+2σ2etx22dB2t+2σ3etx32dB3t=etMdt+2σ1etx12dB1t+2σ2etx22dB2t+2σ3etx32dB3t.
Integrating both sides of (24) from 0 to t and then taking expectations, we get
(25)etEVx1,x2,x3≤Vx10,x20,x30+Met−M.
Hence, we have(26)limt⟶∞EVx1,x2,x3≤M,namely,
(27)limt⟶∞Ex12+x22+x32≤M.
For anyε>0, let H=M/ε. By Chebyshev’s inequality, we can obtain
(28)ℙxt≥H≤Ext2H2.
Then,(29)limt⟶∞supℙxt≥H≤MH2=ε.
The proof of Theorem6 is completed.□


## 5. Existence of Ergodic Stationary Distribution

In this section, we will study the sufficient conditions for the existence and uniqueness of an ergodic stationary distribution of the positive solutions of system (4). The existence of an ergodic stationary distribution means that all species can coexist for a long time and are stochastic weakly persistent [17].Theorem 7.
Assume thatr2>6σ12∨σ22∨σ32,r>4σ12∨σ22∨σ32 and 0≤b<1, then for any initial x10,x20,x30∈ℝ+3, system (4) has a unique ergodic stationary distribution.Proof.
According to Lemma4, we first prove that the condition H1 is true. The diffusion matrix of system (4) is given by
(30)A=σ12x12000σ22x22000σ32x32.
Clearly, the matrixA is strictly positive definite for all x1,x2,x3∈ℝ+3. Hence, the condition H1 in Lemma 4 holds.
Next, we verify that the conditionH2 holds. Let V1=x1+α/r2x2−lnx1. By Ito^ formula, we have
(31)LV1≤αx2−ηx12+αβr2x1−αx2−αx2x1+r1+β+ηx1+m1−bax3+12σ12=−ηx12+αβr2+ηx1−αx2x1+r1+β+m1−bax3+12σ12=−ηx12−1ηαβr2+ηx1+12ηαβr2+η2−14η2αβr2+η2−αx2x1+r1+β+m1−bax3+12σ12≤14ηαβr2+η2−αx2x1+r1+β+max3+12σ12.
DefineV2=V1−4αβ/r2+σ22/22lnx2−lnx3+m+aη1/arx3, then we have
(32)LV2≤14ηαβr2+η2−αx2x1+r1+β+max3+12σ12+r+η1x3+12σ32−4αβ2x1x2r2+σ22/22+4αβr2+σ22/2−m+aη1ax3+kmm+aη1a2rx1x3≤−24α2β2r2+σ22/22+14ηαβr2+η2+r1+β+12σ12+r+12σ32+4αβr2+σ22/2+kmm+aη1a2rx1x3=14ηαβr2+η2+r1+β+12σ12+r+12σ32+kmm+aη1a2rx1x3=λ+kmm+aη1a2rx1x3,where λ=1/4ηαβ/r2+η2+r1+β+1/2σ12+r+1/2σ32>0. Define
(33)V3=−lnx2+Mx3,where M>0 is a sufficiently large number. Then, we obtain
(34)LV3≤−βx1x2+r2+12σ22+kmMax1x3.
Define(35)V4=1θ+1x1+3αr2x2+x3kθ+1,where θ=4. Let
(36)U=x1+3αr2x2+x3k,then
(37)LV4=Uθαx2−r1x1−βx1−ηx12−m1−ba+x11−bx1x3+3αβr2x1−3αx2−rkx3+Uθm1−ba+x11−bx1x3−η1kx32+θ2Uθ−1σ12x12+3αr22σ22x22+1k2σ32x32≤3αr2x1x1+3αr2x2+x3kθ−2α3αr2θx2θ+1−ηx1θ+2−rkθ+1x3θ+1+θ2Uθ+1σ12∨σ22∨σ32=−η2x1θ+2−α3αr2θx2θ+1−r2kθ+1x3θ+1+B1,where
(38)B1=−η2x1θ+2−α3αr2θx2θ+1−r2kθ+1x3θ+1+3αr2x1x1+3αr2x2+x3kθ+θ2Uθ+1σ12∨σ22∨σ32.
From the assumptions thatB1 has an upper bound. Let B=supx1,x2,x3∈ℝ+3B1. Then,
(39)LV4≤−η2x1θ+2−α3αr2θx2θ+1−r2kθ+1x3θ+1+B.
Define aC2-function V¯:ℝ3⟶ℝ,
(40)V¯=NV2x1,x2,x3+V3x2+V4x1,x2,x3,where N<0 is a sufficiently small number satisfying the following condition:
(41)Nλ+B+r2+σ222≤−K−1.
K is a positive constant. Since V¯ is continuous and tends to ∞ as x1,x2,x3 approaches the boundary of ℝ+3. Hence, it must have a lower bound. Now, we assume that it gets this lower bound at point x¯1,x¯2,x¯3∈ℝ+3. Let
(42)V=NV2x1,x2,x3+V3x2+V4x1,x2,x3−V¯x¯1,x¯2,x¯3.
By (32), (34), (39), and (41), we have
(43)LV≤Nλ+kmNaη1+ma2rx1x3−βx1x2−η2x1θ+2−α3αr2θx2θ+1−r2kθ+1x3θ+1+kmMax1x3+B+r2+σ222≤Nλ−βx1x2−η2x1θ+2−α3αr2θx2θ+1−r2kθ+1x3θ+1+kmMax1x3+B+r2+σ222.
Define a bounded open setUε as follows:
(44)Uε=x1,x2,x3∈ℝ+3∣ε<x1<1ε,ε<x3<1ε,ε2<x2<1ε2.
0<ε<1 is a sufficiently small number. In the set ℝ+3\Uε, we can choose ε sufficiently small such that the following conditions hold:
(45)ε≤minaθ+1kmM,arθ+12kθ+2mM,aθ+2kmMθ+1,aθ+2η2kmM,(46)max−βε,−η41εθ+2,−r4kθ+11εθ+1,−α23αr2θε2θ+1≤−D−K,where D is a positive constant which will be determined later. We can divide ℝ+3\Uε into six domains:
(47)U1=x1,x2,x3∈ℝ3∣x1≤ε,U2=x1,x2,x3∈ℝ3∣x3≤ε,U3=x1,x2,x3∈ℝ3∣x1>ε,x2≤ε2,U4=x1,x2,x3∈ℝ3∣x1≥1ε,U5=x1,x2,x3∈ℝ3∣x3>1ε,U6=x1,x2,x3∈ℝ3∣x2>1ε.
Clearly,Uεc=ℝ+3\Uε=U1∪U2∪U3∪U4∪U5∪U6. Next, we will prove that LVx1,x2,x3≤−K for any x1,x2,x3∈Uεc.
For anyx1,x2,x3∈U1, due to x1x3≤εx3≤εθ+x3θ+1/θ+1=εθ/θ+1+εx3θ+1/θ+1, by (43) and (45), we have
(48)LV≤Nλ+kmMεθaθ+1−r2kθ+1−kmMεaθ+1x3θ+1+B+r2+σ222≤−K−1+1=−K<0.
For anyx1,x2,x3∈U2, since x1x3≤εx1≤εθ+1+x1θ+2/θ+2=εθ+1/θ+2+εx1θ+2/θ+2, by (43) and (45), we can get
(49)LV≤Nλ+kmMθ+1εaθ+2−η2−kmMεaθ+2x1θ+2+B+r2+σ222≤−K−1+1=−K<0.
From (43), we have
(50)LV≤−βx1x2−η4x1θ+2−α23αr2θx2θ+1−r4kθ+1x3θ+1−η4x1θ+2−α23αr2θx2θ+1−r4kθ+1x3θ+1+kmMax1x3+B+r2≤−βx1x2−η4x1θ+2−α23αr2θx2θ+1−r4kθ+1x3θ+1+D,(51)D=supx1,x2,x3∈ℝ+3−η4x1θ+2−α23αr2θx2θ+1−r4kθ+1x3θ+1+kmMax1x3+B+r2+σ222.
For anyx1,x2,x3∈U3, by (46) and (50), we obtain
(52)LV≤−βx1x2+D≤−βεε2+D=−βε+D≤−K<0.
For anyx1,x2,x3∈U4, by (46) and (50), we have
(53)LV≤−η4x1θ+2+D≤−η41εθ+2+D≤−K<0.
For anyx1,x2,x3∈U5, by (46) and (50), we obtain
(54)LV≤−r4kθ+1x3θ+1+D≤−r4kθ+11εθ+1+D≤−K<0.
For anyx1,x2,x3∈U6, by (46) and (50), we get
(55)LV≤−α23αr2θx2θ+1+D≤−α23αr2θε2θ+1+D≤−K<0.
Therefore, there is a sufficiently smallε, s.t. LVx1,x2,x3≤−K<0 for any x1,x2,x3∈Uεc. So the condition H2 in Lemma 4 is satisfied. Thus, system (4) has a unique stationary distribution and it has the ergodic property. We have completed the proof.□


## 6. Stochastic Extinction

In this section, we will establish sufficient conditions for the extinction of the prey and predator populations.Theorem 8.
Letx1t,x2t,x3t be the solution of system (4) with an initial value x10,x20,x30∈ℝ+3. If αβ/r2>r1+β and km1−bαβ/r2−r1−β/aη<r+σ3/2, then the predator population will die out exponentially with probability one, i.e.,
(56)limt⟶∞x3t=0a.s.Proof.
Define(57)V=lnx3+Cx1+αr2x2,where C is a positive constant that will be determined later. By Ito^ formula, we have
(58)LV=Cαx2−r1x1−βx1−ηx12−m1−ba+x11−bx1x3+αβr2x1−αx2−r+km1−ba+x11−bx1−η1x3−12σ32≤C−ηx12+αβr2−r1−βx1−r−12σ32+km1−bax1=C−ηx12−2ηαβr2−r1−βx1+1η2αβr2−r1−β2+Cηαβr2−r1−β2−Cαβr2−r1−βx1−r−12σ32+km1−bax1≤Cηαβr2−r1−β2−Cαβr2−r1−βx1−r−12σ32+km1−bax1.
LetC=km1−b/aαβ/r2−r1−β, then
(59)LV≤−r−12σ32+km1−bαβ/r2−r1−βaη.
ByIto^ formula,
(60)dVx1,x2,x3=LVx1,x2,x3dt+σ3dB3t+Cσ1x1dB1t+Cα2σ2x2dB2t.
And hence,(61)dVx1,x2,x3≤−r−12σ32+km1−bαβ/r2−r1−βaηdt+σ3dB3t+Cσ1x1dB1t+Cαr2σ2x2dB2t.
Integrating (61) from 0 to t and then dividing by t on both sides, we get
(62)Vt−V0t≤−r−12σ32+km1−bαβ/r2−r1−βaη+σ3B3tt+Cσ11t∫0tx1sdB1s+Cαr2σ21t∫0tx2sdB2s.
Similar to [38], we can get limt⟶∞1/t∫0tx1sdB1s=0 and limt⟶∞1/t∫0tx2sdB2s=0. Then,
(63)limt⟶∞supVtt≤−r−12σ32+km1−bαβ/r2−r1−βaη<0.
Thus,(64)limt⟶∞suplnx3tt<0.
Then, we obtainlimt⟶∞x3t=0.Remark 9.
Theorem8 shows that the larger noise intensity or high mortality of predators as well as high refuge rate of immature prey can lead to the extinction of the predator population. However, under these conditions, we cannot judge whether the immature and mature prey populations are extinct or not.Theorem 10.
Assume thatx1t,x2t,x3t is the solution of system (4) with initial data x10,x20,x30∈ℝ+3, if 2r1+σ122α−2r2−σ22<−α−r1−r22, then, all the populations will die out exponentially with probability one, i.e.,
(65)limt⟶∞x1t=0,limt⟶∞x2t=0,limt⟶∞x3t=0a.s.Proof.
LetV=lnx1+x2. By Ito^ formula, we have
(66)LVx1,x2=1x1+x2αx2−r1x1−ηx12−r2x2−m1−ba+x11−bx1x3−12σ12x12x1+x22−12σ22x22x1+x22−σ1σ2x1x2x1+x22≤12x1+x222x1+x2αx2−r1x1−r2x2−σ12x12−σ22x22.
Rewriting(67)2x1+x2αx2−r1x1−r2x2−σ12x12−σ22x22as the following form
(68)x1,x2−2r1−σ12α−r1−r2α−r1−r22α−2r2−σ22x1,x2T.
Let the matrix(69)A=−2r1−σ12α−r1−r2α−r1−r22α−2r2−σ22,it is clear that
(70)A1=−2r1−σ12<0,A2=A=−2r1−σ122α−2r2−σ22−α−r1−r22>0.
Thus,A is negative definite. Define λmax as the maximum eigenvalue of matrix A, then we have
(71)LVx1,x2=x1,x2Ax1,x2T≤−λmax2x1+x22x12+x22.
Thus, we obtain(72)dVx1,x2≤−λmax2x1+x22x12+x22dt+σ1x1x1+x2dB1t+σ2x2x1+x2dB2t≤−14λmaxdt+σ1x1x1+x2dB1t+σ2x2x1+x2dB2t.
Integrating (72) from 0 to t and then dividing by t on both sides, we get
(73)lnx1+x2t≤lnx10+x20t−14λmax+1t∫0tσ1x1x1+x2dB1t+1t∫0tσ2x2x1+x2dB2t.
By the strong law of large number of the martingale, we have(74)limt⟶∞suplnx1+x2t≤−14λmax<0.
According to the nonnegativity ofx1 and x2, we can obtain that
(75)limt⟶∞suplnx1t<0,limt⟶∞suplnx2t<0,which implies that limt⟶∞x1t=0,limt⟶∞x2t=0. In other words, the prey population goes to extinction exponentially with probability one. Because of limt⟶∞x1t=0, thus, ∀ε>0,∃T0=T0ω and a set Ωε⊂Ω such that ℙΩε>1−ε and
(76)km1−ba+x11−bx1≤km1−bax1≤km1−baε.
By Itô formula, we have(77)dlnx3=−r+km1−ba+x11−bx1−η1x3−12σ32dt+σ3dB3t≤−r+km1−baε−12σ32dt+σ3dB3t.
Integrating (77) from 0 to t and then dividing by t on both sides, we get
(78)lnx3t−lnx30t≤−r−12σ32+km1−baε+1t∫0tσ3dB3t.
Then, we obtain(79)limt⟶∞lnx3t≤−r−12σ32+km1−baε.
Letε⟶0, then limt⟶∞suplnx3t−lnx30/t≤−r−1/2σ32<0, which implies that limt⟶∞x3=0.□Remark 11.
According to Theorem10, we can find that the extinction of immature prey population leads to the extinction of predator population. This is true, since model (4) shows that the predator population has no additional food source.


## 7. Numerical Simulations

For the sake of confirming our theoretical results established in the previous sections, we numerically simulate stochastic system (4). We use Milstein’s higher order method mentioned in [39] to give numerical simulations. The discretization transformation of system (4) is given as follows:
(80)x1j+1=x1j+αx2j−r1x1j−βx1j−ηx1j2−m1−ba+x1j1−bx1jx3jΔt+σ1x1jΔtε1,j+σ122x1jε1,j2−1Δt,x2j+1=x2j+βx1j−r2x2jΔt+σ2x2jΔtε2,j+σ222x2jε2,j2−1Δt,x3j+1=x3j+x3j−r+km1−ba+x1j1−bx1j−η1x3jΔt+σ3x3jΔtε3,j+σ322x3jε3,j2−1Δt,where the time increment Δt>0.σi2>0i=1,2,3 are the intensities of the white noise, and εi,ji=1,2,3 denote mutually independent Gaussian random variables which follow distribution N0,1.First, we select the parameter values as follows: (i)α=0.9, r1=0.8, β=0.7, η=0.6, r=0.2, m=1, b=0.3, a=1, r2=0.25, k=0.95, η1=0.5. Here, we consider the strengths of white noise as σ1=0.1, σ2=0.1, and σ3=0.2. By simple calculation, we have that r2=0.25>6σ12∨σ22∨σ32=0.24, r=0.2>4σ12∨σ22∨σ32=0.16,0<b=0.3<1. Thus, under these circumstances, Theorem 7 holds. Then, we can obtain that system (4) has a unique ergodic stationary distribution. Figure 1 illustrates this. And from Figure 1, we can directly see that immature prey population, mature prey population, and predator population coexist for a long time and are stochastic persistent. In addition, Figure 1 shows that the trajectories of the solution for the perturbed system (4) will show random fluctuations (red lines) when the undisturbed system are smooth curves (blue lines). These results indicate that the population density does fluctuate more or less due to the influence of white noise.Figure 1
The upper column indicates the paths ofx1, x2, and x3 of system (4) with an initial value x10,x20,x30=0.63,0.8,0.12. The strengths of white noise as σ1=0.1,σ2=0.1,andσ3=0.2 and other parameter values are taken as (i). The red lines represent the solution of system (4), and the blue lines represent the solution of the corresponding undisturbed system (2). The under column displays the histogram of the probability density function of x1, x2, and x3 populations.Next, we choose the parameter values as (ii)α=0.9, r1=0.8, β=0.7, η=0.9, r=0.45, m=1, b=0.3, a=1, k=0.95, η1=0.5, r2=0.3, σ1=0.1, σ2=0.1, σ3=0.7. We can calculate that αβ/r2=2.1>r1+β=1.5, km1−bαβ/r2−r1+β/aη≈0.443<r+σ32/2=0.695. Therefore, the conditions of Theorem 8 are satisfied. By Theorem 8, we get that the predator population goes extinct exponentially with probability one. Figures 2(a)–2(c) show the simulation results. In this case, we can see that the predator population is extinct, but the prey population is persistent. In addition, Theorem 8 shows that large noise intensity will lead to the extinction of the predator population. Therefore, in order to illustrate the effect of noise intensity on predator population, we choose σ3=0.15 and σ3=0.5, respectively, and other parameters are taken as (ii). The simulation result is shown in Figure 2(d). From Figure 2(d), we can see that when the noise intensity is large, that is, when σ3=0.5 (see the cyan line of Figure 2(d)), the predator population is extinct, which can be directly proved by the conclusion of Theorem 8.Figure 2
Trajectories ofx1, x2, and x3 of system (4) with an initial value x10,x20,x30=0.63,0.8,0.12. The red lines and cyan line represent the solution of system (4), and the blue lines represent the solution of the corresponding undisturbed system (2).
(a)(b)(c)(d)In order to illustrate the conclusion of Theorem10, we choose the noise intensities as σ1=0.1, σ2=1.15, and σ3=0.7, respectively. And other parameters are taken as (i). By simple calculation, we have that 2r1+σ122α−2r2−σ22≈−0.0362<−α−r1−r22=−0.0225. Thus, the condition of Theorem 10 holds. By Theorem 10, we obtain that the prey and predator populations die out exponentially with probability one. The numerical simulation results are shown in Figure 3.Figure 3
The paths ofx1, x2, and x3 of system (4) with an initial value x10,x20,x30=0.63,0.8,0.12. The red lines represent the solution of system (4), and the blue lines represent the solution of the corresponding undisturbed system (2).Finally, we use the numerical simulation method to analyze the effect of prey refuge on system (4). The refuge rate b of immature prey varies from 0 to 1. b=0 means that the immature prey population has no refuge, that is to say, the immature prey population is completely exposed to the predatory environment of the predator population. On the contrary, b=1 indicates that the immature prey is totally protected by refuge. Now, we take b=0 and b=1, respectively, and other parameters are as (i). The numerical simulation results are shown in Figure 4. Figure 4 shows that when b=0, the immature prey population, mature prey population, and predator population are persistent, that is, the predator population and the prey population coexist, which can be verified by Theorem 7. But when b=1, the predator population becomes extinct, because the predator population has no additional food source. Theorem 8 also proves this result.Figure 4
Trajectories of system (4) for b=0 and b=1.
(a)(b)(c)


## 8. Conclusion

In this paper, a stochastic Holling-type II predator-prey model with stage structure and refuge for prey is analyzed. Firstly, the existence and uniqueness of the global positive solution for system (4) have been proved. Due to its nonexplosive property, we obtain the conditions for the stochastically ultimate boundedness of positive solution. Next, the sufficient conditions for the existence of a unique ergodic stationary distribution of the positive solution are established. Then, we obtain sufficient conditions for the extinction of predator population in two cases and that of prey population in one case. The theoretical results are proved by numerical simulations.The following conclusions are verified by numerical simulations. Figure1 shows that the existence of ergodic stationary distribution can allow all populations are coexistent and stochastic persistent in a long time. In Figure 2, it is shown that the large noise intensity will lead to the extinction of predator population, as shown in Figures 2(c) and 2(d), but the prey population can survive, as shown in Figures 2(a) and 2(b). Figure 3 shows that the extinction of immature prey population will lead to the extinction of predator population. This is due to the fact that predators do not have additional food sources. Finally, it can be seen from Figure 4(c) that if the immature prey population is completely protected by refuge, it will lead predator population towards extinction.In addition, based on our results, from the perspective of ecological protection, we can provide additional food to the predators to ensure that the prey population and predator population maintain a high-density level. On the one hand, it can reduce the predation rate of the predator population to the prey [30]; on the other hand, it can ensure that the predator population will not be completely extinct when the refuge rate of the prey population is high.



---

*Source: 9479012-2021-07-10.xml*