Conditions for persistence and ergodicity of a stochastic Lotka–Volterra predator–prey model with regime switching

Commun Nonlinear Sci Numer Simulat 29 (2015) 1–11

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Conditions for persistence and ergodicity of a stochastic Lotka–Volterra predator–prey model with regime switching Li Zu a, Daqing Jiang b,∗, Donal O’Regan c,d a

College of Mathematics and Statistics, Hainan Normal University, Hainan 571158, PR China School of Science, China University of Petroleum (East China), Qingdao 266580, PR China c School of Mathematics, Statistics and Applied Mathematics, National University of Ireland, Galway, Ireland d Nonlinear Analysis and Applied Mathematics (NAAM)-Research Group, King Abdulaziz University, Jeddah, Saudi Arabia b

a r t i c l e

i n f o

Article history: Received 17 October 2014 Revised 18 February 2015 Accepted 2 April 2015 Available online 27 April 2015 Keywords: Stochastic Lotka–Volterra predator–prey model Markov switching Ergodic and positive recurrent persistence in mean

a b s t r a c t A stochastic Lotka–Volterra predator–prey model with regime switching is investigated. We examine when the system is persistence in mean and the extinction under some appropriate conditions, and we discuss the threshold between them. Sufficient conditions for the stationary distribution which is ergodic and positive recurrent of the solution is established using stochastic Lyapunov functions. Simulations are also carried out to illustrate our analytical results. © 2015 Elsevier B.V. All rights reserved.

1. Introduction The dynamic relationship between predators and their preys has long been and will continue to be one of the dominant themes in both ecology and mathematical ecology [22]. In addition to mathematics and biology, the application of predator–prey models are also reflected in other areas, for example, physics. Recently, predator–prey interaction have been studied in structured populations in cyclical interactions with alliance-specific heterogeneous invasion rates [34] and in noise-guided evolution within cyclical interactions [35]. Reviews of modern physics pointed out that the potential applicability of the proposed theory extends to statistical mechanics of evolutionary and coevolutionary games, in collective behavior and evolutionary games [36], see [37,38] for more related theories. We focus on the developments and changes of predator–prey models, a vast wealth of literature has investigated on them. From the early simple predator–prey models to the later Gauss-type [29], Leslie–Gower model [25], ratio-dependent predator–prey system [23,39], Holling type [24], Beddington–DeAngelis functional response model [26,27], etc., and based on these models, many scholars have studied the non-autonomous and periodic coefficients model, the models with impulsive perturbations [40] or with time delays [30] or on time scales [28]. In addition, taking into account that nature is full of uncertainty and random phenomena, differential equations disturbed by environment white noise, colored noise have been studied [31–33,41], recently, and the existence and uniqueness of solutions, persistence, stationary distribution and ergodicity, global attractivity and extinction of stochastic systems have been discussed [9–17,20].

∗

Corresponding author. Tel.: +86 53286983361. E-mail address: [email protected], [email protected] (D. Jiang).

http://dx.doi.org/10.1016/j.cnsns.2015.04.008 1007-5704/© 2015 Elsevier B.V. All rights reserved.

2

L. Zu et al. / Commun Nonlinear Sci Numer Simulat 29 (2015) 1–11

From all of the above models, it can be seen that Lotka–Volterra predator–prey models occupy an important position. The classical Lotka–Volterra predator–prey model with white noise is expressed by

 dx(t) = x(t)[a1 − b1 y(t) − c1 x(t)]dt + α x(t)dB1 (t), dy(t) = y(t)[ − a2 + b2 x(t) − c2 y(t)]dt + β y(t)dB2 (t).

(1.1)

The stochastic processes x(t) and y(t) represent, respectively, the prey and predator population. Now ai (i = 1, 2) denotes the intrinsic growth rate of the corresponding population, the coefficient b1 is the capturing rate of the predator, b2 stands for the rate of conversion of nutrients into the reproduction of the predator, and ci represents the density-dependent coefficients of the prey and the predator, respectively. Here Bi (t) is one-dimensional standard Brownian motion, α 2 and β 2 stand for the intensities of the white noises on the prey and predator population, respectively. Note all these coefficients are positive constants. Chessa et al. [4] studied the system (1.1) and obtained the existence, uniqueness and non-extinction property of the solution, whereas Rudnicki [2] proved the system (1.1) has a unique distribution density using the theory of Markov semigroups assuming that the random noise B1 (t) and B2 (t) are correlated. System (1.1) was studied extensively (see [1,3,5,6] and the references therein). Besides environment white noise, in this paper we will also consider another colored noise, say telegraph noise [7,8]. Telegraph noise can be described as a random switching between two or more environmental regimes, which differ in terms of factors such as nutrition or rainfall. Let {r(t), t  0} be a right-continuous Markov chain on the probability space (, F , {Ft }t≥0 , P), taking values in a finite-state space S = {1, . . . , m}. The generator  = (θ ij )1i, jm is given, for δ > 0, by



P {r(t + δ) = j|r(t) = i} =

θij δ + o(δ), if i = j, 1 + θij δ + o(δ), if i = j.

Here θ ij is the transition rate from i to j and θ ij  0 if i  j, while θ ii = −ij θ ij . In this paper we assume θ ij > 0 if i  j. Suppose that the Markov chain r(t) is independent of the Brownian motion B(·) and it is irreducible. Under this condition, the Markov chain has a unique stationary (probability) distribution π = (π 1 , . . . , π m )  R1×m , which can be determined by solving the linear m π = 1 and π > 0, ∀k ∈ S. We refer the reader to [9–11] for details. equation π  = 0, k=1 k k The population stochastic system (1.1) with regime switching can be described by the following model



dx(t) = x(t)[a1 (r(t)) − b1 (r(t))y(t) − c1 (r(t))x(t)]dt + α(r(t))x(t)dB1 (t), dy(t) = y(t)[−a2 (r(t)) + b2 (r(t))x(t) − c2 (r(t))y(t)]dt + β(r(t))y(t)dB2 (t).

(1.2)

