Stationary distribution and extinction of a stochastic ratio-dependent predator–prey system with stage structure for the predator

Journal Pre-proof Stationary distribution and extinction of a stochastic ratio-dependent predator–prey system with stage structure for the predator Xin Zhao, Zhijun Zeng

PII: DOI: Reference:

S0378-4371(19)31854-0 https://doi.org/10.1016/j.physa.2019.123310 PHYSA 123310

To appear in:

Physica A

Received date : 16 May 2019 Revised date : 15 September 2019 Please cite this article as: X. Zhao and Z. Zeng, Stationary distribution and extinction of a stochastic ratio-dependent predator–prey system with stage structure for the predator, Physica A (2019), doi: https://doi.org/10.1016/j.physa.2019.123310. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier B.V.

Journal Pre-proof

Stationary distribution and extinction of a stochastic ratio-dependent predator-prey system with stage structure for the predator Xin Zhao, Zhijun Zeng∗

of

Northeast Normal Univ, Sch Math Stat, 5268 Renmin St, Changchun 130024, Jilin, Peoples R China.

p ro

Abstract

Pr e-

In this paper, we deal with a stochastic predator-prey model with stage structure for predator population and ratio-dependent functional response. The proposed mathematical model consists of a system of three stochastic differential equations to stimulate the interactions between prey population, immature predator and mature predator population. We first establish sufficient conditions for the existence and uniqueness of the positive solutions by constructing an appropriate Lyapunov function. Then we extend the existence of stationary distribution under certain parametric restrictions. We also obtain the sufficient conditions for extinction of the predator populations. Finally, numerical simulations have been carried out to validate our analytical findings. Keywords: Stationary distribution; Extinction; Stochastic predator-prey system; Stage structure 2000 MSC: 34D20, 39A50, 92B99

1. Introduction

Jo

urn

al

The dynamic relationship between predators and their preys has been favored by both biologists and mathematicians due to its worldwide existence and significance in population dynamics since the well-known LotkaVolterra model was brought forward [1, 2, 3, 4]. To describe a predator-prey relationship, functional response has been came up to specify the rate of prey consumption by an average predator, which has been a crucial concept in modeling predator-prey interactions [5, 6, 7]. During the past decades, the investigation of predatorprey models with different functional responses have received great attention. Traditionally, the functional response can be classified as prey-dependent (Michaelis-Menten or Holling type II [8, 9, 10, 11] ) and predatordependent for a long time. The difference lies in the functional response is determined by prey density solely or by both predator and prey populations. It has been universally acknowledged that predator-prey systems with prey-dependent functional response cannot explain the laboratory observations and have been raised doubts by several ecologists [12, 13]. In order to address these issues, Arditi and Ginzburg [14] shifted their attention to seek alternative forms for functional responses and suggested that a ratio-dependent functional response is considered as a reasonable option to model predation. The ratio-dependence is a particular type of predator dependence in which the per capita predator growth rate is a function of the ratio of prey to predator abundance [15]. This kind of interaction has been found to be typical in cases when predators have to search, share, and/or compete for food [16]. In recent years, predator-prey models with ratio-dependent functional response have revealed a wide variety of rich dynamical behaviors and are supported by numerous laboratory experiments which made the plenty of results based on ratio-dependent theory biologically more reliable. Ruan [17] studied the versal unfolding of a predator-prey system with ratio-dependent functional response and obtained all possible ∗

Corresponding author. fax: 043185098237. Email address: [email protected] (Zhijun Zeng)

Preprint submitted to Elsevier

September 15, 2019

Journal Pre-proof

Jo

urn

al

Pr e-

p ro

of

phase portraits for its perturbations. Shi [18] explored the global asymptotic stability of a diffusive predatorprey model with ratio-dependent functional response and showed that the prey and predator will be spatially homogeneously distributed as time converges to infinities. But as we all know, there are many factors that affect dynamical properties of predator-prey systems and influence the outcomes of population evolution such as the functional response, stage structure, refuge, harvesting, etc., especially the joint effect of these factors. Earlier models and approaches primarily absorbed in the simplified situations, assuming that the predators consume the preys of all ages omitting the predator eating preferences within the prey populations or ignoring the fact that only the mature predators are equipped with the ability of predation. The species represented by mammalian population whose generalized life structure is divided into two sub populations such as immature and mature. As a consequence, it is practical to introduce stage structure into the predator-prey models. The predator-prey models with stage structure for prey or predator or even for both prey and predators have been considered by several researchers [19, 20, 21, 22] since the pioneering paper of Aiello and Freedman [23]. Here we provide a ratio-dependent predator-prey system with Holling type II functional response and stage structure for the predator. Many researchers began to study based on the following system such as Deng [24]. !  dx αy2 (t)    , = x(t) r − ax(t) −    dt my2 (t) + x(t)       dy1 βx(t)y2 (t) (1)  = − by1 (t) − dy1 (t),    dt my (t) + x(t) 2      dy    2 = dy1 (t) − cy2 (t), dt where x(t) represents the population density of the prey at time t, y1 (t) and y2 (t) represent the population density of the immature and mature predator at time t, respectively. All the parameters involved in system (1) are assumed to be positive from the viewpoint of ecology, r is the intrinsic growth rate of prey, a is the intra-specific competition rate of the prey, α represents the capturing rate and αβ is the conversion rate of the mature predator, m denotes the half saturation constant, b and c are the death rate of the immature and mature predator, respectively, d is interpreted as the rate of immature predator becoming mature predator. As a matter of fact, in the real world, most natural phenomena can not be explained by deterministic laws but are always affected by the environmental noise which is an inevitable attribute of dynamics of any ecosystems ( see [25]-[30]). So the deterministic models which do not incorporate the effect of environmental noise have plenty of limitations in mathematical modelling of ecological systems and encounter difficulties in fitting data ideally and predicting the future dynamics of the system precisely. It has been shown in large number of literatures that many researchers introduced stochastic environmental variation described by Brownian motion into parameters in the deterministic model to establish stochastic predator-prey models. May [31] pointed out that the birth rates, carrying capacity, competition coefficients or other parameters involved in the model exhibit random fluctuation to a greater or lesser extent due to the environmental fluctuation. Mao [32] revealed that even a tiny environmental noise can suppress a potential population explosion. Liu [33] developed 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 and obtained sufficient conditions for extinction of the predator populations. Deng [24] focused on the Hopf bifurcation analysis for a ratio-dependent predator-prey system with two delays and stage structure for the predator and established sufficient conditions for the local stability and the existence of Hopf bifurcation with respect to both delays. Motivated by the referred works, our objective of this paper is to investigate and analyze a stochastic stage structure predator-prey model with the functional response of ratio-dependent type. There are certain methods to introduce the stochastic fluctuations into population models. Here we follow the approach adopted 2

Journal Pre-proof

in [34]-[38]. We assume the intrinsic growth rate of prey and the death rate of predators are not constants but are subject to environmental noise. Suppose that r, −b and −c are stochastically perturbed with .

r → r + σ1 w1 ,

.

− b → −b + σ2 w2 ,

.

− c → −c + σ3 w3 ,

p ro

of

where σ21 , σ22 and σ23 stand for the intensities of the white noise, w1 and w2 denote the independent standard Brownian motions. And then corresponding to the deterministic system (1), the stochastic version can be rewritten as the following form !   αy (t) 2    dt + σ1 x(t)dw1 (t), dx(t) = x(t) r − ax(t) −    my2 (t) + x(t)    !  βx(t)y2 (t) (2)    − by1 (t) − dy1 (t) dt + σ2 y1 (t)dw2 (t), dy1 (t) =    my2 (t) + x(t)      dy2 (t) = (dy1 (t) − cy2 (t))dt + σ3 y2 (t)dw3 (t). We give some basic theory in stochastic differential equations which is introduced in [25]. Let x(t) be a d-dimensional Itˆo’s process on t ≥ 0 with the stochastic differential

