Stationary distribution of a stochastic predator–prey system with stage structure for prey

Physica A xxx (xxxx) xxx

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Stationary distribution of a stochastic predator–prey system with stage structure for prey Xin Zhao, Zhijun Zeng

∗

School of Mathematics and Statistics, Northeast Normal University, 5268 Renmin Street Changchun, Jilin 130024, PR China

article

info

Article history: Received 15 June 2019 Received in revised form 27 September 2019 Available online xxxx MSC: 34D20 34F05 92B05

a b s t r a c t In this paper, a predator–prey model with stage structure is proposed and analyzed in which preys are divided into juvenile and mature preys from a deterministic framework to a stochastic one. We first prove that the system which we investigate has a unique global positive solution. Furthermore, we obtain the sufficient criteria for the existence of stationary distribution and ergodicity by constructing an appropriate Lyapunov function, which means the species are permanent. Finally, we give the sufficient conditions for extinction of preys. © 2019 Published by Elsevier B.V.

Keywords: Stationary distribution Stage structure Stochastic predator–prey system

1. Introduction Predator-prey systems have been studied enormously in the past decades and will continue to be a dominant theme in theoretical ecology as well as in applied mathematics due to its universal existence and importance [1]. As mathematical models were constructed based on the experiments and observations, so they have been considered as a better way to understand these complex phenomena and were used to make predictions. The general predator–prey system can be described as follows

⎧ dx ⎪ ⎨ 1 = ax1 (t) − ψ (x1 (t), x2 (t))x2 (t), dt

⎪ ⎩ dx2 = c ψ (x1 (t), x2 (t))x2 (t) − bx2 (t), dt

where x1 (t) and x2 (t) represent the population densities of prey and predator respectively at time t, a, b and c are positive constants which denote the intrinsic growth rate of prey, the death rate of predator and the trophic efficiency ranging from 0 to 1. In particular, ψ (x1 (t), x2 (t)) is called the functional response [2–4] which describe the relationship between an individual predator and its prey in term of capture efficiency. It is also one of the familiar nonlinear factors that affect dynamical properties of mathematical models. Note that there are several well-known functional responses in the predator–prey systems in theoretical ecology and mathematical ecology. The function ψ (x1 (t), x2 (t)) includes some special cases ∗ Corresponding author. E-mail address: [email protected] (Z. Zeng). https://doi.org/10.1016/j.physa.2019.123318 0378-4371/© 2019 Published by Elsevier B.V.

Please cite this article as: X. Zhao and Z. Zeng, Stationary distribution of a stochastic predator–prey system with stage structure for prey, Physica A (2019) 123318, https://doi.org/10.1016/j.physa.2019.123318.

2

X. Zhao and Z. Zeng / Physica A xxx (xxxx) xxx

(1) Lotka–Volterra type [5,6] or Holling type I [7]:

ψ (x1 ) = u1 x1 , where u1 is a positive constant. (2) Holling type II [7,8]: u1 x 1 , ψ (x1 ) = x 1 + u2 where u1 and u2 are positive constants and u2 is called the half-saturation constant. (3) Holling type III [7]:

ψ (x1 ) =

u1 x21 x21 + u22

,

where u1 and u2 are positive constants. (4) Beddington–DeAngelis type [9,10]: u3 x 1 ψ (x1 , x2 ) = , u4 x 1 + u5 + x 2 u

u

where u3 is the predator capture rate and u4 is the handling time per prey item. 5 5 To a large extent, however, population also depend on the stage structure which may influence the outcomes of population evolutions. In the real world, almost all of the populations have the life history. Hence, they can be divided into juvenile stage and mature stage where the juvenile populations are raised by mature populations and the juvenile populations are not equipped to independent survival. In consequence, it is unrealistic to assume that all the populations have the same ability. Liu [11] considered a stage structure predator–prey model with Beddington–DeAngelis functional response and obtained the sufficient conditions for global asymptotic stability. Lu [12] investigated a stage-structured predator–prey model with predation over juvenile prey and established the threshold dynamics determined by the net reproductive number of the predator population. We ought to make allowance for the fact that the predator have the chance of just consuming mature prey due to the juvenile prey does not possess the capacity to forage which cause the juvenile prey have little opportunity to encounter the predator. Based on the above assumptions, in our proposed model, we consider stage structure on prey population to make it more realistic. Here, we prescribed a three-dimensional predator–prey model which consists of juvenile prey, mature prey and predator population as follows

⎧ dx1 2 ⎪ ⎪ ⎪ dt = α x2 − r1 x1 − β x1 − bx1 , ⎪ ⎪ ⎨ dx2 = β x1 − r2 x2 − λx2 x3 , ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎩ dx3 = x3 (−r3 + kλx2 − cx3 ),

(1)

dt where x1 = x1 (t) and x2 = x2 (t) denote the densities of juvenile prey population and mature prey population at time t, respectively, x3 = x3 (t) denotes the density of predator population at time t. All the parameters involved in system (1) are assumed to be positive from the viewpoint of ecology, α is interpreted as the rate of mature prey giving birth to newborn juvenile prey, β represents the rate of juvenile prey becoming mature prey, r1 , r2 and r3 are the natural death rates of juvenile prey population, mature prey population and predator population, respectively, b and c are the intra-specific competition rates of the juvenile prey population and predator population, k is the conversion rate of nutrition, λ can be defined as the inter-specific competition coefficient. In the natural world, since populations are inevitably infected by various environmental noises such as white noise, Markovian switching and Lévy noise (see [13–18] and the references therein). Particularly, white noise is a significant component in an ecosystem and extensive attention has been received by some authors. May [19] pointed out the fact that the birth rate, death rate, carrying capacity, competition coefficient and other parameters involved in the system should exhibit random fluctuation to a greater or lesser extent due to the environmental noise. Ji [20] considered a one predator and two independent preys system with stochastic perturbations and investigated the long time behavior of this system. Liu [21] proposed a stochastic predator–prey model with stage structure for predator and Holling type II functional response in which established sufficient conditions for the existence and uniqueness of an ergodic stationary distribution as well as obtained sufficient conditions for extinction of the predator populations. So far as our knowledge is concerned, a very little amount of work have been done with the stochastic predator–prey model with stage structure for prey. Motivated by the referred works, our objective of this paper is to investigate and analyze a three-dimensional predator– prey model which consists of predator population and prey population which is divided into two sub populations, one is juvenile prey and other is mature prey. We consider the environmental fluctuations in population dynamics and assume that the death rates of preys and predators are not constants but are subject to environmental noises. We suppose that r1 , r2 and r3 are stochastically perturbed with .