Assume, for any k ∈ S, that the coefficients a1 (k), a2 (k), b1 (k), b2 (k), c1 (k), c2 (k), α (k) and β (k) are positive. Here Bi (t)(i = 1, 2) are independent standard Brownian motions. The main aim of this paper is to investigate the persistence in mean and the extinction of system (1.2) and find the threshold between them, which is important for assessing the risk of extinction of a population exposed to an environmental toxicant. By constructing suitable Lyapunov functions, we give sufficient conditions for system (1.2) being positive recurrent (and the existence of a unique ergodic stationary distribution). Liu and Wang [12] investigated a logistic model driven by Brownian noise and colored noise and obtained the threshold between weak persistence and extinction. There are only a few results on the model (1.2) and other stochastic predator–prey systems under regime switching. Stationary distribution of stochastic population models were studied in [15–17] and Zhu and Yin [18] studied the ergodicity of a stochastic model with Markov switching. Using the theory in [18], Liu et al. [13] and Settati and Lahrouz [14], investigated a stochastic Lotka–Volterra mutualistic system with regime switching and obtained the existence of the stationary distribution on the stochastic mutualistic system. In this paper we study the ergodic stationary distribution of a stochastic Lotka–Volterra predator–prey model with regime switching using the theory and methods of [18]. The paper is organized as follows. We recall the fundamental theory and give two theorems concerning the existence of a global positive solution and moment estimation in Section 2. Section 3 investigates the persistence in mean and the extinction of the population model (1.2) and we find the threshold between them. We show the existence of a stationary distribution in Section 4. In Section 5, we present an example and figures to illustrate the main results. Finally we present some conclusions and future directions in Section 6. 2. Preliminaries In this paper, we suppose that there is a complete probability space (, F , {Ft }t≥0 , P) with a filtration {Ft }t≥0 satisfying the usual conditions (i.e. it is right continuous and F0 contains all P-null sets). Let R2+ denote the positive zone of R2 , namely t  R2+ = {(x, y) ∈ R2 : x > 0, y > 0}. For any vector φ = (φ (1), . . . , φ (m))T , limt→∞ 1t 0 φ(r(s))ds = k∈S πk φ(k). Let φˆ = mink∈S {φ(k)} and φˇ = maxk∈S {φ(k)}. Lemma 2.1 ([12]). Suppose that z(t)  C[ × [0, +), R+ ]. If there are positive constants λ0 , T and λ  0 such that

log z(t) ≥ λt − λ0



t 0

z(s)ds +

n 

βi Bi (t)

i=1

for t  T, where β i is a constant, 1  i  n, then lim inft→∞

1 t

t 0

z(s)ds ≥ λλ , a.s. 0

L. Zu et al. / Commun Nonlinear Sci Numer Simulat 29 (2015) 1–11

3

Now, we introduce a lemma which gives a criterion for positive recurrence in terms of Lyapunov functions (see e.g. [18], Theorem 3.13, [13], Lemma 2.1). First, consider the diffusion process (z(t), r(t)) described by the equation:

 dz(t) = f (z(t), r(t))dt + g(z(t), r(t))dB(t), z(0) = z, r(0) = r,

(2.1)

where B(·) and r(·) are the d-dimensional Brownian motion and the right-continuous Markov chain in the above discussion, respectively, and f (·, ·) : Rn × S → Rn , g(·, ·) : Rn × S → Rn×d satisfying g(z, k)gT (z, k) = D(z, k). For each k ∈ S, and for any twice continuously differentiable function V(·, k), (z(t), r(t)) has a generator L given as follows LV (z, k) =

n 

fi (z, k)

i=1

n ∂ V (z, k) 1  ∂ 2 V (z, k) + Dij (z, k) + V (z, ·)(k), 2 ∂ zi ∂ zi ∂ zj i,j=1

where

V (z, ·)(k) =

n  i=1

θki V (z, i) =



θki (V (z, i) − V (z, k)), k ∈ S.

i=k∈S

Lemma 2.2. If the following conditions are satisfied: (i) for i  j, θ ij > 0, (ii) for each k ∈ S, D(z, k) is symmetric and satisfies

|ζ |2 ≤ ζ T D(z, k)ζ ≤ −1 |ζ |2 for all ζ ∈ Rn , with some constant ϱ  (0, 1] for all z  Rn , (iii) there exists a bounded open subset D of Rn with a regular (i.e. smooth) boundary satisfying that, for each k ∈ S there exists a nonnegative function V (·, k) : D C → R such that V(·, k) is twice continuously differentiable and that for some ς > 0, LV (z, k) ≤ −ς , for any (z, k) ∈ D C × S, then system (2.1) is ergodic and positive recurrent. That is to say, there exists a unique stationary   density μ(·, ·) and, for any Borel measurable function f (·, ·) : Rn × S → R such that k∈S Rn |f (z, k)|μ(z, k)dz < ∞, we have    1 t P( limt→∞ t 0 f (z(s), r(s))ds = k∈S Rn f (z, k)μ(z, k)dz) = 1. The following two theorems concern the existence and uniqueness of positive solutions and moment boundedness. Since the proof of these two theorems are standard [[10,17], see e.g.], we omit the proofs here. Theorem 2.1. For any given initial value (x(0), y(0), r(0)) in R2+ × S, there is a unique solution to system (1.2) on t  0 and the solution will remain in R2+ × S with probability 1. Theorem 2.2. For any given initial value (x(0), y(0), r(0)) ∈ R2+ × S and for any p > 0, there exists a constant ι(p) > 0, such that the solution (x(t), y(t), r(t)) of system (1.2) satisfies E[x(t) + y(t)]p  ι(p) for all t  0. 3. The threshold between persistence in mean and the extinction In the study of population dynamics, we are interested in whether a population can sustain development or become extinct. Our aim in this section is to study the persistence in mean and the extinction and give the persistence in mean-extinction threshold of system (1.2). This is one of our main conclusions of the paper. Definition 3.1. • •

System (1.2) is said to be extinctive almost surely (a.s.), if limt →  y(t) = 0 a.s.; t System (1.2) is said to be persistence in mean if lim inft→∞ 1t 0 y(s)ds > 0 a.s.