Pr e-

dx(t) = f (x(t), t)dt + g(x(t), t)dw(t),

where f ∈ L1 (R+ ; Rd ) and g ∈ L2 (R+ ; Rd×m ). Let V ∈ C 2,1 (Rd × R+ ; R), where C 2,1 (Rn × R+ ; R) represents for the family of all real-valued nonnegative functions V(x, t) on Rd ×R+ which are continuously twice differentiable in x and once in t. We define a diffusion operator LV : Rd × R+ → R of the Itˆo’s process associated with the C 2,1 -function V by 1 LV(x(t), t) = Vt (x(t), t) + V x (x(t), t) f (x(t), t) + trace(gT (x(t), t)V xx (x(t), t)g(x(t), t)), 2

al

then the Itˆo’s formula can be written as

dV(x(t), t) = LV(x(t), t)dt + V x (x(t), t)g(x(t), t)dw(t) a.s.

Jo

urn

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 w1 (t), w2 (t) and w3 (t) denote the independent standard Brownian motions defined on this probability space. We denote by Rd+ = {x = (x1 , · · · , xd ) ∈ Rd : xi > 0, 1 ≤ i ≤ d}. We also set inf ∅ = ∞. If A is a vector or matrix, its transpose is denoted by AT . This paper is structured in the following way. In Section 2, we demonstrate system (2) has a unique global positive solution. In Section 3, we prove system (2) has a stationary distribution under certain parametric restrictions. In Section 4, we obtain sufficient conditions for extinction of the predator populations in two cases. In Section 5, the results of numerical simulations are given to support the above theoretical results. 2. Global positive solution

In this section, we first verify system (2) has a unique global positive solution in order to investigate the dynamical behavior of system (2). Due to the coefficients of system (2) satisfy the local Lipschitz condition instead of linear growth condition, therefore, there exists a unique maximal local solution to system (2). In other words, the solution may exist from R3+ space at a finite time. We will show that the solution of system (2) is positive and global. This forms the following theorem. 3

Journal Pre-proof

Theorem 1. There is a unique global positive solution X = (x(t), y1 (t), y2(t)) of system (2) on t ≥ 0 and the solution will remain in R3+ for any given initial value (x(0), y1 (0), y2 (0)) ∈ R3+ with probability one.

of