−r1 → r1 + σ1 B1 ,

.

− r2 → −r2 + σ2 B2 ,

.

− r3 → −r3 + σ3 B3 ,

Please cite this article as: X. Zhao and Z. Zeng, Stationary distribution of a stochastic predator–prey system with stage structure for prey, Physica A (2019) 123318, https://doi.org/10.1016/j.physa.2019.123318.

X. Zhao and Z. Zeng / Physica A xxx (xxxx) xxx

3

then the stochastic version corresponding to system (1) takes the following form:

⎧ 2 ⎪ ⎨ dx1 = (α x2 − r1 x1 − β x1 − bx1 )dt + σ1 x1 dB1 (t), dx2 = (β x1 − r2 x2 − λx2 x3 )dt + σ2 x2 dB2 (t), ⎪ ⎩ dx = x (−r + kλx − cx )dt + σ x dB (t), 3 3 3 2 3 3 3 3

(2)

where Bi (t) means the independent standard Brownian motion and σi2 denotes the intensity of the white noise (i = 1, 2, 3). Throughout this paper, unless otherwise specified, let (Ω , F , {Ft }t ≥0 , P) be a complete probability space with a filtration {Ft }t ≥0 satisfying the usual conditions (i.e., it is right continuous and increasing while F0 contains all P-null sets). Let B(t) defined on this complete probability space. If G is a vector or matrix, its transpose is denoted by GT . We define

Rd+ = x = (x1 , . . . , xd ) ∈ Rd : xi > 0, 1 ≤ i ≤ d .

{

}

In general, we consider the d-dimensional stochastic differential equation dx(t) = f (x(t), t)dt + g(x(t), t)dB(t) for t ≥ t0

(3)

d

with initial value x(0) = x0 ∈ R , where B(t) denotes a d-dimensional standard Brownian motion defined on the complete probability space (Ω , F , {Ft }t ≥0 , P). Denote by C 2,1 (Rd × [t0 , ∞); R+ ) the family of all real-valued nonnegative functions V (x, t) defined on Rd × [t0 , ∞) such that they are continuously twice differentiable in x and once in t. The differential operator L of Eq. (3) is defined by [22] L=

d d ∑ ∂ ∂ 1∑ T ∂2 + fi (x, t) + [g (x, t)g(x, t)]ij . ∂t ∂ xi 2 ∂ xi ∂ xj i,j=1

i=1

If L acts on a function V ∈ C 2,1 (Rd × [t0 , ∞]; R+ ), then LV (x, t) = Vt (x, t) + Vx (x, t)f (x, t) +

where Vt =

∂V ∂t

, Vx =

(

∂V ∂V , , . . . , ∂∂xV ∂ x1 ∂ x2 d

)

1 2

trace[g T (x, t)Vxx (x, t)g(x, t)],

, Vxx =

(

∂2V ∂ xi ∂ xj

) d×d

. According to Itoˆ ’s formula [22], if x(t) ∈ Rd , then

dV (x(t), t) = LV (x(t), t)dt + Vx (x(t), t)g(x(t), t)dB(t). This paper is split section wise as follows. In Section 2, we demonstrate system (2) has a unique global positive solution. In Section 3, we establish sufficient conditions for the existence and uniqueness of an ergodic stationary distribution of the positive solutions to system (2). In Section 4, we obtain sufficient conditions for extinction of the prey populations. 2. Global positive solution Theorem 1. For any given initial value (x1 (0), x2 (0), x3 (0)) ∈ R3+ , there is a unique solution X = (x1 (t), x2 (t), x3 (t)) of system (2) on t ≥ 0 and the solution will remain in R3+ with probability one. Proof. Since the coefficients of system (2) satisfy the local Lipschitz condition, then for any initial value (x1 (0), x2 (0), x3 (0)) ∈ R3+ , system (2) has a unique local solution (x1 (t), x2 (t), x3 (t)) ∈ R3+ on [0, τe ) a.s., where τe is the explosion time. We aim to prove this solution is global, i.e. τe = ∞ a.s. To this end, let n0 ≥ 0 be sufficiently large such that x1 (0), x2 (0) and x3 (0) are lying within the interval [ n1 , n0 ]. For each integer n ≥ n0 , we define the stopping time 0

{

τn = inf t ∈ [0, τe ) : min{x1 (t), x2 (t), x3 (t)} ≤

1 n

}

or max{x1 (t), x2 (t), x3 (t)} ≥ n .

Obviously, τn is increasing as n → ∞. Set τ∞ = limn→∞ τn , whence τ∞ ≤ τe a.s. Hence, we only need to prove that τ∞ = ∞ a.s. If the statement is not true, then there exist a pair of constants T > 0 and ϵ ∈ (0, 1) such that

P {τ∞ ≤ T } > ϵ. Thus, there exists an integer n1 ≥ n0 such that

P {τn ≤ T } ≥ ϵ,

∀ n ≥ n1 .

We define a C -function V¯ : R3+ → R+ by 2

V¯ (x1 , x2 , x3 ) = (x1 − 1 − ln x1 ) +

α r2

(x2 − 1 − ln x2 ) +

α r2 k

(x3 − 1 − ln x3 )

Please cite this article as: X. Zhao and Z. Zeng, Stationary distribution of a stochastic predator–prey system with stage structure for prey, Physica A (2019) 123318, https://doi.org/10.1016/j.physa.2019.123318.

4

X. Zhao and Z. Zeng / Physica A xxx (xxxx) xxx