Remark 3.1. Our definition of persistence in mean is taken from [42]. Liu and Wang [20] introduced the definitions that if t t lim supt→∞ 1t 0 y(s)ds > 0 a.s., then y(t) is stochastic weak persistence in the mean, and if lim inft→∞ 1t 0 y(s)ds > 0 a.s., then y(t) is stochastic strong persistence in the mean. Obviously, the condition of stochastic strong persistence in the mean is stronger than stochastic weak persistence in the mean. The definition of persistence in mean in this paper and the definition of stochastic strong persistence in the mean in [20] are equivalent. We now consider the system on the boundary

 d(t) = (t) a1 (r(t)) − c1 (r(t))(t) dt + α(r(t))(t)dB1 (t).

(3.1)

By a comparison theorem, we can check that x(t)  (t) t  0 a.s. provided that x(0) = (0) > 0 and y(0) > 0. If

 k∈S

πk [a1 (k) −

α 2 (k) < 0, lim t →  (t) = 0 a.s. is easy to obtain from Theorem 4.1 of [10]. Hence, limt →  x(t) = 0 a.s., which implies 2 ]  2 lim y(t) = 0 a.s. For this reason, we suppose π (k) − α (k) > 0 throughout the paper. t→

k∈S

k [a1

2

]

4

L. Zu et al. / Commun Nonlinear Sci Numer Simulat 29 (2015) 1–11

Also combined with the results of [14], we know that if distribution ν (·, ·) which is ergodic, and then we have



c1 (k)



R+

k∈S

Let λ =

 k∈S

xν(dx, k) =

b2 (k)





πk a1 (k) −

α 2 (k)

 R+

xν(dx, k) −

 k∈S

k∈S

πk [a1 (k) −

α 2 (k) > 0, system (3.1) has a unique stationary 2 ]



a.s.

2

k∈S



β 2 (k)

πk [a2 (k) +

2

(3.2)

] which is the persistence in mean-extinction threshold we want.

 2 Theorem 3.1. Assume k∈S πk [a1 (k) − α 2(k) ] > 0 hold. If λ < 0, then for any initial value (x(0), y(0), r(0)) ∈ R2+ × S, the predator populations goes to extinction almost surely. Proof. Taking the logarithm of both sides of the second equation of (1.2) and using the generalized Itô’s formula, we get

log y(t) − log y(0) t 

1 t =− a2 (r(s)) + t 0 

1 t ≤− a2 (r(s)) + t 0 where M1 (t) =

β 2 (r(s)) 2

β 2 (r(s)) 2

t

 ds +

1 t



t

b2 (r(s))x(s)ds −

0



1 t

t

c2 (r(s))y(s)ds +

0

  1 t M1 (t) ds + , b2 (r(s))(s)ds + t 0 t

M1 (t) t (3.3)

β(r(s))dB2 (s) is a real-valued continuous local martingale (see e.g. [12] on page 16) whose quadratic variation ˇ 2 t. Making use of the strong law of large numbers for martingales yields limt→∞ M1 (t) = 0 a.s. is M1 , M1 t = 0 β 2 (r(s))ds ≤ (β) t Taking the superior limit on both sides of inequality (3.3) leads to lim sup t→∞

 0t

log y(t)  b2 (k) ≤ t



xν(dx, k) −

R+

k∈S





πk a2 (k) +

β 2 (k)



=λ

2

k∈S

a.s.

(3.4)

Obviously, if λ < 0, then y(t) tends to zero a.s. The proof is complete.   2 Theorem 3.2. Assume k∈S πk [a1 (k) − α 2(k) ] > 0 hold. If λ > 0, then system (1.2) will be persistent in mean, for any initial value (x(0), y(0), r(0)) ∈ R2+ × S. Proof. Applying the generalized Itô’s formula, we get from the first equation of (1.2) that

1 log x(t) − log x(0) = t t

where M2 (t) =

0≥

t 0

 t

0

1 − t



a1 (r(s)) −

t

0

α 2 (r(s)) 2

c1 (r(s))x(s)ds +

  1 t ds − b1 (r(s))y(s)ds t 0

M2 (t) , t

(3.5)

α(r(s))dB1 (s). It then follows from (3.5) that

1 log x(t) − log (t) = t t ≥ cˆ1



t 0

1 t

1 t  t

c1 (r(s))[(s) − x(s)]ds −



t

0

[(s) − x(s)]ds − bˇ 1

1 t

0

 0

t

b1 (r(s))y(s)ds

y(s)ds,

which implies

1 t



t 0

[(s) − x(s)]ds ≤

Now, we find

1 t

t 0

 bˇ 1 1 t y(s)ds. cˆ1 t 0

(3.6)

[(s) − x(s)]ds from system (1.2). Combining (3.3) and (3.6), we calculate

log y(t) − log y(0) t 

1 t =− a2 (r(s)) + t 0 

1 t ≥− a2 (r(s)) + t 0  t

1 ≥− a2 (r(s)) + t 0

β 2 (r(s)) 2

 ds +

1 t

 0

t

b2 (r(s))x(s)ds −

1 t

 0

t

c2 (r(s))y(s)ds +

M1 (t) t

    M1 (t) 1 t 1 t 1 t b2 (r(s))(s)ds − bˇ 2 y(s)ds + ds + [(s) − x(s)]ds − cˇ2 2 t 0 t 0 t 0 t

   t t 2 ˇ ˇ M1 (t) β (r(s)) 1 1 b1 b2 . b2 (r(s))(s)ds − + cˇ2 y(s)ds + ds + 2 t 0 t 0 t cˆ1

β 2 (r(s))

M (t)

(3.7)

In addition, according to the property of inferior limits and the fact limt→∞ 1t = 0 a.s., we know that for arbitrary ε > 0, there t t  2 2   exists t > 0 such that 1t 0 b2 (r(s))(s)ds − 1t 0 [a2 (r(s)) + β (2r(s)) ]ds ≥ k∈S b2 (k) R+ xν(dx, k) − k∈S πk [a2 (k) + β 2(k) ] − ε2 =

L. Zu et al. / Commun Nonlinear Sci Numer Simulat 29 (2015) 1–11

λ−

5