Proof. Since the coefficients of system (2) satisfy the local Lipschitz condition, then for any initial value (x(0), y1 (0), y2 (0)) ∈ R3+ , system (2) has a unique local solution (x(t), y1 (t), y2 (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. Let n0 ≥ 0 be sufficiently large such that x(0), y1 (0) and y2 (0) are lying within the interval [ n10 , n0 ]. For each integer n ≥ n0 , we define the stopping time ) ( 1 τn = inf t ∈ [0, τe ) : min{x(t), y1 (t), y2 (t)} ≤ or max{x(t), y1 (t), y2 (t)} ≥ n . n

p ro

It’s apparently that τn is increasing as n → ∞. Set τ∞ = limn→∞ τn , whence τ∞ ≤ τe a.s. Hence, we only need to prove that τ∞ = ∞ a.s. If it is false, then there exists two constants T > 0 and δ ∈ (0, 1) such that P {τ∞ ≤ T } > δ. Therefore, there exists an integer n1 ≥ n0 such that

Pr e-

P {τn ≤ T } ≥ δ, ∀n ≥ n1 . ˜ R3+ → R+ of the form We define a C 2 -function V:

! α α x(t) ˜ + (y1 (t) − 1 − ln y1 (t)) + (y2 (t) − 1 − ln y2 (t)) . V(x(t), y1 (t), y2 (t)) = x(t) − α − α ln α β β

al

By using Itˆo’s formula, we can get ! ! ! α α 2 1 α 2 1 α α α ˜ dV = 1 − dx(t) + σ1 dt + 1− dy1 (t) + σ2 dt + 1− dy2 (t) + σ23 dt x(t) 2 β y1 (t) 2β β y2 (t) 2β α α ˜ = LV(x(t), y1 (t), y2(t))dt + (x(t) − α)σ1 dw1 (t) + (y1 (t) − 1)σ2 dw2 (t) + (y2 (t) − 1)σ3 dw3 (t), β β where

Jo

urn

! ! ! αy2 (t) 1 αx(t)y2 (t) α α ˜ LV = (x(t) − α) r − ax(t) − + 1− − by1 (t) − dy1 (t) my2 (t) + x(t) y1 (t) my2 (t) + x(t) β β ! ! 2 2 2 ασ1 ασ2 ασ3 α 1 α + 1− dy1 (t) − cy2 (t) + + + y2 (t) β β 2 2β 2β α αx(t)y2 (t) α α2 y2 (t) − by1 (t) − + (b + c + d) = rx(t) − ax2 (t) − αr + aαx + my2 (t) + x(t) β (my2 (t) + x(t))y1 (t) β 2 2 2 αdy1 (t) ασ1 ασ2 ασ3 α + + + − cy2 (t) − β βy2 (t) 2 2β 2β 2 2 2 ασ1 ασ22 ασ23 (aα + r) α α ≤ + + (b + c + d) + + + 4a m β 2 2β 2β := N, where N is a positive constant, it follows that ˜ dV(x(t), y1 (t), y2 (t)) ≤ Ndt + (x(t) − α)σ1 dw1 (t) + 4

α α (y1 (t) − 1)σ2 dw2 (t) + (y2 (t) − 1)σ3 dw3 (t). β β

Journal Pre-proof

By integrating it from 0 to τn ∧ T and taking expectation, we obtain ˜ ˜ ˜ EV(X(τ n ∧ T )) ≤ V(X(0)) + NE(τn ∧ T ) ≤ V(X(0)) + NT.

p ro

of

Let Ωn = {ω ∈ Ω : τn = τn (ω) ≤ T }, so we get P(Ωn ) ≥ δ. Consequently, for all ω ∈ Ωn , x(τn , ω), y1 (τn , ω) or y2 (τn , ω) equals either n or 1n . As a consequence, it infers that   n ˜ ˜ V(X(0)) + NT ≥ E(IΩn V(X(τ ))) ≥ δ min n − α − α ln , n − 1 − ln n , n α or ) ( 1 1 ˜ ˜ − α + α ln(nα), − 1 + ln n , V(X(0)) + NT ≥ E(IΩn V(X(τ n ))) ≥ δ min n n where IΩn denotes the indicator function of Ωn . Let n → ∞, then we see that ˜ ∞ > V(X(0)) + NT = ∞,

which leads to a contradiction. So we have τ∞ = ∞. This completes the proof.

Pr e-

3. The existence of stationary distribution

It was known that random perturbation may destroy the stability of the equilibria existing in the deterministic systems and lead to a stochastic weak stability named stationary distribution. Recently, the study of stationary distribution has attracted many researchers (see [40]-[48] and the references therein) since it is applied extensively in each realm. In this section, we shall consider there exists a unique ergodic stationary distribution of system (2) under certain conditions. Next, we introduce a lemma in order to state and prove the existence and uniqueness of a stationary distribution from Has’minskii [49].

Theorem 2. If

Jo

and

urn

al

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 \ Ω. r− βd

c+

σ23 2

σ21 α > 2 m

! σ22 > b+d+ , 2

(3)

(4)

then system (2) has a unique ergodic stationary distribution for any initial value (x(0), y1 (0), y2 (0)) ∈ R3+ . Proof. Direct computation shows that the diffusion matrix of system (2) is given by  2 2   σ1 x 0 0    σ22 y21 0  . A =  0   0 0 σ23 y22

It is clearly that matrix A is positive definite for any compact subset of R3+ . We have verified that condition A1 in Lemma 1 is satisfied. In the next moment, we define a C 2 − function V(x, y1 , y2 ) such that LV ≤ −1 on 5

Journal Pre-proof

R3+ \ Ω, where Ω is an open bounded set. This guarantees condition A2 in Lemma 1 is satisfied. Then we can show system (2) has a unique ergodic stationary distribution. Firstly, we define

of

V1 (x(t)) = x−θ (t), V2 (y1 (t), y2 (t)) = − ln y1 (t) − k ln y2 (t), !γ+1 α α(b + d) 1 x(t) + y1 (t) + y2 (t) , V3 (x(t), y1 (t), y2 (t)) = γ+1 β 2βd V4 (y1 (t)) = − ln y1 (t), V(x(t), y1 (t), y2 (t)) = V1 + MV2 + V3 + V4 ,

where θ > 0 is a sufficiently small number, M > 0 is a positive number, k =

βd

c+

γ > 0 which satisfies the following condition

σ2 3 2

!2

is also a positive number,

) b+d c . γ3 < min , 4σmax 2σmax We define λ :=

βd σ23 2

σ2 − b+d+ 2 2

and (

(5)

!

Pr e-

c+

p ro

(

γ

!γ ) a γ+2 γ3γ α α(b + d) αγ+1 (b + d) γ+1 αγ+1 c(b + d)γ+1 γ+1 γ+1 − x + σmax x + rx x + y1 + y2 − y1 − γ+3 γ+1 γ+1 y2 , B = sup 2 2 β 2βd 8βγ+1 2 β d (x,y1 ,y2 )∈R3+ n o where σmax := max σ21 , σ22, σ23 . Applying the generalized Itˆo’s formula, we arrive that θ(θ + 1) −θ−2 x (t)(dx(t))2 2 ! θ(θ + 1) −θ−2 αy2 (t) −θ−1 dt − θx−θ−1 (t)σ1 x(t)dw1 (t) + x (t)σ21 x2 (t)dt. = −θx (t)x(t) r − ax(t) − my2 (t) + x(t) 2

al

dV1 (x(t)) = −θx−θ−1 (t)dx(t) +

urn

Hence, it yields that

−θ

LV1 (x(t)) = −θx (t)

= −θx−θ (t)

Jo

≤ −θx−θ (t)

= −θx−θ (t)

! αy2 (t) θ(θ + 1) 2 −θ r − ax(t) − + σ1 x (t) my2 (t) + x(t) 2 ! θ+1 2 αθx−θ (t)y2 (t) r− σ1 + aθx1−θ (t) + 2 my2 (t) + x(t) ! α θ+1 2 σ1 + aθx1−θ (t) + θx−θ (t) r− 2 m ! θ+1 2 α + aθx1−θ (t). σ1 − r− 2 m

Similarly, it can be derived that

1 k k 1 dy1 (t) + 2 (dy1 (t))2 − dy2 (t) + 2 (dy2 (t))2 dV2 (y1 (t), y2 (t)) = − y1 (t) y2 (t) 2y1 (t) 2y2 (t) " # ! βx(t)y2 (t) y1 (t) 1 k 2 2 2 2 = − + b + d + 2 σ2 y1 (t) − kd − ck + 2 σ3 y2 (t) dt (my2 (t) + x(t))y1 (t) y2 (t) 2y1 (t) 2y2 (t) − σ2 dw2 (t) − cσ3 dw3 (t). 6

Journal Pre-proof

Thus, it follows that

=

=

≤

of

=

p ro

≤

Pr e-

LV2 (y1 (t), y2 (t)) =

! σ23 σ22 y1 βx(t)y2 (t) − +b+d+ − kd + k c + (my2 (t) + x(t))y1 (t) 2 y2 2 s ! σ23 σ22 βkdx(t) +b+d+ +k c+ −2 my2 (t) + x(t) 2 2 s ! ! p p σ23 σ22 βkdx(t) +k c+ + 2 βkd − 2 −2 βkd + b + d + 2 2 my2 (t) + x(t) s   ! !   p p σ23 σ22 x(t)  +k c+ + 2 βkd 1 − −2 βkd + b + d + 2 2 my2 (t) + x(t)  ! ! x(t) 1 − my2 (t)+x(t) p p σ23 σ22 −2 βkd + b + d + +k c+ + 2 βkd q 2 2 x(t) 1 + my2 (t)+x(t) ! ! p p σ23 σ22 my2 (t) +k c+ + 2 βkd , −2 βkd + b + d + 2 2 my2 (t) + x(t)

√ where in the second inequality, we have received the result with the help of basic inequality A + B ≥ 2 AB. If and only if A = B the equality holds. We are going to plug k = βdσ2 !2 into the above inequality, we get c+

3 2

c+

σ23 2

! σ22 2βd my2 (t) + . + b+d+ 2 σ 2 c + 23 my2 (t) + x(t)

al

LV2(y1 (t), y2 (t)) ≤ −

βd

By means of Itˆo’s formula, we obtain that

Jo

urn

!γ ! α α(b + d) α α(b + d) dV3 (x(t), y1 (t), y2 (t)) = x(t) + y1 (t) + y2 (t) dx(t) + dy1 (t) + dy2 (t) β 2βd β 2βd !γ−1 !2 α α(b + d) α α(b + d) γ y2 (t) dx(t) + dy1 (t) + dy2 (t) . + x(t) + y1 (t) + 2 β 2βd β 2βd

7

Journal Pre-proof

Then we can derive that

≤

=

of

Jo

=

urn

al

≤

p ro

≤

Pr e-

LV3(x(t), y1 (t), y2 (t)) =

!γ ! α α(b + d) α(b + d) cα(b + d) 2 x(t) + y1 (t) + y2 (t) rx(t) − ax (t) − y1 (t) − y2 (t) β 2βd 2β 2βd !γ−1 ! α α(b + d) α2 2 2 α2 (b + d)2 2 2 γ 2 2 y2 (t) σ1 x (t) + 2 σ2 y1 (t) + σ3 y2 (t) + x(t) + y1 (t) + 2 β 2βd β 4β2d 2 !γ ! α(b + d) α(b + d) cα(b + d) α 2 y2 (t) rx(t) − ax (t) − y1 (t) − y2 (t) x(t) + y1 (t) + β 2βd 2β 2βd !γ+1 γ α α(b + d) + σmax x(t) + y1 (t) + y2 (t) 2 β 2βd !γ α(b + d) αγ+1 (b + d) γ+1 α y2 (t) − axγ+2 (t) − y1 (t) rx(t) x(t) + y1 (t) + β 2βd 2βγ+1 !γ+1 cαγ+1 (b + d)γ+1 γ+1 γ α α(b + d) − γ+1 γ+1 γ+1 y2 (t) + σmax x(t) + y1 (t) + y2 (t) 2 β d 2 β 2βd αγ+1 (b + d) γ+1 cαγ+1 (b + d)γ+1 γ+1 a y (t) − y (t) − xγ+2 (t) − 1 2 4βγ+1 2γ+2 βγ+1 d γ+1 2 a cαγ+1 (b + d)γ+1 γ+1 αγ+1 (b + d) γ+1 − xγ+2 (t) − y (t) − y (t) 1 2 4βγ+1 2γ+2 βγ+1 d γ+1 2 !γ !γ+1 α α(b + d) γ α α(b + d) +rx(t) x(t) + y1 (t) + y2 (t) + σmax x(t) + y1 (t) + y2 (t) β 2βd 2 β 2βd a αγ+1 (b + d) γ+1 cαγ+1 (b + d)γ+1 γ+1 − xγ+2 (t) − y (t) − y (t) 1 2 4βγ+1 2γ+2 βγ+1 d γ+1 2 a αγ+1 (b + d) γ+1 cαγ+1 (b + d)γ+1 γ+1 − xγ+2 (t) − y (t) − y (t) 1 2 4βγ+1 2γ+2 βγ+1 d γ+1 2  !γ+1 !γ+1   γ+1  α α(b + d) γ3γ  σmax  x (t) + y1 (t) + y2 (t) + 2 β 2βd !γ α α(b + d) +rx(t) x(t) + y1 (t) + y2 (t) β 2βd αγ+1 (b + d) γ+1 cαγ+1 (b + d)γ+1 γ+1 a γ+2 γ3γ a y (t) − y (t) − x (t) + σmax xγ+1 (t) − xγ+2 (t) − 1 2 4βγ+1 2γ+2 βγ+1 d γ+1 2 2 2 !γ α α(b + d) αγ+1 (b + d) γ+1 cαγ+1 (b + d)γ+1 γ+1 +rx(t) x(t) + y1 (t) + y2 (t) − y (t) − y (t) 1 β 2βd 8βγ+1 2γ+3 βγ+1 d γ+1 2 ! ! αγ+1 (b + d) γ3γ αγ+1 σmax γ+1 cαγ+1 (b + d)γ+1 γ3γ αγ+1 (b + d)γ+1 σmax γ+1 − y1 − y2 − − 8βγ+1 2βγ+1 2γ+3 βγ+1 d γ+1 2γ+2 βγ+1 d γ+1 a αγ+1 (b + d) γ+1 cαγ+1 (b + d)γ+1 γ+1 a γ+2 γ3γ − xγ+2 (t) − y (t) − y (t) − x (t) + σmax xγ+1 (t) 1 2 4βγ+1 2γ+2 βγ+1 d γ+1 2 2 2 !γ α(b + d) αγ+1 (b + d) γ+1 cαγ+1 (b + d)γ+1 γ+1 α y2 (t) − y (t) − y (t) +rx(t) x(t) + y1 (t) + 1 β 2βd 8βγ+1 2γ+3 βγ+1 d γ+1 2 a αγ+1 (b + d) γ+1 cαγ+1 (b + d)γ+1 γ+1 − xγ+2 (t) − y (t) − y (t) + B, 1 2 4βγ+1 2γ+2 βγ+1 d γ+1 2

≤

≤

8

Journal Pre-proof

where in the last equality, we can obtain αγ+1 (b + d) γ3γ αγ+1 σmax − >0 8βγ+1 2βγ+1 and

of

cαγ+1 (b + d)γ+1 γ3γ αγ+1 (b + d)γ+1 σmax − > 0, 2γ+3 βγ+1 d γ+1 2γ+2 βγ+1 d γ+1

which follows from condition (5). According to Itˆo’s formula, it holds that

σ2 βx(t)y2 (t) + b + d + 2. (my2 (t) + x(t))y1 (t) 2

p ro

LV4 (y1 (t)) = −

Therefore, the following inequality can be estimated

Pr e-

  2 !    2Mβd my2 (t) σ θ+1 2 α  βd 2   + + aθx1−θ (t) − M  σ1 − LV ≤ −θx−θ (t) r − − b + d +  σ23 σ2 2 m 2 c+ 2 c + 23 my2 (t) + x(t) !

σ22 a γ+2 αγ+1 (b + d) γ+1 cαγ+1 (b + d)γ+1 γ+1 βx(t)y2 (t) − x (t) − y1 (t) − γ+2 γ+1 γ+1 y2 (t) + B − +b+d+ 2 4βγ+1 2 β d (my2 (t) + x(t))y1 (t) 2 ! γ+1 2Mβd my2 (t) θ+1 2 α a α (b + d) γ+1 + aθx1−θ (t) + ≤ −θx−θ (t) r − σ1 − − xγ+2 (t) − y1 (t) 2 σ 2 m 4βγ+1 c + 3 my2 (t) + x(t) 2 2

σ22 βx(t)y2 (t) cα (b + d) γ+1 y (t) − − Mλ + B + b + d + . 2γ+2 βγ+1 d γ+1 2 (my2 (t) + x(t))y1 (t) 2 γ+1

−

γ+1

Jo

urn

al

We aim to prove that LV ≤ −1 on R3+ \ Ω, once this is established, the Theorem 2 follows from Lemma 1, where Ω is an open bounded set defined as ) ( 1 2 1 1 4 3 Ω = (x, y1 , y2 ) ∈ R+ : ε < x < , ε < y1 < 4 , ε < y2 < 2 . ε ε ε S S S S S Next, we divide Ωc into six domains: Ωc = Ωc1 Ωc2 Ωc3 Ωc4 Ωc5 Ωc6 to show LV ≤ −1 actually holds, in which n o 1. Ωc1 = (x, y1 , y2 ) ∈ R3+ : 0 < x ≤ ε , o n 2. Ωc2 = (x, y1 , y2 ) ∈ R3+ : x > ε, 0 < y2 ≤ ε2 , o n 3. Ωc3 = (x, y1 , y2 ) ∈ R3+ : x > ε, 0 < y1 ≤ ε4 , y2 ≥ ε2 , n o 4. Ωc4 = (x, y1 , y2 ) ∈ R3+ : x > 1ε , o n 5. Ωc5 = (x, y1 , y2 ) ∈ R3+ : y1 > ε14 , n o 6. Ωc6 = (x, y1 , y2 ) ∈ R3+ : y2 > ε12 , and 0 < ε < 1 is a constant satisfying the following conditions   !  2     σ θ θ+1 2 α 2Mβd  a γ+2 2 1−θ +C1 ≤ −1, C1 = sup  − θ r− σ1 − x (t) + aθx (t) + − . (6) − Mλ + B + b + d +  2  σ  ε 2 m 2  2  x∈Ωc1  c + 23 9

Journal Pre-proof

2Mβdmε c+

σ23 2

+ C2 − Mλ ≤ −1,

C2 = sup

x,y2 ∈Ωc2

(

) σ22 a γ+2 1−θ . − x (t) + aθx (t) + B + b + d + 2 2

(7)

a cαγ+1 (b + d)γ+1 ε2γ+2 1 − εγ+2 − − + C3 ≤ −1, 4 2γ+2 βγ+1 d γ+1 mε3 + ε2    2     σ 2Mβd  a γ+2 2 1−θ x (t) + aθx (t) + . − Mλ + B + b + d + C3 = sup  −   σ2  2   4 x,y1 ,y2 ∈Ωc3  c+ 3 a 1 4 ε

+ C3 ≤ −1.

p ro

−

of

2

!γ+2

(8)

Pr e-

!γ+1 cαγ+1 (b + d)γ+1 1 + C1 ≤ −1. − γ+2 γ+1 γ+1 2 β d ε2 !γ+1 αγ+1 (b + d) 1 − + C1 ≤ −1. 4βγ+1 ε4 Now we have to start splitting in the following situation. I. If (x, y1 , y2 ) ∈ Ωc1 ,

(9)

(10)

(11)

! σ22 2Mβd my2 (t) θ+1 2 α a γ+2 + aθx1−θ (t) + LV ≤ −θx (t) r − σ1 − − x (t) − Mλ + B + b + d + σ2 2 m 2 c + 23 my2 (t) + x(t) 2 ! θ+1 2 α −θ + C1 σ1 − ≤ −θx (t) r − 2 m ! θ θ+1 2 α ≤− θ r− + C1 σ1 − ε 2 m ≤ −1,

urn

which follows from (6). II. If (x, y1 , y2 ) ∈ Ωc2 ,

al

−θ

LV ≤ aθx1−θ (t) +

c+

σ23 2

σ2 a my2 (t) − xγ+2 (t) − Mλ + B + b + d + 2 my2 (t) + x(t) 2 2

2Mβd mε2 + C2 − Mλ σ2 c+ 3 ε

Jo

≤

2Mβd

2

≤ −1,

which follows from (7).

10

Journal Pre-proof

III. If (x, y1 , y2 ) ∈ Ωc3 , LV ≤ aθx

(t) +

2Mβd c+

σ23 2

a γ+2 cαγ+1 (b + d)γ+1 γ+1 βx(t)y2 (t) my2 (t) − x (t) − γ+2 γ+1 γ+1 y2 (t) − my2 (t) + x(t) 2 2 β d (my2 (t) + x(t))y1 (t)

σ22 − Mλ + B + b + d + 2 a γ+2 cαγ+1 (b + d)γ+1 2γ+2 1 ≤ − ε − γ+2 γ+1 γ+1 ε − + C3 3 4 2 β d mε + ε2 ≤ −1,

LV ≤ aθx

a 1 ≤− 4 ε ≤ −1,

(t) + !γ+2

2Mβd c+

σ23 2

σ22 a γ+2 a γ+2 my2 (t) − x (t) − x (t) − Mλ + B + b + d + my2 (t) + x(t) 4 4 2

+ C3

Pr e-

1−θ

p ro

which follows from (8). IV. If (x, y1 , y2 ) ∈ Ωc4 ,

which follows from (9). V. If (x, y1 , y2 ) ∈ Ωc5 , 2Mβd σ23 2 γ+1

c+

cαγ+1 (b + d) 2γ+2 βγ+1 d γ+1 ≤ −1,

≤−

σ22 a cαγ+1 (b + d)γ+1 my2 (t) − xγ+2 (t) − γ+2 γ+1 γ+1 yγ+1 (t) − Mλ + B + b + d + 2 my2 (t) + x(t) 2 2 β d 2 !γ+1 1 + C1 ε2

al

LV ≤ aθx1−θ (t) +

of

1−θ

urn

which follows from (10). VI. If (x, y1 , y2 ) ∈ Ωc6 ,

LV ≤ aθx1−θ (t) +

2Mβd

c+

σ23 2

Jo

αγ+1 (b + d) 1 ≤− 4βγ+1 ε4 ≤ −1,

σ22 a αγ+1 (b + d) γ+1 my2 (t) − xγ+2 (t) − y − Mλ + B + b + d + 1 my2 (t) + x(t) 2 4βγ+1 2

!γ+1

+ C1

which follows from (11). Based on the analysis above-mentioned, we come to the conclusion that LV ≤ −1 on R3+ \ Ω. According to Lemma 1, we know system (2) admits a stationary distribution. This completes the proof. 4. Extinction of the predator populations Theorem 3. Let (x(t), y1 (t), y2 (t)) be the solution of system (2) with any initial value (x(0), y1 (0), y2 (0)) ∈ R3+ . If q     √ √ σ2 −2 −1 r > 21 and min{b + d, c} R − 1 I{ √R≤1} + max{b + d, c} R − 1 I{ √R>1} + cσr 1 Rr − (2(σ−2 < 0, 2 + σ3 )) 2 11

Journal Pre-proof

where R =

βd , (b+d)c

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

lim y2 (t) = 0.

t→∞

t→∞

Moreover, the distribution of x(t) converges weakly a.s. to the measure which has the density −2+ 2r2 − 2a2 x

χ(x) = Gσ−2 1 x "

σ

1

e

σ

1

, x ∈ (0, ∞),

of

 2  2r2 −1  #−1 R∞ σ σ 2r 1 1 where G = σ−2 − 1 is a constant satisfying χ(x)dx = 1. This statement highlights the Γ 2 1 2a σ 0 1

lim

t→∞

R+

f (x)ϑt (x)dx =

p ro

relationship of the probability density ϑt (x) of the prey population at time t and any sufficiently smooth function f (x) is Z Z f (x)χ(x)dx.

R+

Pr e-

Proof. Since we have obtained the solution of system (2) is positive for any initial value (x(0), y1 (0), y2 (0)) ∈ R3+ , so we get dx(t) ≤ x(t)(r − ax(t))dt + σ1 x(t)dw1 (t). We consider the following 1-dimensional stochastic differential equation with the initial value X(0) = x(0) > 0 dX(t) = X(t)(r − aX(t))dt + σ1 X(t)dw1 (t),

al

then according to the comparison theorem of stochastic differential equation one can get x(t) ≤ X(t) for any t ≥ 0 a.s. Let √ R(ω1 , ω2 ) = (ω1 , ω2 )D, √  β where (ω1 , ω2 ) = R, b+d and ! β 0 b+d D= d . 0 c

urn

We define a C 2 − function V¯ : R2+ → R+ as

¯ 1 (t), y2 (t)) = py1 (t) + qy2 (t), V(y

where p =

ω1 , b+d

q=

ω2 . c

By using Itˆo’s formula, we can derive that

Jo

1 ¯ 1 ¯ ¯ 2 d(ln V(t)) = dV(t) − (dV(t)) ¯ ¯ V(t) 2V 2(t) ! p p βx(t)y2 (t) = − by1 (t) − dy1 (t) dt + σ2 y1 (t)dw2 (t) ¯ ¯ my2 (t) + x(t) V(t) V(t) q 1 q (dy1 (t) − cy2 (t))dt + σ3 y2 (t)dw3 (t) − (p2 σ22 y21 (t) + q2 σ23 y22 (t))dt + ¯ ¯ ¯ V(t) V(t) 2V 2 (t) p q ¯ = LV(t)dt + σ2 y1 (t)dw2 (t) + σ3 y2 (t)dw3 (t), ¯ ¯ V(t) V(t) where ! p 1 βx(t)y2 (t) q ¯ LV(t) = (dy1 (t) − cy2 (t)) − (p2 σ22 y21 (t) + q2 σ23 y22 (t)). − by1 (t) − dy1 (t) + ¯ ¯ my2 (t) + x(t) V(t) V(t) 2V¯ 2(t) 12

Journal Pre-proof

In addition, one can get 1 1 V (t) = pσ2 y1 (t) + qσ3y2 (t) σ2 σ3 ¯2

!2

≤ (p

2

σ22 y21 (t)

+q

2

σ23 y22 (t))

! 1 1 + , σ22 σ23

and

Pr e-

p ro

of

( ! ) 1 βx(t)y2 (t) p − by1 (t) − dy1 (t) + q(dy1 (t) − cy2 (t)) ¯ my2 (t) + x(t) V(t) ( ! ) ! β ar βx(t) βr py2 (t) 1 p − y2 (t) − (b + d)y1 (t) + q(dy1 (t) − cy2 (t)) = + ¯ ¯ my2 (t) + x(t) my2 (t) + ar may2 (t) + r V(t) V(t)   ! ) ( my2 (t) x(t) − ar pβy2 (t) 1 βr ω2 ω1 = + y2 (t) − (b + d)y1 (t) + (dy1 (t) − cy2 (t)) ¯ ¯ (my2 (t) + x(t))(my2 (t) + ar ) V(t) b + d may2 (t) + r c V(t) r   1 pβy2 (t) X(t) − a ω2 ω1 + ≤ (βy2 (t) − (b + d)y1 (t)) + (dy1 (t) − cy2 (t)) r ¯ ¯ c V(t) V(t) b + d a   pβa r 1 X(t) − + ≤ (ω1 , ω2) D(y1 (t), y2 (t))T − (y1 (t), y2 (t))T ¯ qr a V(t)  r 1 √ pβa X(t) − + = R − 1 (ω1 y1 (t) + ω2 y2 (t)) ¯ qr a V(t) √  r 1 pβa X(t) − + R − 1 [(b + d)py1 (t) + cqy2 (t)] = qr a py1 (t) + qy2 (t) √   √ r pβa X(t) − + min{b + d, c} R − 1 I{ √R≤1} + max{b + d, c} R − 1 I{ √R>1} . ≤ qr a

Hence, we have

al

√ √ pβa X(t) − ¯ ≤ min{b + d, c}( R − 1)I{ √R≤1} + max{b + d, c}( R − 1)I{ √R>1} + LV(t) qr

r −2 −1 − (2(σ−2 2 + σ3 )) . a

urn

It follows that " √ √ pβa √ √ ¯ X(t) − d ln V(t) ≤ min{b + d, c}( R − 1)I{ R≤1} + max{b + d, c}( R − 1)I{ R>1} + qr p q + σ2 y1 (t)dw2 (t) + σ3 y2 (t)dw3 (t). ¯ ¯ V(t) V(t)

# r −2 −2 −1 − (2(σ2 + σ3 )) dt a (12)

By integrating it from 0 to t and dividing by t on both sides of (12), we get

Jo

Z √ √ ¯ − ln V(0) ¯ ln V(t) pβa t √ √ X(s) − ≤ min{b + d, c}( R − 1)I{ R≤1} + max{b + d, c}( R − 1)I{ R>1} + t qrt 0 Z Z 1 t pσ2 y1 (s) 1 t qσ3 y2 (s) −2 −2 −1 − (2(σ2 + σ3 )) + dw2 (s) + dw3 (s) ¯ ¯ t 0 V(s) t 0 V(s) Z √ √ pβa t √ √ X(s) − ≤ min{b + d, c}( R − 1)I{ R≤1} + max{b + d, c}( R − 1)I{ R>1} + qrt 0 ¯ ¯ M(t) N(t) −2 −1 − (2(σ−2 + σ )) + + , 2 3 t t 13

r ds a r ds a (13)

Journal Pre-proof

¯ where M(t) := follows

Rt

pσ2 y1 (s) dw2 (s) ¯ V(s) 0

¯ ¯ h M(t), M(t)i t =

σ22

¯ := and N(t)

Z

t 0

py1 (s) ¯ V(s)

Rt 0

!2

qσ3 y2 (s) dw3 (s) ¯ V(s)

ds ≤

are local martingales with quadratic variations as

¯ ¯ t= hN(t), N(t)i

σ22 t,

σ23

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

t→∞

0

Z t X(s) − 0

ds ≤ σ23 t.

¯ N(t) = 0 a.s. t→∞ t lim

χ(x)dx < ∞, we have 1 lim t→∞ t

!2

Z ∞ r x − r χ(x)dx ds = a a 0 # 12 "Z ∞  r 2 χ(x)dx ≤ x− a 0 !1 rσ21 2 = , 2a2

p ro

R∞

0

qy2 (s) ¯ V(s)

Pr e-

Since X(t) is ergodic and

t

of

¯ M(t) = 0, t

lim

Z

where a detailed calculation of this integral has been accomplished by Liu [33]. Taking the superior limit on both sides of (13), we have

al

1 2 !2 √ √ ¯ rσ pβa ln V(t) 1 −2 −1 − (2(σ−2 ≤ min{b + d, c}( R − 1)I{ √R≤1} + max{b + d, c}( R − 1)I{ √R>1} + lim sup 2 + σ3 )) 2 t qr 2a t→+∞ r √ √ cσ1 Rr −2 −1 √ √ ≤ min{b + d, c}( R − 1)I{ R≤1} + max{b + d, c}( R − 1)I{ R>1} + − (2(σ−2 2 + σ3 )) r 2 < 0,

which is equivalent to

lim sup It implies that

urn

t→+∞

ln y1 (t) < 0, t

lim y1 (t) = 0,

t→∞

lim sup t→+∞

ln y2 (t) < 0. t

lim y2 (t) = 0 a.s.

t→∞

This completes the proof.

Jo

Theorem 4. Let (x(t), y1 (t), y2 (t)) be the solution of system (2) with any initial value (x(0), y1 (0), y2 (0)) ∈ R3+ . If σ2 r < 21 and β < min{b, c}, then the prey and predator populations will die out, that is to say, lim x(t) = 0,

t→∞

lim y1 (t) = 0,

t→∞

lim y2 (t) = 0.

t→∞

Proof. By using Itˆo’s formula to the first equation of system (2), we can derive that ! σ21 αy2 (t) dt + σ1 dw1 (t). − d(ln x(t)) = r − ax(t) − my2 (t) + x(t) 2

14

(14)

Journal Pre-proof

By integrating it from 0 to t and dividing by t on both sides of (14), we get Z Z σ21 a t y2 (s) α t σ1 w1 (t) ln x(t) − ln x(0) x(s)ds − =r− − ds + t 2 t 0 t 0 my2 (s) + x(s) t 2 σ σ1 w1 (t) . ≤r− 1 + 2 t

(15)

lim

t→∞

w1 (t) = 0, a.s. t

lim sup t→+∞

p ro

Taking the superior limit on both sides of (15), we have

of

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

σ2 ln x(t) ≤ r − 1 < 0 a.s., t 2

which implies that

lim x(t) = 0.

t→∞

Pr e-

Making use of Itˆo’s formula to ln(y1 (t) + y2 (t)), one can get  βx(t)y2 (t)   my (t)+x(t) − by1 (t) − cy1 (t) σ2 y2 (t) + σ2 y2 (t)  3 3  dt + σ2 y1 (t) dw2 (t) + σ3 y2 (t) dw3 (t) d ln(y1 (t) + y2 (t)) =  2 − 2 2 2 y1 (t) + y2 (t) 2(y1 (t) + y2 (t))  y1 (t) + y2 (t) y1 (t) + y2 (t) ! βx(t)y2 (t) σ3 y2 (t) σ2 y1 (t) ≤ − min{b, c} dt + dw2 (t) + dw3 (t). (my2 (t) + x(t))(y1 (t) + y2 (t)) y1 (t) + y2 (t) y1 (t) + y2 (t)

urn

al

According to lim x(t) = 0, we know there exists T > 0 such that x(t) < ǫ when t > T . Hence, a direct calculation t→∞ shows that ! βǫ σ3 y2 (t) σ2 y1 (t) d ln(y1 (t) + y2 (t)) ≤ − min{b, c} dt + dw2 (t) + dw3 (t) my2 (t) + ǫ y1 (t) + y2 (t) y1 (t) + y2 (t) (16) σ3 y2 (t) σ2 y1 (t) dw2 (t) + dw3 (t). ≤ (β − min{b, c}) dt + y1 (t) + y2 (t) y1 (t) + y2 (t) Rt R t σ y (s) 2 y1 (s) 3 2 Noting that M(t) := 0 y1σ(s)+y dw (s) and N(t) := dw3 (s) are local martingales with quadratic varia2 (s) 0 y1 (s)+y2 (s) 2 tions as follows !2 !2 Z t Z t y2 (s) y1 (s) 2 2 2 ds ≤ σ2 t, hN(t), N(t)it = σ3 ds ≤ σ23 t. hM(t), M(t)it = σ2 y1 (s) + y2 (s) y1 (s) + y2 (s) 0 0

Jo

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

t→∞

M(t) = 0, t

N(t) = 0 a.s. t→∞ t lim

By integrating it from 0 to t and dividing by t on both sides of (16), we get Z Z 1 t σ2 y1 (s) 1 t σ3 y2 (s) ln(y1 (t) + y2 (t)) − ln(y1 (0) + y2 (0)) ≤ β − min{b, c} + dw2 (s) + dw3 (s) t t 0 y1 (s) + y2 (s) t 0 y1 (s) + y2 (s) M(t) N(t) + . ≤ β − min{b, c} + t t (17) 15

Journal Pre-proof

Taking the superior limit on both sides of (17), we have lim sup t→+∞

ln(y1 (t) + y2 (t)) ≤ β − min{b, c} < 0 a.s., t

which implies that lim y1 (t) = 0,

t→∞

lim y2 (t) = 0 a.s.

t→∞

of

This completes the proof. 5. Numerical test

Pr e-

p ro

In this section, we will carry out some numerical simulations to illustrate our theoretical results. We shall corroborate all the analytical findings with the help of numerical simulations accomplished with MATLAB. The initial value of system (2) is given by (0.2, 0.2, 0.2). We use Milstein’s higher-order method mentioned in Higham [50] to obtain the following discretization transformation of system (2) !  i √  αy σ21 i 2  2 i i i+1 i i   ∆t + σ1 x ∆tε1,i + x (ε1,i − 1)∆t, x = x + x r − ax − i    2 my2 + xi     !   √  i+1 βxi yi2 σ22 i 2 i i i i (18)  y1 (ε2,i − 1)∆t, ∆tε + y = y + − by − dy ∆t + σ y 2,i 2 1  1 1 1 1 i  i  2 my + x  2      √ σ23 i 2   i i i i   yi+1 y (ε − 1)∆t, = y + (dy − cy )∆t + σ y ∆tε + 3 2 3,i 2 2 1 2 2 2 3,i

where the time increment ∆t is positive and εk,i are the Gaussian random variables which follow the distribution N(0, 1), k = 1, 2, 3.

find r = 0.6 > 0.5 =

σ21 2

+

α m

and

al

Example 1. We numerically simulate the model for the following set of parameter values: σ21 = 0.04, σ22 = 0.34, σ23 = 0.08. Other values of the system parameters can be seen from Table 1. By simple computations, we βd

σ2 c+ 23

= 4 > 1 = b+d+

σ22 . 2

It is obvious that this group of parameters satisfy

urn

condition (3) and (4). From Theorem 2, there exists an ergodic distribution ν(·) of system (2) which is clearly shown in Figure 1. Table 1: List of parameters

Description Intrinsic growth rate of the prey Intra-specific competition rate of the prey Death rate of immature predator Death rate of mature predator The rate of immature predator becoming mature predator Capture rate of the predator Conversion rate of the mature predator Half saturation constant

Jo

Parameters r a b c d α β m

16

Values 0.6 0.2 0.03 0.04 0.8 0.6 0.4 1.25

Journal Pre-proof

4

8

5 4

6 x(t)

x 10

3 4 2 2 0

1 0

1000 2000 3000 4000 5000 Time t

0

0

1

2

3

4

4

4

10

2

5

x 10

0

1000

2000 3000 Time t

4000

5000

0

0

0.5

1

1.5

2

2.5

p ro

0

of

15

y1(t)

6

y2(t)

5

20

4

15

3

10

2

5

1

0

1000

2000 3000 Time t

4000

5000

0

0

100

200

300

Pr e-

0

x 10

Figure 1: The left column shows the paths of x(t), y1 (t) and y2 (t) of system (2) with initial value (x(0), y1(0), y2 (0)) = (0.2, 0.2, 0.2) under the noise intensities σ21 = 0.04, σ22 = 0.34 and σ23 = 0.08. The blue lines represent the solution to system (2) and the red lines represent the solution to the corresponding undisturbed system (1). The right column displays the histogram of the probability density functions of x, y1 , y2 populations.

urn

al

Example 2. We numerically simulate the model for the following set of parameter values: σ21 = 0.04, σ22 = 15, σ23 = 15. Other values of the system parameters can be seen from Table 2. By simple computations, we find q √ σ21 βd cσ1 Rr √ r = 0.5 > 0.02 = 2 , R = (b+d)c = 16 > 1, max{b + d, c} = 1 and max{b + d, c}( R − 1)I{ R>1} + r − 2 −2 −2 −1 (2(σ2 + σ3 )) = 3.008 − 3.75 < 0. Obviously this group of parameters obey Theorem 3, y1 (t) and y2 (t) will become extinct as t → ∞ with probability 1. This is consistent with the results in Theorem 3. It also explains the big noise can cause the extinction of species which is clearly shown in Figure 2. Table 2: List of parameters

Description Intrinsic growth rate of the prey Intra-specific competition rate of the prey Death rate of immature predator Death rate of mature predator The rate of immature predator becoming mature predator Capture rate of the predator Conversion rate of the mature predator Half saturation constant

Jo

Parameters r a b c d α β m

17

Values 0.5 0.2 0.6 0.01 0.4 0.6 0.4 1.25

Journal Pre-proof

8

6

6 y1(t)

x(t)

4 4

2 2 0

0

1000

2000 3000 Time t

4000

5000

0

1000

2000 3000 Time t

4000

5000

0

0

1000 2000 3000 4000 5000 Time t

of

20

y2(t)

15 10

Pr e-

0

p ro

5

Figure 2: The paths of x(t), y1 (t) and y2 (t) of system (2) with initial value (x(0), y1(0), y2 (0)) = (0.2, 0.2, 0.2) under the noise intensities σ21 = 0.04, σ22 = 15 and σ23 = 15. The blue lines represent the solution to system (2) and the red lines represent the solution to the corresponding undisturbed system (1).

Example 3. We numerically simulate the model for the following set of parameter values: σ21 = 2, σ22 = 15, σ23 = 15. Other values of the system parameters can be seen from Table 3. By simple computations, we σ2

urn

al

find r = 0.5 < 1 = 21 and β = 0.4 < 0.5 = min{b, c}. It is obvious that this group of parameters satisfy the conditions of Theorem 4. It indicates that the big noise can cause both species of system (2) to die out ultimately and it may happen when a serious disease or severe weather occurs in the real life. This is clearly shown in Figure 3. Table 3: List of parameters

Description Intrinsic growth rate of the prey Intra-specific competition rate of the prey Death rate of immature predator Death rate of mature predator The rate of immature predator becoming mature predator Capture rate of the predator Conversion rate of the mature predator Half saturation constant

Jo

Parameters r a b c d α β m

18

Values 0.5 0.2 0.6 0.5 0.7 0.6 0.4 1.25

Journal Pre-proof

8

6

6 y1(t)

x(t)

4 4

2 2 0

0

1000 2000 3000 4000 5000 Time t

0

0

1000 2000 3000 4000 5000 Time t

of

20

y2(t)

15 10

0

1000 2000 3000 4000 5000 Time t

Pr e-

0

p ro

5

Figure 3: The paths of x(t), y1 (t) and y2 (t) of system (2) with initial value (x(0), y1(0), y2 (0)) = (0.2, 0.2, 0.2) under the noise intensities σ21 = 2, σ22 = 15 and σ23 = 15. The blue lines represent the solution to system (2) and the red lines represent the solution to the corresponding undisturbed system (1).

6. Acknowledgment

urn

References

al

This work was supported by the Fundamental Research Funds for the Central Universities (Grant No. 2412017FZ004).

[1] A.A. Berryman, The origin and evolution of predator-prey theory, Ecology 73 (1992) 1530-1535. [2] A.J. Lotka, Elements of Physical Biology, William and Wilkins, Baltimore (1925), Reissued as Elements of Mathematical Biology, Dover, New York, 1956.

Jo

[3] V. Volterra, Fluctuations in the abundance of a species considered mathematically, Nature 118 (1926) 558-560. [4] D. Jiang, C. Ji, X. Li, D. O’Regan, Analysis of autonomous Lotka-Volterra competition systems with random perturbation, J. Math. Anal. Appl. 390 (2) (2012) 582-595. [5] Y. Cai, C. Zhao, W. Wang, Dynamics of a Leslie-Gower predator-prey model with additive Allee effect, Appl. Math. Model. 39 (7) (2015) 2092-2106. [6] R.M. May, Stability and Complexity in Model Ecosystems, Princeton Univ. Press, Princeton, 1974. [7] J.D. Murray, Mathematical Biology, Springer-Verlag, Berlin, 1989. 19

Journal Pre-proof

[8] H.I. Freedman, Deterministic Mathematical Models in Population Ecology, Monogr. Textbooks Pure Appl. Math., vol. 57, Marcel Dekker, New York, 1980. [9] C.S. Holling, The functional response of predator to prey density and its role in mimicry and population regulation, Mem. Entomol. Soc. Can. 97 (1965) 5-60.

of

[10] 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. [11] M. Liu, Dynamics of a stochastic regime-switching predator-prey model with modified Leslie-Gower Holling-type II schemes and prey harvesting, Nonlinear Dyn. 96 (2019) 417-442.

p ro

[12] C.B. Huffaker, Experimental studies on predation: dispersion factors and predator-prey oscillations, Hilgardia 27 (1958) 343-383. [13] R. Arditi, L.R. Ginzburg, H.R. Akcakaya, Variation in plankton densities among lakes: a case for ratiodependent models, Am. Nat. 138 (1991) 1287-1296.

Pr e-

[14] R. Arditi, L.R. Ginzburg, Coupling in predator-prey dynamics: ratio-dependence, J. Theor. Biol. 139 (1989) 311-326. [15] S.B. Hsu, T.W. Hwang, Y. Kuang, Global analysis of the Michelis-Menten type ratio-dependent predatorprey system, J. Math. Biol. 42 (2001) 489-506. [16] A.P. Gutierrez, The physiological basis of ratio-dependent predator-prey theory: a metabolic pool model of Nicholson’s blowflies as an example, Ecology 73 (1992) 1552-1563. [17] S. Ruan, Y. Tang, W. Zhang, Versal unfoldings of predator-prey systems with ratio-dependent functional response, J. Differential Equations 249 (2010) 1410-1435.

al

[18] H. Shi, Y. Li, Global asymptotic stability of a diffusive predator-prey model with ratio-dependent functional response, Appl. Math. Comput. 250 (2015) 71-77.

urn

[19] R. Xu, M.A.J. Chaplain, F.A. Davidson, Persistence and global stability of a ratio-dependent predator-prey model with stage structure, App. Math. Comput. 158 (2004) 729-744. [20] S. Khajanchi, Modeling the dynamics of stage-structure predator-prey system with Monod-Haldane type functional response, Appl. Math. Comput. 302 (2017) 122-143.

Jo

[21] M. Liu, K. Wang, Global stability of stage-structure predator-prey models with Beddington-Deangelis functional response, Commun. Nonlinear Sci. Numer. Simul. 16 (2011) 3792-3797. [22] X. Shi, X. Zhou, X. Song, Analysis of a stage-structured predator prey model with Crowley-Martin function, J. Appl. Math. Comput. 36 (2011) 459-472. [23] W.G. Aiello, H.I. Freedman, A time-delay model of single-growth with stage structure, Math. Biosci. 101 (1990) 139-153. [24] L. Deng, X. Wang, M. Peng, Hopf bifurcation analysis for a ratio-dependent predator-prey system with two delays and stage structure for the predator, Appl. Math. Comput. 231 (2014) 214-230. [25] X. Mao, Stochastic Differential Equations, Chichester, Horwood Publishing, 1997. 20

Journal Pre-proof

[26] M. Wang, W. Li, Stability of random impulsive coupled systems on networks with Markovian switching, Stoch. Anal. Appl. DOI: 10.1080/07362994.2019.1643247 (2019). [27] H. Zhou, Y. Zhang, W. Li, Synchronization of stochastic L´evy 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.

of

[28] H. Zhou, W. Li, Synchronisation of stochastic-coupled intermittent control systems with delays and L´evy noise on networks without strong connectedness, IET Contr. Theory Appl. 13 (1) (2019) 36-49.

p ro

[29] S. Li, B. Zhang, W. Li, Stabilisation of multiweights stochastic complex networks with time-varying delay driven by G-Brownian motion via aperiodically intermittent adaptive control, Int. J. Control DOI: 10.1080/00207179.2019.1577562 (2019). [30] P. Wang, B. Zhang, H. Su, Stabilization of stochastic uncertain complex-valued delayed networks via aperiodically intermittent nonlinear control, IEEE Trans. Syst. Man Cybern. -Syst. 49 (3) (2019) 649-662. [31] R.M. May, Stability and complexity in model ecosystems, Princeton University Press, NJ, 2001.

Pr e-

[32] X. Mao, G. Marion, E. Renshaw, Environmental Brownian noise suppresses explosions in population dynamics, Stoch. Process. Appl. 97 (2002) 95-110. [33] 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. [34] C. Ji, D. Jiang, N. Shi, A note on a predator-prey model with modified Leslie-Gower and Holling-type II schemes with stochastic perturbation, J. Math. Anal. Appl. 377 (2011) 435-440. [35] X. Mao, Stationary distribution of stochastic population systems, Syst. Control Lett. 60 (2011) 398-405.

al

[36] Y. Wu, J. Zhu, W. Li, Intermittent discrete observation control for synchronization of stochastic neural networks, IEEE T. Cybern. DOI: 10.1109/TCYB.2019.2930579 (2019).

urn

[37] M. Liu, M. Fan, Stability in distribution of a three-species stochastic cascade predator-prey system with time delays, IMA J. Appl. Math. 82 (2017) 396-423. [38] M. Liu, J. Yu, P.S. Mandal, Dynamics of a stochastic delay competitive model with harvesting and Markovian switching, Appl. Math. Comput. 337 (2018) 335-349. [39] Y. Wu, C. Wang, W. Li, Generalized quantized intermittent control with adaptive strategy on finite-time synchronization of delayed coupled systems and applications, Nonlinear Dyn. 95 (2) (2019) 1361-1377.

Jo

[40] Q. Liu, D. Jiang, Stationary distribution and extinction of a stochastic predator-prey model with distributed delay, Appl. Math. Lett. 78 (2018) 79-87. [41] 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. [42] 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. [43] D. Jiang, W. Zuo, Stationary distribution and periodic solutions for stochastic Holling-Leslie predator-prey systems, Physica A 460 (2016) 16-18. 21

Journal Pre-proof

[44] Y. Cai, X. Mao, Stochastic prey-predator system with foraging arena scheme, Appl. Math. Model. 64 (2018) 357-371. [45] M. Liu, Y. Zhu, Stationary distribution and ergodicity of a stochastic hybrid competition model with L´evy jumps, Nonlinear Anal. Hybrid Syst. 30 (2018) 225-239.

of

[46] 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.

p ro

[47] 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. [48] Y. Liu, H. Xu, W. Li, Intermittent control to stationary distribution and exponential stability for hybrid multi-stochastic-weight coupled networks based on aperiodicity, J. Franklin Inst. 356 (13) (2019) 72637289.

Pr e-

[49] R.Z. Hasminskii, Stochastic stability of differential equations. Sijthoff and Noordhoff, Alphen aan denRijn (1980).

Jo

urn

al

[50] D.J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev. 43 (2001) 525-546.

22

Journal Pre-proof

A stochastic ratio-dependent predator-prey system with stage structure for the predator is studied.

of

We establish sufficient conditions for the existence of a unique

p ro

ergodic stationary distribution.

We obtain sufficient conditions for extinction which is discussed

Jo

urn

al

Pr e-

in two cases.