The nonnegativity of this function can be seen from u − 1 − ln u ≥ 0,

∀u ≥ 0.

By using Itoˆ ’s formula, we can derive that LV¯ =

( 1−

α

+

1

)

x1

(α x2 − r1 x1 − β x1 − bx21 ) +

(

1

1−

r2 k

x3

)

+α+

2

x3 (−r3 + kλx2 − cx3 ) +

= α x2 − r1 x1 − β x1 − bx21 − αλ

σ12

ασ22

α r3

α x2 x1

α

+

( 1−

r2

x2

ασ22

(β x1 − r2 x2 − λx2 x3 ) +

2r2

ασ32 2r2 k

+ r1 + β + bx1 + αλ

)

1

σ12 2

αc

αβ

+

r2

α r3

x1 − α x2 −

αλ

αλ r2

x2 x3 −

αβ x1 r 2 x2

ασ32

αc

x2 + − x2 + x3 + r2 k 3 r2 k r2 r2 k 2r2 k ( ) ( ) 2 ασ32 σ ασ 2 αλ αβ αc 2 αc α r3 ≤ −bx21 + b + x1 − x3 + + x3 + r 1 + β + α + + 1 + 2 + r2 r2 k r2 r2 k r2 k 2 2r2 2r2 k

(

b+

≤

r2

αβ

x3 +

)2

r2 k

r2

−

2r2

(

αλ r2

+

4b

r2 k

+

x3 +

αc r2 k

x2 x3 −

)2 + r1 + β + α +

4α c

:= C0 ,

r2

α r3 r2 k

+

σ12 2

+

ασ22 2r2

+

ασ32 2r2 k

where C0 is a positive constant and

{ sup

x1 ∈(0,+∞)

(

−bx21 + b +

αβ

( ) }

r2

x1

b+

=

αβ

)2

r2

,

4b

{ sup

x3 ∈(0,+∞)

−

αc r2 k

x23 +

(

αλ r2

+

αc r2 k

r2 k

) } x3

(

=

αλ r2

+

4α c

αc r2 k

)2 .

The rest of the proof is similar to Theorem 2 of Liu [23], so we omit it here. This completes the proof. □ 3. Existence of ergodic stationary distribution As is well-known that a majority of stochastic models do not have the positive equilibria traditionally when the random factors are took into account in the deterministic systems. Therefore, a stochastic weak stability named stationary distribution has attracted many researchers [24–33] and it has been applied extensively in many fields. In this section, we shall consider there exists a unique ergodic stationary distribution of system (2) under certain conditions. To deduce the existence and uniqueness of a stationary distribution, the following lemma [34] is given which will be used in main results. Lemma 1. The Markov process X (t) has a unique ergodic stationary distribution ν (·) if there exists a bounded open set Ω of Rd with regular boundary U such that the following properties hold: A1. The diffusion matrix A(x) is strictly positive definite for all x ∈ Ω . A2. There exists a nonnegative C 2 − function V such that LV is negative on Rd \ Ω . Theorem 2. r2 >

If

σ22

(4)

2

and r3 +

σ32 2

<

kλβ (r1 + β ) br2

(

kλ(r1 + β ) r1 +

−

σ12 2

+β

)

αb

(5)

then system (2) has a unique ergodic stationary distribution for any initial value (x1 (0), x2 (0), x3 (0)) ∈ R3+ . Proof. Direct computation shows that the diffusion matrix of system (2) is given by

⎛ 2 2 σ1 x1 A=⎝ 0 0

0 σ22 x22 0

⎞

0 0 ⎠. σ32 x23

Apparently, the matrix A is positive definite for any compact subset of R3+ . We have verified that condition A1 in Lemma 1 is satisfied. Next, we define a C 2 − function V (x1 , x2 , x3 ) such that LV ≤ −1 on R3+ \ Ω , where Ω is an open bounded set. Please cite this article as: X. Zhao and Z. Zeng, Stationary distribution of a stochastic predator–prey system with stage structure for prey, Physica A (2019) 123318, https://doi.org/10.1016/j.physa.2019.123318.

X. Zhao and Z. Zeng / Physica A xxx (xxxx) xxx

5

Then we can show system (2) has a unique ergodic stationary distribution according to Lemma 1. Firstly, we define b

V1 = − ln x1 +

r1 + β

c + mλ

V4 = V3 +

x1 , V2 = − ln x3 +

kλ(r1 + β )

αb

lα

(

1

x3 , V5 = − ln x1 , V6 =

V1 , V3 = V2 − m ln x2 ,

x1 +

x2 +

lα

)θ +2 x3

,

r3 θ +2 r2 r2 k then V is defined as V = MV4 + V5 + V6 , where m > 0 will be defined later, l > 0 is a sufficient large number, M is a positive number, θ > 0 satisfies the following condition σ22

r2 −

0<θ <

2

σ22

r2 +

2

−

r2 l

−

r2 l

.

By using Itoˆ ’s formula, it yields that L V1 = −

≤−

α x2

σ12

+ r1 + β + bx1 +

x1

α x2

+ r1 + β +

x1

σ12 2

+

2

b r1 + β

bα

+

r1 + β

(α x2 − (r1 + β )x1 − bx21 )

x2 .

σ32

[ ( ) ] σ2 b2 x21 kλ(r1 + β ) α x2 αb − kλx2 + cx3 + − + r1 + β + 1 + x2 − 2 αb x1 2 r1 + β r1 + β ( ) σ12 kλ(r1 + β ) r1 + 2 + β σ2 kλ(r1 + β )x2 ≤ r3 + 3 + − + cx3 . 2 αb bx1

LV2 = r3 +

LV3 = r3 +

σ

√ ≤ −2 We plug m =

(

2 3

2

kλ(r1 + β ) r1 +

+

2

+β

)

b

+ r3 +

kλ(r1 + β )x2

−

αb

kλ(r1 + β )mβ

kβλ(r1 +β )

σ12

σ32 2

+ cx3 −

bx1