ε and M1 (t) ≥ − ε . Therefore inequality (3.7) yields t 2 2

log y(t) − log y(0) ≥λ−ε− t



bˇ 1 bˇ 2 + cˇ2 cˆ1

1 t



t

0

y(s)ds.

(3.8)

Then, by Lemma 2.1 and letting ε → 0, yields

lim inf t→∞

1 t



t 0

λ

y(s)ds ≥  ˇ

b1 bˇ 2 cˆ1

+ cˇ2

 > 0 a.s.

(3.9)

Therefore the desired assertion follows from (3.9) immediately.   2 Remark 3.2. Assume the condition k∈S πk [a1 (k) − α 2(k) ] > 0 holds. If λ < 0, the predator will go to extinction and the prey has a stationary distribution, and if λ > 0, the system (1.2) will be persistent in mean. λ is the persistence in mean-extinction threshold. If λ  0, then λ > 0 is the sufficient and necessary condition for system (1.2) persistence in mean. 4. Ergodic property of positive recurrence The aim of this section is to investigate the ergodic property of system (1.2) by using the Lyapunov function method. To begin with, let us impose the condition

(H) λ¯ = bˆ2





πk a1 (k) −

α 2 (k) 2

k∈S

 − cˇ1





πk a2 (k) +

β 2 (k) 2

k∈S

 > 0.

Recall that by putting ξ (t) = log x(t) and η(t) = log y(t) for t  0, we have

⎧

⎪ ⎪ ⎪ ξ ( t ) = a1 (r(t)) − d ⎨ ⎪ ⎪ ⎪ ⎩dη(t) =



α 2 (r(t))

− a2 (r(t)) −

2

 − b1 (r(t))eη(t) − c1 (r(t))eξ (t) dt + α(r(t))dB1 (t),

β 2 (r(t)) 2

 + b2 (r(t))eξ (t) − c2 (r(t))eη(t) dt + β(r(t))dB2 (t).

(4.1)

From the proof of Lemma 3.2 in [13], we know that the ergodic property of system (1.2) is equivalent to that of system (4.1). The following theorem is to verify that system (4.1) satisfies the three conditions of Lemma 2.2. Theorem 4.1. Let hypothesis (H) hold. Then for any k ∈ S and for any initial value (x(0), y(0), r(0)) ∈ R2+ × S be given, (x(t), y(t), r(t)) of system (1.2) is ergodic and has a unique stationary distribution in R2+ × S. Proof. We assume θ ij > 0, i  j in Section 1 and then condition (i) in Lemma 2.2 holds. To verify condition (ii), consider the bounded open subset D = {(ξ , η) : |ξ | ≤ log  −1 , |η| ≤ log  −1 , (ξ , η) ∈ R2 }, where 0 < ε < 1 is a sufficiently small number. Let g(z, k) = diag(α(k), β(k))(k ∈ S), then D(z, k) = g(z, k)gT (z, k) = diag(α 2 (k), β 2 (k)) is positive definite. Therefore, condition (ii) in Lemma 2.2 is satisfied. We now need to verify condition (iii) in Lemma 2.2. First, define a C2 -function

[eξ + peη ]2 bˇ 1 bˆ 2 + cˇ1 cˇ2 η h(ξ , η) = M[−bˆ 2 ξ − cˇ1 η + e ]+ , 2 aˆ 2 here, p =

bˆ 1 bˇ 2

¯ max {2, sup and M = (2/λ) (ξ ,η)∈R2 {−

cˆ1 3ξ 2 e

−

p2 cˆ2 3η 2 e

(4.2)

ˇ2 2 + (aˇ 1 + αˇ2 )e2ξ − p2 (aˆ 2 − β2 )e2η + pˇa1 eξ +η }}. From the partial

derivative equations of h(ξ , η), we know that

−

Mbˆ 2 pˇc1 eξ + eξ + =0 ξ ˇ ˆ c1 cˇ2 ) ξ e pbˆ 2 + (b1 b2aˆ+ˇ e 2

has a unique solution ξ 0 , which can be seen from the monotonically property of the left function. The following equation

pbˆ 2 cˇ1 bˇ 1 bˆ 2 + cˇ1 cˇ2 − ξ =− η + e aˆ 2 e has a unique solution (ξ 0 , η0 ) which is the minimum point of h(ξ , η), here η0 = log h(ξ 0 , η0 )  0.



cˇ1 bˇ 1 bˆ 2 +ˇc1 cˇ2 aˆ 2

+pbˆ 2 e−ξ0



. Therefore we have h(ξ , η) −

6

L. Zu et al. / Commun Nonlinear Sci Numer Simulat 29 (2015) 1–11

Define a C2 -function V : R2 × S → R+ by



 cˇ1 cˇ2 + bˇ 1 bˆ 2 η [eξ + peη ]2 ˆ − h(ξ0 , η0 ) + M(k + | |) V (ξ , η, k) = M −b2 ξ − cˇ1 η + e + 2 aˆ 2 = V1 (ξ , η) + V2 (ξ , η) − h(ξ0 , η0 ) + V3 (k), 

)T , | |

(4.3)

where  = (1 , . . . , m = and k (k ∈ S) will be determined below. We put |ϖ| here in order to make ϖk + |ϖ| non-negative. Therefore, from the above analysis, we know that V(ξ , η, k) is non-negative. An application of the Itô’s formula, yields



12



LV1 (ξ , η) ≤ M −bˆ 2 a1 (k) − and LV3 (k) = M



+ · · · + m2

α 2 (k) 2



+ cˇ1 a2 (k) +

β 2 (k) 2



+

bˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 ) ξ +η e aˆ 2

(4.4)

θkl (l − k ).

(4.5)

l=k∈S

We now define a vector  = (1 , . . . , m )T with k = bˆ 2 [a1 (k) − α 2(k) ] − cˇ1 [a2 (k) + β 2(k) ]. Since the generator matrix  is irreducible, then for k , there exists ϖ = (ϖ1 , . . . , ϖm )T a solution of the poisson system (see [19], Lemma 2.3]), such that 2

2

k − k = −

m 

πj j .

(4.6)

j=1

Thus, we have







θkl (l − k ) − bˆ2 a1 (k) −

l=k∈S





= − bˆ 2 πk a1 (k) −

α 2 (k)

k∈S

2

α 2 (k) 2

 − cˇ1



− cˇ1 a2 (k) +





πk a2 (k) +

k∈S

β 2 (k)



2

β 2 (k)



2

¯. = −λ ¯ + Consequently, from (4.4), (4.5) and (4.7), one gets L(V1 + V3 ) ≤ M(−λ

(4.7) bˇ 2 (bˇ 1 bˆ 2 +ˇc1 cˇ2 ) ξ +η e ). aˆ 2

We note

LV2 (ξ , η)

  α 2 (k) β 2 (k) 2η e2ξ + p2 e = (eξ + peη ) a1 (k)eξ − b1 (k)eξ +η − c1 (k)e2ξ − pa2 (k)eη + pb2 (k)eξ +η − pc2 (k)e2η + 2 2

 2 2 αˇ βˇ e2ξ − p2 aˆ 2 − e2η + pˇa1 eξ +η + pbˇ 2 e2ξ +η + p2 bˇ 2 eξ +2η − bˆ 1 e2ξ +η −pbˆ 1 eξ +2η ≤ −ˆc1 e3ξ − p2 cˆ2 e3η + aˇ 1 + 2 2

  αˇ 2 2ξ βˇ 2 2η 3ξ 2 3η 2 e + pˇa1 eξ +η . = −ˆc1 e − p cˆ2 e + aˇ 1 + e − p aˆ 2 − 2 2

It therefore follows that

  bˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 ) ξ +η αˇ 2 2ξ βˇ 2 2η 3ξ 2 3η 2 ¯ LV (ξ , η, k) ≤ −Mλ + M e − p aˆ 2 − e + pˇa1 eξ +η . e − cˆ1 e − p cˆ2 e + aˇ 1 + 2 2 aˆ 2

In the set D C × S, we choose sufficiently small ε such that

0<<

¯ aˆ 2 λ , 4bˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 )

(4.8)

0<<

p2 aˆ 2 cˆ2 , 2Mbˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 )

(4.9)

0<<

aˆ 2 cˆ1 , ˇ ˇ 2Mb2 (b1 bˆ 2 + cˇ1 cˇ2 )

(4.10)

¯ − −Mλ

cˆ1 + K1 ≤ −1, 2 3

(4.11)

¯ − −Mλ

p2 cˆ2 + K2 ≤ −1, 2 3

(4.12)

L. Zu et al. / Commun Nonlinear Sci Numer Simulat 29 (2015) 1–11

7

where K1 and K2 are positive constants which can be found from (4.13) and (4.14). Denote that D1 = {(ξ , η) ∈ R2 : −∞ < ξ ≤ log }, D3 = {(ξ , η) ∈ R2 : ξ ≥ log  −1 },

D2 = {(ξ , η) ∈ R2 : −∞ < η ≤ log }, D4 = {(ξ , η) ∈ R2 : η ≥ log  −1 }.

Obviously, D C = D1 ∪ D2 ∪ D3 ∪ D4 . Next, we will prove LV (ξ , η, k) ≤ −1 on D C , which is equivalent to proving it on the above four domains. Case 1. For any (ξ , η, k) ∈ D1 × S, since eξ + η  ε eη  ε (1 + e3η ), we have



¯ ¯ Mλ Mbˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 ) Mbˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 ) Mλ p2 cˆ2 cˆ1 3ξ + − + + LV (ξ , η, k) ≤ − − e + − e3η 4 4 2 2 aˆ 2 aˆ 2

    ¯ cˆ1 3ξ Mλ p2 cˆ2 3η αˇ 2 2ξ βˇ 2 2η 2 ξ +η − e − e + aˇ 1 + e + pˇa1 e + − e − p aˆ 2 − 2 2 2 2 2

¯ ¯ Mλ Mbˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 ) Mbˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 ) cˆ1 p2 cˆ2 Mλ + − + + − e3ξ + − e3η ≤− 4 4 2 2 aˆ 2 aˆ 2

     ¯ Mλ p2 cˆ2 3η cˆ1 3ξ αˇ 2 2ξ βˇ 2 2η 2 ξ +η + sup − e − e + aˇ 1 + e − p aˆ 2 − e + pˇa1 e + − . 2 2 2 2 2 (ξ ,η)∈R2

¯ max {2, sup By the definition of M = (2/λ) (ξ ,η)∈R2 {− ¯ Mλ 4

cˆ1 3ξ 2 e

−

p2 cˆ2 3η 2 e

ˇ2 2 + (aˇ 1 + αˇ2 )e2ξ − p2 (aˆ 2 − β2 )e2η + pˇa1 eξ +η }}, we know that

≥ 1. From (4.8) and (4.9), we obtain LV (ξ , η, k) ≤ −

¯ ¯ cˆ1 Mλ Mλ − e3ξ ≤ − ≤ −1. 4 2 4

Therefore, LV ≤ −1 on D1 × S. Case 2. On D2 × S, since eξ + η  ε eξ  ε (1 + e3ξ ), we have





¯ ¯ Mλ Mbˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 ) Mbˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 ) Mλ p2 cˆ2 3η cˆ1 + − + e + − + e3ξ − 4 4 2 2 aˆ 2 aˆ 2

     ¯ cˆ1 Mλ p2 cˆ2 3η αˇ 2 2ξ βˇ 2 2η + sup − e3ξ − e + aˇ 1 + e + pˇa1 eξ +η . + − e − p2 aˆ 2 − 2 2 2 2 2 (ξ ,η)∈R2

LV (ξ , η, k) ≤ −

Combining (4.8) and (4.10), we get LV (ξ , η, k) ≤ −

¯ ¯ p2 cˆ2 3η Mλ Mλ − e ≤− ≤ −1. 4 2 4