(

kλ(r1 + β ) r1 +

+ mr2 +

σ12 2

m β x1 x2

+β

+ mr2 + mλx3

)

αb

+ (c + mλ)x3 .

into the above inequality and we also set

br22

kβλ2 (r1 + β ) + bcr22

µ0 := c + mλ =

br22

,

then we obtain L V3 ≤ −

kβλ(r1 + β )

+ r3 +

br2

σ32 2

(

kλ(r1 + β ) r1 +

+

σ12 2

+β

)

αb

+ µ 0 x3 .

In the same way, one can get

L V4 = −

kβλ(r1 + β ) br2

+ r3 +

σ32 2

(

kλ(r1 + β ) r1 +

+

⎛ σ2 ⎜ kβλ(r1 + β ) − r3 − 3 − ≤ −⎝ br2

σ12

+β

2

αb ( kλ(r1 + β ) r1 +

c + mλ + (c + mλ)x3 + x3 (−r3 + kλx2 − cx3 ) r3 )⎞ +β ⎟ kλµ0 x2 x3 . ⎠+

σ12 2

αb

2

)

r3

We define

ρ :=

kβλ(r1 + β ) br2

− r3 −

σ32 2

(

kλ(r1 + β ) r1 +

−

αb

σ12 2

+β

) ,

noting that ρ > 0 from the condition (5). Therefore, we arrive at LV4 ≤ −ρ +

kλµ0 r3

x2 x3 .

Please cite this article as: X. Zhao and Z. Zeng, Stationary distribution of a stochastic predator–prey system with stage structure for prey, Physica A (2019) 123318, https://doi.org/10.1016/j.physa.2019.123318.

6

X. Zhao and Z. Zeng / Physica A xxx (xxxx) xxx

Applying the generalized Itoˆ ’s formula, we have

α x2

L V5 = −

σ12

+ r1 + β + bx1 +

x1

.

2

Similarly, it is easy to obtain

( LV6 ≤

x1 +

+

r2

+

x1

x1 +

θ +1

lα

x1 +

(

r2

lα r2

x1 +

b

= − xθ1+3 −

r2

+

)θ +1

r2 k

x3

)θ

lα r2 k

θ +3

x3

+

lα

( x1 +

r2

lαβ r2

l2 α 2 r22

xθ2+2 −

σ22 x22 +

r2

x2 +

lα

( x1

x1 +

lα

)θ

r2 k

r2

σ22 x22 +

θ +1

xθ3+3 2r θ +2 kθ +2

l2 α 2

x3

lα

σ

r22

r2 k

2 2 2 x2

+

x1 +

lα r2

x2 +

( x1 +

2

)θ

lα r2 k

x3

l2 α 2 r22 k2

lα r2

x2 +

σ32 x23 −

lθ +1 α θ +2 (l − 1)(1 − θ )

2

r2θ +1

θ +1

+

2

( x1 +

2

θ +1

b

)θ +1 θ +1

)

x2θ+2 +

− xθ1+3 −

x3

)

r22 k