Obviously, LV ≤ −1 for all (ξ , η, k) ∈ D2 × S. Case 3. On D3 × S, LV (ξ , η, k)



   cˆ1 3ξ cˆ1 3ξ bˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 ) ξ +η αˇ 2 2ξ βˇ 2 2η 2 3η 2 ¯ e − p aˆ 2 − e + pˇa1 + M ≤ −Mλ − e + − e − p cˆ2 e + aˇ 1 + e 2 2 2 2 aˆ 2   2 ¯ − cˆ1 e3ξ + sup − cˆ1 e3ξ − p2 cˆ2 e3η + aˇ 1 + αˇ e2ξ ≤ −Mλ 2 2 2 (ξ ,η)∈R2



 βˇ 2 2η bˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 ) ξ +η 2 e + pˇa1 + M e − p aˆ 2 − 2 aˆ 2 ¯ − cˆ1 + K1 , ≤ −Mλ 2 3

where K1 = sup(ξ ,η)∈R2 {−

cˆ1 3ξ 2 e

(4.13) ˇ2 2 bˇ (bˇ bˆ +ˇc cˇ ) − p2 cˆ2 e3η + (aˇ 1 + αˇ2 )e2ξ − p2 (aˆ 2 − β2 )e2η + (pˇa1 + M 2 1 aˆ2 1 2 )eξ +η }, which gives LV ≤ −1

in this domain, in view of (4.11).

2

8

L. Zu et al. / Commun Nonlinear Sci Numer Simulat 29 (2015) 1–11

Case 4. When (ξ , η, k) ∈ D4 × S,

  αˇ 2 2ξ p2 cˆ2 3η p2 cˆ2 3η e + −ˆc1 e3ξ − e + aˇ 1 + e 2 2 2

 βˇ 2 2η bˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 ) ξ +η e + (pˇa1 + M −p2 aˆ 2 − )e 2 aˆ 2   2 2 2 ¯ − p cˆ2 e3η + sup −ˆc1 e3ξ − p cˆ2 e3η + aˇ 1 + αˇ ≤ −Mλ e2ξ 2 2 2 (ξ ,η)∈R2



 bˇ 2 (bˇ 1 bˆ 2 + cˇ1 cˇ2 ) ξ +η βˇ 2 2η e + pˇa1 + M − p2 aˆ 2 − e 2 aˆ 2

¯ − LV (ξ , η, k) ≤ −Mλ

¯ − ≤ −Mλ

p2 cˆ2 + K2 , 2 3

(4.14)

p2 cˆ2 3η 2 e

ˇ2 2 bˇ (bˇ bˆ +ˇc cˇ ) + (aˇ 1 + αˇ2 )e2ξ − p2 (aˆ 2 − β2 )e2η + (pˇa1 + M 2 1 aˆ2 1 2 )eξ +η }. Note from (4.12), we