(

2

lθ +2 α θ +2 c

x2 +

σ 2 x2 2 3 3

r2θ +1

2

lαβ

r2 k

x23

l2 α 2

lθ +1 α θ +2 (l − 1)

− bxθ1+3 −

r22

lα c

x1 − α (l − 1)x2 −

σ12 x21 +

l2 α 2

x3

r2

2r2θ +2 kθ +2 2

r2 k

x3

θ +1

lθ +2 α θ +2 c

θ +1

)θ (

lα

lα

x2 +

−bx21 +

lθ +1 α θ +2 (l − 1)θ

2

−

x2 +

x2 +

lα

(

2

x3

r2 k

(

2

)θ +1 (

lα

x2 +

r2

θ +1

lαβ

≤

lα

lα r2

( x1 +

x2 +

lα r2

x2 +

)θ

lα r2 k

x3

)θ

lα r2 k

l2 α 2 r22 k2

x3

lα r2 k

)θ x3

σ12 x21

lθ+2 α θ +2 c r2θ +2 kθ +2

x3θ+3

xθ2+2

σ12 x21

σ32 x23 .

We notice that 0<θ <

r2 − r2 +

σ22 2

σ22 2

−

r2 l

−

r2 l

,

that is to say, 1 1−θ

<

lσ22 + 2r2 l − 2r2 2lσ22

.

It is equivalent to 1+θ 1−θ

<

2r2 (l − 1) lσ22

.

Through simple calculation, we get lθ+1 α θ +2 (l − 1)(1 − θ ) θ +1

r2

>

θ + 1 lθ α θ l2 α 2 σ22 r2θ

2

r22

,

(6)

then we have b θ +3 lθ +1 α θ +2 (l − 1)θ θ +2 lθ +2 α θ +2 c θ +3 L V6 ≤ − x 1 − x2 − θ +2 x3 + B, θ + 1 2 r2 2r2 kθ +2 which follows from condition (6), where

{ B=

θ +3

− x1

sup

2

(x1 ,x2 ,x3 )∈R3+

G(x1 , x2 , x3 ) =

b

lαβ r2

+

−

( x1

θ +1 2

x1 +

lθ +1 α θ +2 (l − 1)(1 − θ ) r2θ +1 lα r2

( x1 +

x2 +

lα r2

lα r2 k

x2 +

)θ +1 x3

lα r2 k

+

)θ x3

x2

2

lθ +2 α θ +2 c

−

θ +1

l2 α 2 r22

θ +2

2r2θ +2 kθ +2

( x1 +

σ22 x22 +

lα r2

θ +1 2

} θ +3

x3

x2 +

+ G(x1 , x2 , x3 ) ,

lα r2 k

( x1 +

)θ x3

lα r2

x2 +

σ12 x21 lα r2 k

)θ x3

l2 α 2 r22 k2

σ32 x23 .

Please cite this article as: X. Zhao and Z. Zeng, Stationary distribution of a stochastic predator–prey system with stage structure for prey, Physica A (2019) 123318, https://doi.org/10.1016/j.physa.2019.123318.

X. Zhao and Z. Zeng / Physica A xxx (xxxx) xxx

7

As a result, it leads to

( ) σ2 kλµ0 α x2 LV ≤ M −ρ + x2 x3 − + r1 + β + bx1 + 1 r3

x1

θ +1 θ +2

α

l

b

− xθ1+3 −

(l − 1)θ

r2θ +1

2

2

θ +2 θ+2

l

xθ2+2 −

α

c

xθ3+3 + B.

2r2θ +2 kθ +2

We aim to prove that LV ≤ −1 on R3+ \ Ω , once this is established, Theorem 2 follows from Lemma 1, where Ω is an open bounded set defined as

} { 1 1 1 . Ω = (x1 , x2 , x3 ) ∈ R3+ : ε 2 < x1 < 2 , ε < x2 < , ε < x3 < ε ε ε ⋃ c⋃ c⋃ c⋃ c⋃ c Next, we divide Ω c into six domains: Ω c = Ω1c Ω2 Ω3 Ω4 Ω5 Ω6 to show LV ≤ −1 actually holds, in which { } c 3 2 1. Ω1 = {(x1 , x2 , x3 ) ∈ R+ : 0 < x1 ≤ } ε , x2 > ε , 2. Ω2c = (x1 , x2 , x3 ) ∈ R3+ : x1 ≥ ε12 , { } 3. Ω3c = (x1 , x2 , x3 ) ∈ R3+ : 0 < x2 ≤ ε , { } 4. Ω4c = (x1 , x2 , x3 ) ∈ R3+ : x2 ≥ 1ε , { } 5. Ω5c = (x1 , x2 , x3 ) ∈ R3+ : 0 < x3 ≤ ε , { } 6. Ω6c = (x1 , x2 , x3 ) ∈ R3+ : x3 ≥ 1ε . While for convenience, we set K := r1 + β +

− M ρ + S1 −

θ +1 θ +2

α l − ε

α

(l − 1)θ

θ +1

2r2

σ12 2

+ B and 0 < ε < 1 is a constant satisfying the following conditions

ε θ+2 < −1,

(7)

b

− M ρ + S2 − ε −2θ−6 < −1,

(8)

4

− M ρ + S3 + lθ +2 α θ +2 c 2r θ +2 kθ+2

>

2

− M ρ + S1 − − M ρ + S4 +

Mkλµ0 ε θ + 2

Mkλµ0 ε r3 (θ + 3)

Mkλµ0 ε r3 (θ + 2)

lθ +2 α θ +2 c 4r2θ +2 kθ +2

sup (x1 ,x2 ,x3 )∈R3+

sup

b

θ +3

− x1 4

{ sup

θ +3

2

(x1 ,x2 ,x3 )∈R3+

(x1 ,x2 ,x3 )∈R3+

b

b

θ +3

− x1 2

ε −θ −2 < −1,

(11)

< −1,

(12)

,

(13)

ε −θ −3 < −1.

− x1

{

S3 =

θ +3

r3

{

S2 =

(10)

Mkλµ0 ε θ + 2

>

(9)

,

2r2θ +1

r2θ +1

S1 =

< −1,

lθ +1 α θ +2 (l − 1)θ

lθ +1 α θ +2 (l − 1)θ

− M ρ + S5 −

θ +3

r3

+ bx1 −

+ bx1 −

+ bx1 −

(14)

lθ +1 α θ +2 (l − 1)θ 2r2θ +1 lθ +1 α θ +2 (l − 1)θ r2θ +1 lθ +2 α θ +2 c 2r2θ +2 kθ +2

θ +3

x3

θ +2

x2

θ+2

x2

+

+

Mkλµ0

+

Mkλµ0

2r3

2r3

x22

x22

Mkλµ0 ε xθ3+3 r3

θ +3

−

−

lθ +2 α θ +2 c 2r2θ +2 kθ +2 lθ +2 α θ +2 c 2r2θ +2 kθ +2

θ +3

x3

θ +3

x3

+

Mkλµ0

+

Mkλµ0

2r3

2r3

} x23

+K ,

} x23

+K ,

} +K ,

Please cite this article as: X. Zhao and Z. Zeng, Stationary distribution of a stochastic predator–prey system with stage structure for prey, Physica A (2019) 123318, https://doi.org/10.1016/j.physa.2019.123318.

8

X. Zhao and Z. Zeng / Physica A xxx (xxxx) xxx

{ S4 =

θ +3

− x1

sup

{ sup (x1 ,x2 ,x3 )∈R3+

b

lθ +1 α θ +2 (l − 1)θ

+ bx1 −

2

(x1 ,x2 ,x3 )∈R3+

S5 =

b

θ +3

− x1

r2θ +1 lθ +1 α θ +2 (l − 1)θ

+ bx1 −

2

r2θ +1

θ +2

x2

θ +2

x2

+

Mkλµ0 ε xθ2+3

+

Mkλµ0

}

θ +3

r3

2r3

x22

−

+K ,

lθ +2 α θ +2 c 4r2θ +2 kθ +2

θ +3

x3

+

Mkλµ0 2r3

} x23

+K .

I. If (x1 , x2 , x3 ) ∈ Ω1c , LV ≤ −M ρ + K +

≤ −M ρ + K + ≤ −M ρ + S1 −

Mkλµ0 r3

α x2

x2 x3 −

x1

Mkλµ0 x22 + x23 r3

2 θ +1 θ +2

α l − ε

α

−

b

lθ+1 α θ +2 (l − 1)θ

2

r2θ +1

+ bx1 − xθ1+3 −

xθ2+2 −

lθ +2 α θ +2 c 2r2θ +2 kθ +2

xθ3+3

2lθ +1 α θ +2 (l − 1)θ θ +2 lθ +2 α θ +2 c θ +3 α b − − xθ1+3 + bx1 − x x3 2 ε 2 2r2θ +1 2r2θ +2 kθ +2

(l − 1)θ

2r2θ+1

ε θ +2

< −1, which follows from (7). Thus, we get LV < −1 for any (x1 , x2 , x3 ) ∈ Ω1c . II. If (x1 , x2 , x3 ) ∈ Ω2c , LV ≤ −M ρ + K +

≤ −M ρ + K +

Mkλµ0 r3

x2 x3 + bx1 −

Mkλµ0 x22 + x23 r3

2

b 2

xθ1+3 −

lθ +1 α θ +2 (l − 1)θ θ +1

r2

b

b

+ bx1 − xθ1+3 − xθ1+3 − 4

xθ2+2 −

lθ +2 α θ +2 c 2r2θ +2 kθ+2

lθ +1 α θ +2 (l − 1)θ

4

θ +1

r2

xθ3+3

xθ2+2 −

lθ +2 α θ +2 c 2r2θ +2 kθ +2

xθ3+3

b

≤ −M ρ + S2 − ε −2θ −6 < −1,

4

which follows from (8). Then we can obtain that LV < −1 for any (x1 , x2 , x3 ) ∈ Ω2c . III. If (x1 , x2 , x3 ) ∈ Ω3c , LV ≤ −M ρ + K +

Mkλµ0

ε x3 + bx1 − r3 ( Mkλµ0 ε θ + 2 ≤ −M ρ + K + + r3 θ +3 ≤ −M ρ + S3 +

b 2

xθ1+3 −

xθ3+3

θ +3

lθ +2 α θ +2 c 2r2θ +2 kθ +2

)

xθ3+3

b

lθ +2 α θ +2 c

2

2r2θ +2 kθ +2

+ bx1 − xθ1+3 −

xθ3+3

Mkλµ0 ε θ + 2

θ +3

r3

< −1,

which follows from (9) and (10). It is easy to obtain LV < −1 for any (x1 , x2 , x3 ) ∈ Ω3c . IV. If (x1 , x2 , x3 ) ∈ Ω4c , LV ≤ −M ρ + K +

≤ −M ρ + K + ≤ −M ρ + S1 −

Mkλµ0 r3

x2 x3 + bx1 −

Mkλµ0 x22 + x23 r3

2

2

xθ1+3 − b

lθ +1 α θ +2 (l − 1)θ

+ bx1 − xθ1+3 −

lθ +1 α θ +2 (l − 1)θ 2r2θ +1

b

2

θ +1

r2

xθ2+2 −

2lθ +1 α θ +2 (l − 1)θ θ +1

2r2

lθ +2 α θ +2 c 2r2θ +2 kθ+2

xθ2+2 −

xθ3+3

lθ +2 α θ +2 c 2r2θ +2 kθ+2

xθ3+3

ε −θ −2

< −1, which follows from (11). It can be derived that LV < −1 for any (x1 , x2 , x3 ) ∈ Ω4c . Please cite this article as: X. Zhao and Z. Zeng, Stationary distribution of a stochastic predator–prey system with stage structure for prey, Physica A (2019) 123318, https://doi.org/10.1016/j.physa.2019.123318.

X. Zhao and Z. Zeng / Physica A xxx (xxxx) xxx

9

V. If (x1 , x2 , x3 ) ∈ Ω5c , LV ≤ −M ρ + K +

Mkλµ0

ε x2 + bx1 − r3 ( Mkλµ0 ε θ + 1 ≤ −M ρ + K + + r3 θ +2 ≤ −M ρ + S4 +

b 2

xθ1+3 −

xθ2+2

lθ +1 α θ +2 (l − 1)θ

)

θ +2

r2θ +1

xθ2+2

b

lθ +1 α θ +2 (l − 1)θ

2

r2θ +1

− xθ1+3 + bx1 −

xθ2+2

Mkλµ0 ε θ + 1

θ +2

r3

< −1,

which follows from (12) and (13). Therefore, we have LV < −1 for any (x1 , x2 , x3 ) ∈ Ω5c .

VI. If (x1 , x2 , x3 ) ∈ Ω6c , LV ≤ −M ρ + K +

≤ −M ρ + K + ≤ −M ρ + S5 −

Mkλµ0 r3

x2 x3 + bx1 −

Mkλµ0 x22 + x23 r3

2

lθ +2 α θ+2 c 4r2θ +2 kθ +2

b 2

3 xθ+ − 1

lθ +1 α θ +2 (l − 1)θ r2θ+1

xθ2+2 −

b

lθ +1 α θ +2 (l − 1)θ

2

r2θ +1

− xθ1+3 + bx1 −

lθ +2 α θ +2 c 2r2θ +2 kθ +2

xθ2+2 −

xθ3+3

2lθ +2 α θ +2 c 4r2θ +2 kθ +2

xθ3+3

ε −θ −3

< −1, which follows from (14). Then it holds that LV < −1 for any (x1 , x2 , x3 ) ∈ Ω6c .

Clearly, we came to the conclusion that LV ≤ −1 on R3+ \ Ω for a sufficiently small ε . According to Lemma 1, we can get that system (2) is ergodic and has a unique stationary distribution. This completes the proof. □

4. Extinction Theorem 3.

Let (x1 (t), x2 (t), x3 (t)) be the solution of system (2) with any initial value (x1 (0), x2 (0), x3 (0)) ∈ R3+ . If

√

√

min{r1 + β, r2 }( R − 1)I{√R≤1} + max{r1 + β, r2 }( R − 1)I{√R>1} − (2(σ1−2 + σ2−2 ))−1 < 0, where R=

αβ r2 (r1 + β )

,

then the prey populations will die out, that is to say, lim x1 (t) = 0,

lim x2 (t) = 0.

t →∞

t →∞

Proof. The train of thought in this section is inspired by Liu [21]. Let

√

R(ω1 , ω2 ) = (ω1 , ω2 )Q ,

where (ω1 , ω2 ) =

( Q =

0 β r2

(√

R,

α r1 +β

0

)

α r1 +β

)

and

.

We define a C 2 − function V˜ : R2+ → R+ as V˜ (x1 , x2 ) = px1 + qx2 , Please cite this article as: X. Zhao and Z. Zeng, Stationary distribution of a stochastic predator–prey system with stage structure for prey, Physica A (2019) 123318, https://doi.org/10.1016/j.physa.2019.123318.

10

X. Zhao and Z. Zeng / Physica A xxx (xxxx) xxx

where p =

ω1 r1 +β

, q= p

d(ln V˜ ) =

V˜

+

ω2 r2

. By using Itoˆ ’s formula, we can derive that

[α x2 − (r1 + β )x1 − bx21 ]dt + q V˜

(β x1 − r2 x2 − λx2 x3 )dt + p

= LV˜ dt +

V˜

σ1 x1 dB1 (t) +

q

q V˜

p V˜

σ1 x1 dB1 (t)

σ2 x2 dB2 (t) −

1 2V˜ 2

(p2 σ12 x21 + q2 σ22 x22 )dt

σ2 x2 dB2 (t),

V˜

where LV˜ =

p

[α x2 − (r1 + β )x1 − bx21 ] +

V˜

q V˜

(β x1 − r2 x2 − λx2 x3 ) −

1 2V˜ 2

(p2 σ12 x21 + q2 σ22 x22 ).

In addition, one can get V˜ 2 =

(

p σ1 x1

1

σ1

)2

1

+ qσ2 x2

≤ (p2 σ12 x21 + q2 σ22 x22 )

σ2

(

1

σ12

+

1

)

σ22

,

and 1 V˜

[pα x2 − p(r1 + β )x1 − pbx21 + qβ x1 − qr2 x2 − qλx2 x3 ]

= ≤ = = = =

1 V˜ 1

(−pbx21 − qλx2 x3 ) +

1 V˜

[pα x2 − p(r1 + β )x1 + qβ x1 − qr2 x2 ]

[pα x2 − p(r1 + β )x1 + qβ x1 − qr2 x2 ] V˜ ( ) ω2 β 1 ω1 α x2 − ω1 x1 + x1 − ω2 x2 r2 V˜ r1 + β ( ) 1 (ω , ω ) Q (x1 , x2 )T − (x1 , x2 )T ˜V 1 2 ) 1 (√ R − 1 (ω1 x1 + ω2 x2 ) V˜ ) (√ 1 R − 1 (p(r1 + β )x1 + qr2 x2 ) px1 + qx2

≤ min{r1 + β, r2 }

(√

)

R − 1 I{√R≤1} + max{r1 + β, r2 }

(√

)

R − 1 I{√R>1} .

This implies that LV˜ ≤ min{r1 + β, r2 }

(√

)

R − 1 I{√R≤1} + max{r1 + β, r2 }

It follows that

[

d(ln V˜ ) ≤ min{r1 + β, r2 }

+

p

(√

σ1 x1 dB1 (t) +

q

)

(√

)

R − 1 I{√R>1} − (2(σ2−2 + σ3−2 ))−1 .

R − 1 I{√R≤1} + max{r1 + β, r2 }

(√

)

]

R − 1 I{√R>1} − (2(σ1−2 + σ2−2 ))−1 dt

σ2 x2 dB2 (t).

(15)

V˜ V˜ By integrating it from 0 to t and dividing by t on both sides of (15), we get ln V˜ (t) − ln V˜ (0) t

(√ (√ ) ) R − 1 I{√R≤1} + max{r1 + β, r2 } R − 1 I{√R>1} − (2(σ1−2 + σ2−2 ))−1 ≤ min{r1 + β, r2 } ∫ ∫ 1 t pσ1 x1 (s) 1 t qσ2 x2 (s) + dB1 (s) + dB2 (s) t 0 t V˜ (s) V˜ (s) (√ ) 0 (√ ) ≤ min{r1 + β, r2 } R − 1 I{√R≤1} + max{r1 + β, r2 } R − 1 I{√R>1} − (2(σ1−2 + σ2−2 ))−1 +

˜ M(t) t

+

˜ N(t) t

, (16)

˜ where M(t) :=

∫t 0

pσ1 x1 (s) V˜ (s)

˜ dB1 (s) and N(t) :=

˜ , M(t) ˜ ⟩t = σ12 ⟨M(t)

∫ t( 0

px1 (s) V˜ (s)

)2

∫t 0

qσ2 x2 (s) V˜ (s)

ds ≤ σ12 t ,

dB2 (s) are local martingales with quadratic variations as follows

˜ , N(t) ˜ ⟩t = σ22 ⟨N(t)

∫ t( 0

x2 (s) V˜ (s)

)2

ds ≤ σ22 t .

Please cite this article as: X. Zhao and Z. Zeng, Stationary distribution of a stochastic predator–prey system with stage structure for prey, Physica A (2019) 123318, https://doi.org/10.1016/j.physa.2019.123318.

X. Zhao and Z. Zeng / Physica A xxx (xxxx) xxx

11

Applying the strong law of large numbers for local martingales, we obtain lim

˜ M(t)

t →∞

t

= 0,

lim

t →∞

˜ N(t) t

= 0 a.s.

Taking the superior limit on both sides of (16), one can get ln V˜ (t)

lim sup

t

t →+∞

≤ min{r1 + β, r2 }

(√

)

R − 1 I{√R≤1} + max{r1 + β, r2 }

(√

)

R − 1 I{√R>1} − (2(σ1−2 + σ2−2 ))−1

< 0, which is equivalent to lim sup

ln x1 (t)

t →+∞

t

< 0,

lim sup

ln x2 (t)

t →+∞

t

< 0.

It implies that lim x1 (t) = 0,

t →∞

lim x2 (t) = 0,

t →∞

which means the extinction of preys. This completes the proof. □ Acknowledgment This work was supported by the Fundamental Research Funds for the Central Universities, People’s Republic of China (Grant No. 2412017FZ004). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26]