where K2 = sup(ξ ,η)∈R2 {−ˆc1 e3ξ −

2

have LV ≤ −1 on D4 × S. Therefore, we deduce that LV (ξ , η, k) ≤ −1,

for all (ξ , η, k) ∈ DC × S.

(4.15)

Therefore according to Lemma 2.2, we know that (ξ (t), η(t), r(t)) is ergodic and positive recurrent, and so the system (1.2) is positive recurrent and admits a unique stationary distribution by Lemma 3.2 in [13]. The stationary density is in R2+ × S which can be established by a proof similar to that in Lemma 3.1 in [13].  Remark 4.1. (a) We prove system (4.1) is positive recurrent, by verifying the three conditions of Lemma 2.2, thereby proving system (1.2) is positive recurrent and has a unique stationary distribution by using Lemma 3.2 of [13]. Here we use the equivalent theory of system (1.2) and system (4.1). (b) The Lyapunov function V(ξ , η, k) consists of three parts. In order to guarantee V(ξ , η, k) is positive, we analyze the ¯ in LV and λ ¯ reflects the impact of environmenfunction h(ξ , η) and introduce |ϖ|. In addition, with the help of ϖ(k), we get the λ tal white noise and the probability distribution of the Markov Chain. (c) The condition λ¯ > 0 is stronger than λ > 0. From Theorem 4.1, we know that system (1.2) has a unique stationary distribution in R2+ × S, expressed as μ(·, ·). We get the following theorem. Theorem 4.2. Let us assume hypothesis (H) holds. For any k ∈ S, we have



b1 (k)

 R+

k∈S

yμ(dy, k) +



c1 (k)

 R+

k∈S

xμ(dx, k) =

 k∈S

and



b2 (k)



k∈S

R+

xμ(dx, k) −



c2 (k)



k∈S

R+



πk a1 (k) −

yμ(dy, k) =





πk a2 (k) +

k∈S

α 2 (k)

 (4.16)

2

β 2 (k) 2

 .

(4.17)

Proof. Theorem 2.2 states that, for any p > 0, there exists a positive constant ι(p), such that E[x(t) + y(t)]p  ι(p), and then we have E[x(t)]  ι(1), and similarly, E[y(t)]  ι(1). Then, by the ergodic property of (x(t), y(t), r(t)) and for every integer n > 0, we have



1 P lim t→∞ t



t 0

(x(s) ∧ n)ds =

 k∈S

R+

(x ∧ n)μ(x, k)dx = 1.

Using the dominated convergence theorem, we get

     1 t 1 t E lim (x(s) ∧ n)ds = lim E(x(s) ∧ n)ds = (x ∧ n)μ(x, k)dx. t→∞ t 0 t→∞ t 0 R+ k∈S

It is easy to see that

 k∈S

R+

xμ(x, k)dx = lim

n→∞

 k∈S

R+

(x ∧ n)μ(x, k)dx ≤ lim sup E[x(t)] < +∞, t→∞

which implies the function f(x) = x is integrable with respect to the measure μ. Similarly, one can obtain

(4.18)

L. Zu et al. / Commun Nonlinear Sci Numer Simulat 29 (2015) 1–11

(a)

9

(b)

16

10 k=1 k=2

14

k=1 k=2

9 8

12 7 10

6

8

5 4

6

3 4 2 2 0

1 0

0.5 1 1.5 2 The density functions of x(t)

0

0

0.5 1 1.5 2 The density functions of y(t)

Fig. 1. The curves on the subgraphs (a) and (b) are the density functions of x(t) and y(t) of stochastic system (1.2) in k = 1, k = 2, respectively. The initial value is x(0) = 0.9, y(0) = 0.7. (α (1), α (2)) = (0.07, 0.09), (β (1), β (2)) = (0.1, 0.05). The density functions of x(t) and y(t) concentrated in smaller intervals. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

 R+

k∈S

yμ(y, k)dy ≤ lim sup E[y(t)] < +∞. t→∞

Now, let us return to (3.5). By the ergodicity of (x(t), y(t), r(t)), one can see that limt→∞ logtx(t) = 0 a.s. Otherwise, if limt→∞ logtx(t) > 0 or < 0, we conclude that, limt→ x(t) = + or limt→ x(t) = 0, which contradicts the fact that its stationary density lies in R+ . Thus, we have shown that





πk a1 (k) −

k∈S

α 2 (k) 2



−



b1 (k)

k∈S



R+

yμ(dy, k) −

 k∈S

c1 (k)



R+

xμ(dx, k) = 0,

(4.19)

which proves the first part (4.16). The second result also holds in a similar way.  5. Examples and computer simulations In this section, we will illustrate our analytical results by an example with the help of the method mentioned in [21] and Matlab software. Example 5.1. For the sake of convenience we assume that the system (1.2) switches from one to the other according to the movement of the Markov chain r(t) on the state space S = {1, 2} with the following coefficients: when k = 1,

(a1 (1), a2 (1)) = (6, 4), (b1 (1), b2 (1)) = (6, 5), (c1 (1), c2 (1)) = (4, 3), and when k = 2,

(a1 (2), a2 (2)) = (7, 4.5), (b1 (2), b2 (2)) = (4.5, 4), (c1 (2), c2 (2)) = (5, 4). The Markov chain has the generator

=



 −2 2 . 1 −1

Then its stationary distribution is π = (π1 , π2 ) = ( 13 , 23 ). Next we will discuss two cases according to the intensity of white noise. Case 1. Let (α (1), α (2)) = (0.07, 0.09), (β (1), β (2)) = (0.1, 0.05). Relatively small white noise act on predator and prey. 2   β 2 (k) α 2 (k) α 2 (k) ¯ ˆ Note k∈S πk [a1 (k) − 2 ] − cˇ1 k∈S πk [a2 (k) + k=1 πk [a1 (k) − 2 ] = 26.6526 > 0, λ = b2 2 ] = 4.9734 > 0. Then, the conditions of Theorem 4.1 hold, and system (1.2) has a unique stationary distribution. We now use Figs. 1–2 to illustrate the results. Fig. 1 describes the density function images of the stationary distribution of x(t) and y(t) in k = 1 and k = 2, respectively. The red , blue ◦ and black + represent the phase portrait of x(t) and y(t) when there is only one state k = 1, k = 2 and switching

10

L. Zu et al. / Commun Nonlinear Sci Numer Simulat 29 (2015) 1–11

(a)

(b)

2

2 Alternative k=1 k=2

1.6

1.6

1.4

1.4

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

0.5

1 x

1.5

Alternative k=1 k=2

1.8

y

y

1.8

0

2

0

0.5

1 x

1.5

2

Fig. 2. The subgraphs (a) and (b) denote sample phase portrait of stochastic system (1.2) and the corresponding deterministic system with initial value x(0) = 0.9, y(0) = 0.7. The red , blue ◦ and black + area denote x(t), y(t) in k = 1, k = 2 and converting between the above two states, respectively. (α (1), α (2)) = (0.07, 0.09), (β (1), β (2)) = (0.1, 0.05). The black area located between the areas of k = 1 and k = 2. The system (1.2) has a stationary distribution. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

(a)

(b)

2

2 Alternative k=1 k=2

1.6

1.6

1.4

1.4

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

0.5

1 x

1.5

Alternative k=1 k=2

1.8

y

y

1.8

2

0

0

0.5

1 x

1.5

2

Fig. 3. The subgraphs (a) and (b) have the same definitions as in Fig. 2. (α (1), α (2)) = (1.1, 1), (β (1), β (2)) = (1.6, 1.5). Large white noise can lead to population extinction, though the deterministic system is persistent.

back and forth from one state k = 1 to another state k = 2 according to the movement of r(t), respectively, in Fig. 2. In Fig. 2(a), the system is affected by white noise and the corresponding subgraph Fig. 2(b) has no random disturbance. Clearly, the black area is located between the red region and the blue region, and the red area and the blue area is similar to the two limit state of the black region. That is to say, the switching state is located between the states k = 1 and k = 2, and thus it can be seen that the existence of colored noise reduce the risk of species extinction. Case 2. Let (α (1), α (2)) = (1.1, 1), (β (1), β (2)) = (1.6, 1.5). Predator and prey are affected by relatively big white noise. 2

πk [a1 (k) − α 2(k) ] = 6.1317 > 0 and λ¯ = −3.0233 < 0. The hypothesis (H) does not hold. From Remark 4.1 (c), we ¯ = −3.0233 < 0, that is, the predator will become extinct by Theorem 3.1, and Fig. 3(a) confirms this. Large white know that λ < λ noise will lead to population extinction, even though the corresponding deterministic model is persistent (see Fig. 3(b)). Note

2

k=1

6. Conclusions and future directions This article discusses a basic predator–prey Lotka–Volterra model disturbed by both white noise and colored noise. We obtain sufficient conditions for persistence in mean and the extinction of system (1.2) and get the threshold between them. We prove

L. Zu et al. / Commun Nonlinear Sci Numer Simulat 29 (2015) 1–11

11

the stochastic system under regime switching has a unique stationary distribution which is ergodic and positive recurrent by using the Lyapunov function method. The main difference between the deterministic and stochastic model is that large noise can also cause extinction. The existence of a Markov chain in this stochastic system can reduce the risk of species becoming extinct which can be seen from the example and figures. In future studies, we will try to develop this approach to other predator–prey models. Acknowledgements The work was supported by Program for NSFC of China (No: 11371085), the Fundamental Research Funds for the Central Universities (No:15CX08011A ), and the Foundation of Hainan Provincial Science and Technology Department (No: 20151009). 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] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42]

Mandal PS, Banerjee M. Stochastic persistence and stationary distribution in a Holling–Tanner type prey–predator model. Phys A 2012;391:1216–33. Rudnicki R. Long-time behaviour of a stochastic prey–predator model. Stochastic Process Appl 2003;108:93–107. Rudnicki R, Pichór K. Influence of stochastic perturbation on prey–predator systems. Math Biosci 2007;206:108–19. Chessa S, Fujita Yashima H. The stochastic equation of predator–prey population dynamics (Italian). Boll Unione Mat Ital Sez B Artic Ric Mat 2002;5(8):789– 804. Wu RH, Zou XL, Wang K. Asymptotic behavior of a stochastic non-autonomous predator–prey model with impulsive perturbations. Commun Nonlinear Sci Numer Simul 2015;20:965–74. Yagi A, Ton TV. Dynamic of a stochastic predator–prey population. Appl Math Comput 2011;218(31):3100–9. Mao XR, Sabanis S, Renshaw E. Asymptotic behaviour of the stochastic Lotka–Volterra model. J Math Anal Appl 2003;287:141–56. Li XY, Jiang DQ, Mao XR. Population dynamical behavior of Lotka–Volterra system under regime switching. J Comput Appl Math 2009;232:427–48. Liu M, Wang K. Asymptotic properties and simulations of a stochastic logistic model under regime switching II. Math Comput Model 2012;55:405–18. Li XY, Gray A, Jiang DQ, Mao XR. Sufficient and necessary conditions of stochastic permanence and extinction for stochastic logistic populations under regime switching. J Math Anal Appl 2011;376:11–28. Mao XR. Stochastic differential equations and applications. Horwood; 1997. Liu M, Wang K. Asymptotic properties and simulations of a stochastic logistic model under regime switching. Math Comput Model 2011;54:2139–54. Liu H, Li XX, QS Yang. The ergodic property and positive recurrence of a multi-group Lotka–Volterra mutualistic system with regime switching. Syst Control Lett 2013;62:805–10. Settati A, Lahrouz A. Stationary distribution of stochastic population systems under regime switching. Appl Math Comput 2014;244:235–43. Ji CY, Jiang DQ, Shi NZ. A note on a predator–prey model with modified Leslie–Gower and Holling-type II schemes with stochastic perturbation. J Math Anal Appl 2011;377:435–40. Mao XR. Stationary distribution of stochastic population systems. Syst Control Lett 2011;60:398–405. Zu L, Jiang DQ, O’Regan D. Stochastic permanence, stationary distribution and extinction of a single-species nonlinear diffusion system with random perturbation. Abstr Appl Anal 2012;2012:1–14. Zhu C, Yin G. Asymptotic properties of hybrid diffusion systems. SIAM J Control Optim 2007;46:1155–79. Khasminskii RZ, Zhu C, Yin G. Stability of regime-switching diffusions. Stochastic Process Appl 2007;117:1037–51. Liu M, Wang K. Persistence and extinction of a stochastic single-specie model under regime switching in a polluted environment. J Theor Biol 2010;264:934– 44. Higham DJ. An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Rev 2001;43:525–46. Berryman AA. The origins and evolution of predator–prey theory. Ecology 1992;73:1530–5. Arditi R, Ginzburg LR. Coupling in predator–prey dynamics: ratio-dependence. J Theor Biol 1989;139:311–26. Holling CS. The functional response of predators to prey density and its role in mimicry and population regulation. Mem Entomol Soc Canada 1965;46:1–60. Hsu SB, Huang TW. Global stability for a class of predator–prey systems. SIAM J Appl Math 1995;55:763–83. Beddington JR. Mutual interference between parasites or predators and its effect on searching efficiency. J Anim Ecol 1975;44:331–40. DeAngelis DL, Goldstein RA, ONeill RV. A model for trophic interaction. Ecology 1975;56:881–92. Zhang JM, Fan M, Zhu HP. Periodic solution of single population models on time scales. Math Comput Model 2010;52:515–21. Liu GR, Yan JR. Existence of positive periodic solutions for neutral delay Gause-type predator–prey system. Appl Math Model 2011;35:5741–50. Sahoo B, Poria S. Effects of additional food in a delayed predator–prey model. Math Biosci 2015;261:62–73. Zu L, Lin XN, Jiang DQ. Existence theory for single and multiple solutions to singular boundary value problems for second order impulsive differential equations. Topol Method Nonlinear Anal 2007;30(1):171–91. ´ noise. J Math Anal Appl 2012;391:363–75. Bao JH, Yuan CG. Stochastic population dynamics driven by Levy CY Ji, Jiang DQ, Shi NZ. Analysis of a predator–prey model with modified Leslie–Gower and Holling-type II schemes with stochastic perturbation. J Math Anal Appl 2009;359:482–98. Perc M, Szolnoki A, Szab G. Cyclical interactions with alliance-specific heterogeneous invasion rates. Phys Rev E 2007;75:052102. Perc M, Szolnoki A. Noise-guided evolution within cyclical interactions. New J Phys 2007;9:267. Perc M, Grigolini P. Collective behavior and evolutionary games—an introduction. Chaos Solitons Fractals 2013;56:1–5. Szolnoki A, Mobilia M, Jiang LL, Szczesny B, Rucklidge AM, Perc M. Cyclic dominance in evolutionary games: a review. J R Soc Interface 2014;11:20140735. Castellano C, Fortunato S, Loreto V. Statistical physics of social dynamics. Rev Mod Phys 2009;81(2):591–646. Shi HB, Li Y. Global asymptotic stability of a diffusive predator–prey model with ratio-dependent functional response. Appl Math Comput 2015;250:71–7. Wu RH, Zou XL, Wang K. Asymptotic behavior of a stochastic non-autonomous predator–prey model with impulsive perturbations. Commun Nonlinear Sci Numer Simul 2015;20:965–74. Zhang XH, Li WX, Liu M, Wang K. Dynamics of a stochastic Holling II one-predator two-prey system with jumps. Phys A 2015;421:571–82. Chen L, Chen J. Nonlinear biological dynamical system. Beijing: Science Press; 1993.