A.A. Berryman, The origin and evolution of predator–prey theory, Ecology 73 (1992) 1530–1535. R.M. May, Stability and Complexity in Model Ecosystems, Princeton University Press, Princeton, 1974. J.D. Murray, Mathematical Biology, Springer-Verlag, Berlin, 1989. R.M. May, Stability and Complexity in Model Ecosystems, Princeton University Press, Princeton, 1974. A.J. Lotka, Elements of Physical Biology, William and Wilkins, Baltimore, 1925, Reissued as Elements of Mathematical Biology, Dover, New York, 1956. V. Volterra, Fluctuations in the abundance of a species considered mathematically, Nature 118 (1926) 558–560. C.S. Holling, The functional response of predator to prey density and its role in mimicry and population regulation, Mem. Entomol. Soc. Can. 45 (1965) 1–60. Q. Liu, L. Zu, D. Jiang, Dynamics of stochastic predator–prey models with Holling II functional response, Commun. Nonlinear Sci. Numer. Simul. 37 (2016) 62–76. D.L. DeAngelis, R.A. Goldsten, R. Neill, A model for trophic interaction, Ecology 56 (1975) 881–892. J.R. Beddington, Mutual interference between parasites or predators and its effect on searching efficiency, J. Animal. Ecol. 44 (1975) 331–341. M. Liu, K. Wang, Global stability of stage-structured predator–prey models with Beddington–DeAngelis functional response, Commun. Nonlinear Sci. Numer. Simul. 16 (9) (2011) 3792–3797. Y. Lu, K.A. Pawelek, S.Q. Liu, A stage-structured predator–prey model with predation over juvenile prey, Appl. Math. Comput. 297 (2017) 115–130. P. Wang, J. Feng, H. Su, Stabilization of stochastic delayed networks with Markovian switching and hybrid nonlinear coupling via aperiodically intermittent control, Nonlinear Anal.-Hybrid Syst. 32 (2019) 115–130. Y. Xu, H. Zhou, W. Li, Stabilisation of stochastic delayed systems with lévy noise on networks via periodically intermittent control, Internat. J. Control (2018) http://dx.doi.org/10.1080/00207179.2018.1479538. Y. Wu, B. Chen, W. Li, Synchronization of stochastic coupled systems via feedback control based on discrete-time state observations, Nonlinear Anal.-Hybrid Syst. 26 (2017) 68–85. M. Wang, W. Li, Stability of random impulsive coupled systems on networks with Markovian switching, Stoch. Anal. Appl. (2019) http: //dx.doi.org/10.1080/07362994.2019.1643247. H. Zhou, W. Li, Synchronisation of stochastic-coupled intermittent control systems with delays and Lévy noise on networks without strong connectedness, IET Control Theory Appl. 13 (1) (2019) 36–49. H. Zhou, Y. Zhang, W. Li, Synchronization of stochastic Lévy noise systems on a multi-weights network and its applications of Chua’s circuits, IEEE Trans. Circuits Syst. I-Regul. Pap. 66 (7) (2019) 2709–2722. R.M. May, Stability and Complexity in Model Ecosystems, Princeton University Press, NJ, 2001. C. Ji, D. Jiang, D. Lei, Dynamical behavior of a one predator and two independent preys system with stochastic perturbations, Physica A 515 (2019) 649–664. Q. Liu, D. Jiang, T. Hayat, A. Alsaedi, Dynamics of a stochastic predator–prey model with stage structure for predator and Holling type II functional response, J. Nonlinear Sci. 28 (2018) 1151–1187. X. Mao, Stochastic Differential Equations, Horwood Publishing, Chichester, 1997. Q. Liu, D. Jiang, T. Hayat, A. Alsaedi, Stationary distribution and extinction of a stochastic predator–prey model with herd behavior, J. Franklin Inst. 355 (2018) 8177–8193. Q. Liu, D. Jiang, T. Hayat, A. Alsaedi, Stationary distribution of a regime-switching predator–prey model with anti-predator behaviour and higher-order perturbations, Physica A 515 (2019) 199–210. C. Lu, X. Ding, Periodic solutions and stationary distribution for a stochastic predator–prey system with impulsive perturbations, Appl. Math. Comput. 350 (2019) 313–322. W. Guo, Y. Cai, Q. Zhang, W. Wang, Stochastic persistence and stationary distribution in an SIS epidemic model with media coverage, Physica A 492 (2018) 2220–2236.

Please cite this article as: X. Zhao and Z. Zeng, Stationary distribution of a stochastic predator–prey system with stage structure for prey, Physica A (2019) 123318, https://doi.org/10.1016/j.physa.2019.123318.

12

X. Zhao and Z. Zeng / Physica A xxx (xxxx) xxx

[27] C. Ji, D. Jiang, D. Lei, Dynamical behavior of a one predator and two independent preys system with stochastic perturbations, Physica A 515 (2019) 649–664. [28] Z. Cao, W. Feng, X. Wen, L. Zu, Stationary distribution of a stochastic predator–prey model with distributed delay and higher order perturbations, Physica A 521 (2019) 467–475. [29] H. Li, F. Cong, Dynamics of a stochastic Holling-Tanner predator–prey model, Physica A 531 (2019) 121761. [30] Y. Liu, P. Yu, D. Chu, H. Su, Stationary distribution of stochastic multi-group models with dispersal and telegraph noise, Nonlinear Anal.-Hybrid Syst. 33 (2019) 93–103. [31] S. Li, H. Su, X. Ding, Synchronized stationary distribution of hybrid stochastic coupled systems with applications to coupled oscillators and a Chua’s circuits network, J. Frankl. Inst.-Eng. Appl. Math. 355 (2018) 8743–8765. [32] M. Liu, M. Deng, Analysis of a stochastic hybrid population model with allee effect, Appl. Math. Comput. 364 (2020) 124582. [33] H. Wang, M. Liu, Stationary distribution of a stochastic hybrid phytoplankton-zooplankton model with toxin-producing phytoplankton, Appl. Math. Lett. 101 (2020) 106077. [34] R.Z. Hasminskii, Stochastic Stability of Differential Equations, Sijthoff and Noordhoff, Alphen aan denRijn, 1980.

Please cite this article as: X. Zhao and Z. Zeng, Stationary distribution of a stochastic predator–prey system with stage structure for prey, Physica A (2019) 123318, https://doi.org/10.1016/j.physa.2019.123318.