Stability of a stochastic one-predator-two-prey population model with time delays

Accepted Manuscript

Stability of a stochastic one-predator-two-prey population model with time delays Jing Geng, Meng Liu, Youqiang Zhang PII: DOI: Reference:

S1007-5704(17)30128-4 10.1016/j.cnsns.2017.04.022 CNSNS 4172

To appear in:

Communications in Nonlinear Science and Numerical Simulation

Received date: Revised date: Accepted date:

25 October 2016 6 March 2017 23 April 2017

Please cite this article as: Jing Geng, Meng Liu, Youqiang Zhang, Stability of a stochastic one-predatortwo-prey population model with time delays, Communications in Nonlinear Science and Numerical Simulation (2017), doi: 10.1016/j.cnsns.2017.04.022

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.

ACCEPTED MANUSCRIPT

Highlights • A stochastic one-predator-two-prey model with delays is proposed. • An asymptotic approach is used to study the stability in distributions of the model.

AC

CE

PT

ED

M

AN US

CR IP T

• Sharp sufficient criteria for the stability in distributions are established.

1

ACCEPTED MANUSCRIPT

Stability of a stochastic one-predator-two-prey population model with time delays Jing Genga , Meng Liua,∗, Youqiang Zhangb of Mathematical Science, Huaiyin Normal University, Huaian 223300, PR China of Mathematics, Harbin Institute of Technology (Weihai), Weihai, 264209, PR China

CR IP T

a School b Department

Abstract

This paper is concerned with the stability in distribution of a delay stochastic population model with two competing preys (X1 and X2 ) and one predator (X3 ). Under some assumptions we prove that there are three numbers γ1 > γ2 > γ3 which have the following properties: if γ1 < 1, then all the populations go to extinction lim Xi (t) = 0 a.s., i = 1, 2, 3; If γi > 1 > γi+1 , i = 1, 2, then the distribution

t→+∞

AN US

almost surely (a.s.), i.e.,

of (X1 (t), ..., Xi (t))T converges weakly to a unique ergodic invariant distribution and lim Xj (t) = 0 a.s., j = t→+∞

i+1, ..., 3; If γ3 > 1, then the distribution of (X1 (t), X2 (t), X3 (t))T converges weakly to a unique ergodic invariant distribution a.s.. The influence of random perturbations on the stability are discussed and some numerical simulations are given to illustrate the main results.

M

Keywords: Two-prey-one-predator model; random perturbations; delays; stability in distribution.

1. Introduction

ED

In the natural world, it is a common phenomenon that a predator feeds on some competing preys, for example, Hyenas, Zebra and Connochaetes in Africa. On the other hand, time delays should be taken into account in biological models ([1]). Therefore in the past few decades, delay population models with two competing preys

PT

and one predator have received great attention and have been studied extensively. Among various types of delay population models with two competing preys and one predator, we should specially mention the following

AC

CE

Lotka-Volterra system, which has been widely studied in [2]-[6]:      dX1 (t) = X1 (t) b1 − c11 X1 (t) − c12 X2 (t − τ12 ) − c13 X3 (t − τ13 ) dt,        dX2 (t) = X2 (t) b2 − c21 X1 (t − τ21 ) − c22 X2 (t) − c23 X3 (t − τ23 ) dt,          dX3 (t) = X3 (t) − b3 + c31 X1 (t − τ31 ) + c32 X2 (t − τ32 ) − c33 X3 (t) dt, with initial data

Xi (θ) = ξi (θ), θ ∈ [−¯ τ , 0], τ¯ = max{τij },

(1)

where Xi (t) is the size of the prey i, i = 1, 2, and X3 (t) is the size of the predator. bi is the growth rate of species i, i = 1, 2, b3 is the death rate of the predator. cii is the intra-specific competition rate, i = 1, 2, 3. c12 and c21 are the inter-specific competition rates between species 1 and 2, c13 and c23 are the capture rates, c31 and c32 stand for the efficiency of food conversion. All coefficients mentioned above are positive constants. τij ≥ 0 represents ∗ Corresponding

author. Tel.: +86 0517 84183732. Email address: [email protected] (Meng Liu)

Preprint submitted to Elsevier

May 3, 2017

ACCEPTED MANUSCRIPT

the time delay. ξ(θ) = (ξ1 (θ), ξ2 (θ), ξ3 (θ))T ∈ U, U represents the space of all the continue functions from [−¯ τ , 0]

3 to R+ = {x = (x1 , x2 , x3 ) ∈ R3 |xi > 0, i = 1, 2, 3}.

On the other hand, the natural growth of species in the real world is inevitably subject to environmental fluctuations ([7]). It is therefore important to study the stochastic population models and reveal the effects of stochasticity on the dynamics of the models. There are several approaches to introduce the stochastic fluctuations

mainly affect the growth/death rates of the species, with ˙ i (t), i = 1, 2, 3, bi → bi + βi W

CR IP T

into population models. Here we follow the approach adopted in [8]-[30], i.e., we assume that the stochastic factors

where {W (t)}t≥0 = {W1 (t), W2 (t), W3 (t)}t≥0 is a three-dimensional Brownian motion defined on a complete

probability space (Ω, F, {F}t≥0 , P) with a filtration {F}t≥0 satisfying the usual conditions, and βi2 is the intensity

AN US

of the noise, i = 1, 2, 3. Thus we derive the following stochastic delay model:      dX (t) = X (t) b − c X (t) − c X (t − τ ) − c X (t − τ ) dt + β1 X1 (t)dW1 (t),  1 1 1 11 1 12 2 12 13 3 13       dX2 (t) = X2 (t) b2 − c21 X1 (t − τ21 ) − c22 X2 (t) − c23 X3 (t − τ23 ) dt + β2 X2 (t)dW2 (t),          dX3 (t) = X3 (t) − b3 + c31 X1 (t − τ31 ) + c32 X2 (t − τ32 ) − c33 X3 (t) dt + β3 X3 (t)dW3 (t), with initial value (1).

(2)

M

In the investigation of population model, an interesting question is to consider the stability of positive equilibrium state ([1]). However, lots of stochastic population models (for example, model (2)) do not have the

ED

traditional positive equilibrium state. Several researchers (see e.g., [8]-[16]) have suggested that one can consider the stability in distribution of the solutions (SDS) of stochastic population models. As far as we are concerned, no results relative to the SDS of model (2) have been reported. One reason is that the classical approaches could

PT

not be used.

One classical approach to investigate the SDS of stochastic population models is to solve the corresponding Fokker-Planck equations explicitly (see, e.g. [8]). However, this is almost impossible for stochastic population

CE

models with time delays. Another approach is the Markov semigroup theory, see, e.g., [9, 10, 27]. However, this approach requires some standard measure. The standard measure of the phase space of ordinary differential

AC

equations is the Lebesgue measure. But it is difficult to decide the standard measure of the phase space of delay differential equations (communicated with Professor Ryszard Rudnick). The third approach is the Lyapunov function method proposed in [31] which has been widely used in literatures, for example, by using this approach, Ji et al. [11, 12] studied the SDS of predator-prey models perturbed by white noises; Mao [13] explored the SDS of a n-species Lotka-Volterra mutualism system with random perturbations; Jiang et al. [14] considered the SDS of a n-species stochastic Lotka-Volterra competitive model; Zou et al. [15] investigated the SDS of a predator-prey model perturbed by white noises; Zhao et al. [16] analyzed the SDS of a three-species stochastic Lotka-Volterra competitive system in polluted environments. However, the Lyapunov function approach needs the Markov property of the solutions which the stochastic delay models do not own. The fourth approach is the basic theory of SDS for stochastic functional differential equations (see, e.g., Hu and Wang [32]). However, this approach requires that the coefficients of the models must obey the linear growth condition, while most population systems do not obey. 3

ACCEPTED MANUSCRIPT

In this paper, we use an asymptotic approach ([27]) to analyze the SDS of model (2). Under some assumptions we show that there are three numbers γ1 > γ2 > γ3 which are represented by the coefficients of the model. We prove that: if γ1 < 1, then lim Xi (t) = 0 almost surely (a.s.), i = 1, 2, 3; If γi > 1 > γi+1 , i = 1, 2, then t→+∞

the distribution of (X1 (t), ..., Xi (t))T converges weakly to a unique ergodic invariant distribution (UEID) and lim Xj (t) = 0 a.s., j = i + 1, ..., 3; If γ3 > 1, then the distribution of (X1 (t), X2 (t), X3 (t))T converges weakly

t→+∞

to a UEID a.s..

CR IP T

The organization of this paper is outlined as follows. In Section 2, we give the main theorem and its proof. In Section 3, we introduce some numerical simulations to illustrate the main results. In Section 4, we discuss the theoretical results and their biological interpretations.

2. Main results

AN US

For the sake of simplification, we define some notations:

β2 ¯3 βi2 3 , i = 1, 2, a3 = b3 + 3 , R + = {x = (x1 , x2 , x3 ) ∈ R |xi ≥ 0, i = 1, 2, 3}, 2 2 c11 c11 c12 c13 b1 β12 /2 Γ = c21 c22 c23 , b2 β22 /2 , C = c21 −c31 −c32 c33 −c31 −b3 β32 /2 b1 c11 c11 c12 c13 b1 c13 c12 b1 C1 = b2 c22 c23 , C2 = c21 b2 c23 , C3 = c21 c22 b2 , −b3 −c32 c33 −c31 −b3 c33 −c31 −c32 −b3 2 β1 /2 c12 c13 c11 β12 /2 c13 c11 c12 β12 /2 C˜1 = β22 /2 c22 c23 , C˜2 = c21 β22 /2 c23 , C˜3 = c21 c22 β22 /2 , 2 β3 /2 −c32 c33 −c31 β32 /2 c33 −c31 −c32 β32 /2

PT

ED

M

ai = bi −

Θ1 = c22 b1 − c12 b2 , Θ2 = c11 b2 − c21 b1 , Θ3 = c31 b1 − c11 b3 ,

CE

˜ 1 = c22 β12 /2 − c12 β22 /2, Θ ˜ 2 = c11 β22 /2 − c21 β12 /2, Θ ˜ 3 = c11 β32 /2 + c31 β12 /2. Θ

AC

For a function g, we denote the following notations: Z Z 1 t 1 t g(s)ds, hg(t)i∗ = lim sup g(s)ds, hg(t)i = t 0 t→+∞ t 0

hg(t)i∗ = lim inf t→+∞

1 t

Z

t

g(s)ds.

0

Lemma 1. For any initial value ξ(θ) = (ξ1 (θ), ξ2 (θ), ξ3 (θ))T ∈ U, there is a unique global positive solution X(t) to model (2) a.s.. Moreover, there is a positive constant K such that lim sup E(Xi (t)) ≤ K, i = 1, 2, 3.

(3)

t→∞

Proof. Consider the following system:     N1 (t) N2 (t−τ12 ) N3 (t−τ13 )  dN (t) = a − c e − c e − c e dt + β1 dW1 (t),  1 1 11 12 13       dN2 (t) = a2 − c21 eN1 (t−τ21 ) − c22 eN2 (t) − c23 eN3 (t−τ23 ) dt + β2 dW2 (t),          dN3 (t) = − a3 + c31 eN1 (t−τ31 ) + c32 eN2 (t−τ32 ) − c33 eN3 (t) dt + β3 dW3 (t), 4

(4)

ACCEPTED MANUSCRIPT

with initial value Ni (θ) = ln ξi (θ), i = 1, 2, 3. Since the coefficients of (4) obey the local Lipschitz condition, hence (4) has a unique local solution N (t) on [0, τe ), where τe stands for the explosion time. It then follows from Itˆo’s formula that (2) has a unique positive local solution X(t) = (X1 (t) = eN1 (t) , X2 (t) = eN2 (t) , X3 (t) = eN3 (t) )T .

CR IP T

Now we are in the position to show that X(t) is global, i.e., τe = +∞. Consider the following system:      dY (t) = Y (t) b − c Y (t) dt + β1 Y1 (t)dW1 (t),  1 1 1 11 1       dY2 (t) = Y2 (t) b2 − c22 Y2 (t) dt + β2 Y2 (t)dW2 (t),          dY3 (t) = Y3 (t) − b3 + c31 Y1 (t − τ31 ) + c32 Y2 (t − τ32 ) − c33 Y3 (t) dt + β3 Y3 (t)dW3 (t),

(5)

with initial value Yi (θ) = ξi (θ), i = 1, 2, 3. By the stochastic comparison theorem ([33]), one can see that for

AN US

t ∈ [0, τe ),

Xi (t) ≤ Yi (t) a.s., i = 1, 2, 3.

(6)

Thanks to Theorem 4.2 in [17], system (5) can be explicitly solved as follows

ED

M

 exp{a1 t + β1 W1 (t)}  , Y1 (t) = −1  Rt   Y1 (0) + c11 0 exp{a1 s + β1 W1 (s)}ds     exp{a2 t + β2 W2 (t)} , Y2 (t) = −1 Rt Y2 (0) + c22 0 exp{a2 s + β2 W2 (s)}ds   R  t   exp{−a3 t + 0 (c31 Y1 (s − τ31 ) + c32 Y2 (s − τ32 ))ds + β3 W3 (t)}   . Rt Rs  Y3 (t) = −1 Y3 (0) + c33 0 exp{−a3 s + 0 (c31 Y1 (u − τ31 ) + c32 Y2 (u − τ32 ))du + β3 W3 (s)}ds

Note that Y1 (t), Y2 (t) and Y3 (t) are existent on t ≥ 0, hence τe = +∞. The proof of (3) is standard and hence is omitted (see e.g. [18]).

PT

Before we state the main theorem of this paper, we need to introduce several hypotheses. First of all, without loss of generality, we assume that b1 /β12 > b2 /β22 , which implies that the persistence ability of species 1 is greater

CE

than that of species 2.

Hypothesis 1. C > 0, Ci > 0, i = 1, 2, 3, which imply that all the populations could coexist if model (2) frees

AC

from stochastic noises. C33 > 0, Θi > 0, i = 1, 2, which imply that the two prey populations could coexist if both environmental noises and the predator are absent, where Cij represents the complement minor of cij in the determinant C, i, j = 1, 2, 3. In addition, we also need the following two technical hypotheses.

Hypothesis 2. C˜i > 0, C23 > 0, C32 > 0, C31 < 0, Hypothesis 3. c11 > c12 + c13 , c22 > c21 + c23 , c33 > c31 + c32 . Theorem 1. For model (2), define ˜ 2 , γ3 = C3 /C˜3 , π2 = Θ3 /Θ ˜ 3 , π3 = C2 /C˜2 . γ1 = 2b1 /β12 , γ2 = Θ2 /Θ Under Hypotheses 1, 2 and 3, we have the following results: 5

ACCEPTED MANUSCRIPT

(I) if Γ > 0, then γ1 > γ2 > γ3 and, (i) if γ1 < 1, then lim Xi (t) = 0 a.s., i = 1, 2, 3; t→+∞

(ii) If γ1 > 1 > γ2 , then lim Xi (t) = 0 a.s., i = 2, 3, and there is a unique ergodic invariant distribution t→+∞

(UEID) µ1 such that the distributions of X1 (t) converge weakly to µ1 and Z Z a1 1 t X1 (s)ds = x1 µ1 (dx1 ) = lim , a.s.. t→+∞ t 0 c 11 R+ lim X3 (t) = 0 a.s., and there is a UEID µ2 such that the distributions of

t→+∞

T

(X1 (t), X2 (t)) converge weakly to µ2 and Z

1 t→+∞ t lim

t

Xi (s)ds =

0

Z

2 R+

xi µ2 (dx1 , dx2 ) =

CR IP T

(iii) If γ2 > 1 > γ3 , then

˜i Θi − Θ , a.s., i = 1, 2. C33

(iv) If γ3 > 1, then there is a UEID µ3 such that the distributions of (X1 (t), X2 (t), X3 (t))T converge weakly

1 t→+∞ t lim

Z

t

Xi (s)ds =

0

Z

3 R+

xi µ3 (dx1 , dx2 , dx3 ) =

(II) If Γ < 0, then γ1 > π2 > π3 and, (v) if γ1 < 1, then (i) holds; (vi) If γ1 > 1 > π2 , then (ii) holds;

(7)

lim X2 (t) = 0 a.s., and there is a UEID µ4 such that the distributions of

t→+∞

T

M

(vii) If π2 > 1 > π3 , then

Ci − C˜i , a.s., i = 1, 2, 3. C

AN US

to µ3 and

(X1 (t), X3 (t)) converge weakly to µ4 and Z

t

Xi (s)ds =

ED

1 t→+∞ t lim

0

Z

2 R+

xi µ4 (dx1 , dx3 ) =

˜i Θi − Θ , a.s., i = 1, 3. C22

(viii) If π3 > 1, then (iv) holds.

PT

Before we prove Theorem 1, let us recall some definitions. Definition 1. A population, x(t), is said to be extinctive (respectively, stable in the mean ([14])) if lim x(t) = 0 t→+∞

CE

(respectively, lim hx(t)i = a positive constant). t→+∞

ei (t)| = 0 a.s., i = 1, 2, 3, where Definition 2. System (2) is said to be globally attractive if lim E|Xi (t) − X t→+∞

AC

e e1 (t), X e2 (t), X e3 (t))T are two arbitrary solutions of system (2) with X(t) = (X1 (t), X2 (t), X3 (t))T and X(t) = (X e initial conditions X(θ) ∈ U and X(θ) ∈ U respectively.

3 ¯+ Definition 3. If there is a unique probability measure µ with nowhere zero density in R such that for arbitrary

X(θ) = ξ(θ) ∈ U, the transition probability p(t, ξ, ·) of X(t) converges weakly to µ with t → +∞, then system (2) is said to be stable in distribution. We divide the whole proof of Theorem 1 into several parts. 2.1. Extinction and stability in the mean Lemma 2. For arbitrarily τ ≥ 0, 1 t→+∞ t lim

Z

t

Yi (s)ds = 0, a.s., i = 1, 2.

t−τ

6

(8)

ACCEPTED MANUSCRIPT

Proof. Consider the first two equations in model (5), by [14], we have  Z  ai β2 1 t   lim Yi (s)ds = a.s., if bi ≥ i , i = 1, 2; t→+∞ t 0 cii 2 2  β  lim Y (t) = 0 a.s., if b < i , i = 1, 2.  i i t→+∞ 2

If bi ≥ βi2 /2, then 1 lim t→+∞ t

Z

t

1 Yi (s)ds = lim t→+∞ t t−τ

Z

t

0

Yi (s)ds −

Z

t−τ

0

CR IP T

Consequently, if bi < βi2 /2, then  Z t Z t−τ Z 1 1 t Yi (s)ds = 0. Yi (s)ds = lim Yi (s)ds − lim t→+∞ t t→+∞ t t−τ 0 0

 ai ai Yi (s)ds = − = 0. cii cii

Lemma 3. ([34]). Let h(t) be a continuous function from Ω × [0, +∞) to R+ . (I) If there exist three constants T > 0, λ0 > 0 and λ such that

AN US

ln h(t) ≤ λt − λ0

Z

t

h(s)ds +

0

n X

βi Wi (t)

i=1

for all t ≥ T , where Wi (t), 1 ≤ i ≤ n, are independent Brownian motions and βi , 1 ≤ i ≤ n, are constants, then    hh(t)i∗ ≤ λ/λ0 a.s., if λ ≥ 0, lim h(t) = 0

t→+∞

a.s.,

if

λ < 0.

M

 

ED

(II) If there exist three positive constants T, λ and λ0 such that ln h(t) ≥ λt − λ0

Z

t

h(s)ds +

0

n X

βi Wi (t)

i=1

PT

for all t ≥ T , where Wi (t), 1 ≤ i ≤ n, are independent Brownian motions and βi , 1 ≤ i ≤ n, are constants, then

 h(t) ∗ ≥ λ/λ0 a.s.

CE

Lemma 4. For model (5),

(a) If a1 < 0 and a2 < 0, then lim Yi (t) = 0, a.s., i = 1, 2, 3; t→+∞

AC

(b) If a1 ≥ 0, a2 < 0 and c31 a1 /c11 − a3 < 0, then Z 1 t a1 lim Y1 (s)ds = , t→+∞ t 0 c11

lim Y2 (t) = 0 , lim Y3 (t) = 0, a.s..

t→+∞

(c) If a1 ≥ 0, a2 < 0 and c31 a1 /c11 − a3 ≥ 0, then Z 1 t a1 lim Y1 (s)ds = , lim Y2 (t) = 0, t→+∞ t 0 c11 t→+∞

t→+∞

1 t→+∞ t lim

Z

7

1 t→+∞ t lim

Y3 (s)ds =

0

(d) If a1 < 0, a2 ≥ 0 and c32 a2 /c22 − a3 < 0, then Z 1 t a2 lim Y1 (t) = 0, lim Y2 (s)ds = , t→+∞ t→+∞ t 0 c22 (e) If a1 < 0, a2 ≥ 0 and c32 a2 /c22 − a3 ≥ 0, then Z 1 t a2 lim Y1 (t) = 0, lim Y2 (s)ds = , t→+∞ t→+∞ t 0 c22

t

c31 a1 − c11 a3 , a.s.. c11 c33

lim Y3 (t) = 0, a.s..

t→+∞

Z

0

t

Y3 (s)ds =

c32 a2 − c22 a3 , a.s.. c22 c33

ACCEPTED MANUSCRIPT

(f ) If a1 ≥ 0, a2 ≥ 0 and c31 a1 /c11 + c32 a2 /c22 − a3 < 0, then Z Z a1 a2 1 t 1 t Y1 (s)ds = Y2 (s)ds = lim , lim , t→+∞ t 0 c11 t→+∞ t 0 c22

lim Y3 (t) = 0, a.s..

t→+∞

(g) If a1 ≥ 0, a2 ≥ 0 and c31 a1 /c11 + c32 a2 /c22 − a3 ≥ 0, then Z Z 1 t a1 a2 1 t lim Y1 (s)ds = Y2 (s)ds = , lim , t→+∞ t 0 c11 t→+∞ t 0 c22 Z

t

Y3 (s)ds =

0

c22 c31 a1 + c11 c32 a2 − c11 c22 a3 , a.s.. c11 c22 c33

CR IP T

1 t→+∞ t lim

Proof. Applying Itˆ o’s formula to model (5) results in ln Y1 (t) − ln Y1 (0) = a1 t − c11

t

Y1 (s)ds + β1 W1 (t),

(9)

Y2 (s)ds + β2 W2 (t),

(10)

0

Z

t

0

AN US

ln Y2 (t) − ln Y2 (0) = a2 t − c22

Z

ln Y3 (t) − ln Y3 (0) Z t Z t Z t = − a3 t + c31 Y1 (s − τ31 )ds + c32 Y2 (s − τ32 )ds − c33 Y3 (s)ds + β3 W3 (t) 0

− c32

Z

t

Y1 (s)ds + c32

0

0

t

Y2 (s)ds − c33

0

t

t−τ32

Z

Y2 (s)ds −

Z

0

−τ32

Z

0

t

Y3 (s)ds − c31

0

 Y2 (s)ds + β3 W3 (t).

Z

t

t−τ31

Y1 (s)ds −

Z

0

−τ31

 Y1 (s)ds

(11)

M

= − a3 t + c31

Z

ED

Dividing both sides of (9), (10) and (11) by t, one can see that β1 W1 (t) 1 Y1 (t) ln = a1 − c11 hY1 (t)i + , t Y1 (0) t

CE

PT

β2 W2 (t) 1 Y2 (t) ln = a2 − c22 hY2 (t)i + , t Y2 (0) t Z t  Z 0 1 Y3 (t) c31 ln = − a3 + c31 hY1 (t)i + c32 hY2 (t)i − c33 hY3 (t)i − Y1 (s)ds − Y1 (s)ds t Y3 (0) t t−τ31 −τ31 Z t  Z 0 β3 W3 (t) c32 Y2 (s)ds − Y2 (s)ds + . − t t t−τ32 −τ32

(12) (13)

(14)

AC

And it is well known that

lim

t→+∞

Wi (t) = 0, a.s., i = 1, 2, 3. t

(15)

Now let us prove (a) firstly. Since a1 = b1 − β12 /2 < 0 and a2 = b2 − β22 /2 < 0, making use of Lemma 3, it is

not difficult to show that lim Y1 (t) = 0, a.s.

(16)

lim Y2 (t) = 0, a.s..

(17)

t→+∞

and t→+∞

When (16) and (17) are used in (11), one can easily see that for sufficiently large t, Z t ln Y3 (t) − ln Y3 (0) ≤ −a3 t + εt − c33 Y3 (s)ds + β3 W3 (t). 0

8

ACCEPTED MANUSCRIPT

It then follows from a3 = b3 + β32 /2 > 0 and (I) in Lemma 3 that lim Y3 (t) = 0 a.s.

t→+∞

and lim Y2 (t) = 0, a.s..

t→+∞

CR IP T

Now let us prove (b) and (c). Using Lemma 3, a1 ≥ 0 and a2 < 0, it is easy to verify that Z a1 1 t Y1 (s)ds = , a.s. lim t→+∞ t 0 c11

(18)

(19)

AN US

As a consequence, model (5) reduces to following model:       dY1 (t) = Y1 (t) b1 − c11 Y1 (t) dt + β1 Y1 (t)dW1 (t),      dY3 (t) = Y3 (t) − b3 + c31 Y1 (t − τ31 ) − c33 Y3 (t) dt + β3 Y3 (t)dW3 (t). Then with the similar proof of that in [35], it is not difficult to show that   lim Y3 (t) = 0, a.s., if c31 a1 /c11 − a3 < 0,  t→+∞ Z 1 t c31 a1 − c11 a3   lim Y3 (s)ds = , if c31 a1 /c11 − a3 ≥ 0. t→+∞ t 0 c11 c33

The proofs of (d) and (e) are similar to these of (b) and (c) respectively, and hence are omitted.

M

Now let us prove (f) and (g). Utilizing Lemma 3, a1 = b1 − β12 /2 ≥ 0 and a2 = b2 − β22 /2 ≥ 0, it is easy to 1 lim t→+∞ t

Z

1 lim t→+∞ t

Z

ED

derive that

PT

and

t

Y1 (s)ds =

a1 , a.s. c11

(20)

Y2 (s)ds =

a2 , a.s.. c22

(21)

0 t

0

lim t−1 ln Y1 (t) = 0, a.s.

(22)

lim t−1 ln Y2 (t) = 0, a.s..

(23)

t→+∞

t→+∞

AC

and

CE

Plugging (20), (21) into (12), (13), respectively, then combining (15) leads to

Besides, computing (12) × (−c22 c31 ) + (13) × (−c11 c32 ) + (14) × c11 c22 , one can derive that Y3 (t) c22 c31 Y1 (t) c11 c31 Y2 (t) c11 c22 ln − ln − ln t Y3 (0) t Y1 (0) t Y2 (0)

Z t  Z 0 c11 c22 c31 =c11 c22 a3 − c22 c31 a1 − c11 c32 a2 − c11 c22 c33 hY3 (t)i + Y1 (s)ds − Y1 (s)ds t t−τ31 −τ31 Z t  Z 0 c11 c22 c32 −c22 c31 β1 W1 (t) − c11 c32 β2 W2 (t) + c11 c22 β3 W3 (t) Y2 (s)ds − Y2 (s)ds + . + t t t−τ32 −τ32

In view of (8), (22), (23), (24) and Lemma 3, we could gain that   lim Y3 (t) = 0, a.s., if c31 a1 /c11 + c32 a2 /c22 − a3 < 0,  t→+∞ Z 1 t c22 c31 a1 + c11 c32 a2 − c11 c22 a3   lim Y3 (s)ds = , a.s., if c31 a1 /c11 + c32 a2 /c22 − a3 ≥ 0. t→+∞ t 0 c11 c22 c33 9

(24)

ACCEPTED MANUSCRIPT

Lemma 5. For arbitrarily τ ≥ 0,

1 t→+∞ t lim

Z

t

Y3 (s)ds = 0, a.s.

(25)

t−τ

1 t→+∞ t

Proof. From Lemma 4 we can see that either lim Y3 (t) = 0 or lim t→+∞

When lim Y3 (t) = 0, clearly,

Z

t

Y3 (s)ds = a constant.

0

t→+∞

1 t→+∞ t

When lim

Z

Z

t

1 Y3 (s)ds = lim t→+∞ t t−τ

Z

t

Y3 (s)ds −

0

Z

t−τ

0

t

Y3 (s)ds = a constant, it easy to see that

0

1 t→+∞ t lim

Z

t

1 t→+∞ t

Y3 (s)ds = lim

t−τ

Z

t

0

Y3 (s)ds −

Z

t−τ

0

Lemma 6. The solution X(t) of (2) satisfies

 Y3 (s)ds = 0, a.s..

AN US

This completes the proof.

 Y3 (s)ds = 0, a.s..

CR IP T

1 lim t→+∞ t

lim sup ln Xi (t)/t ≤ 0, a.s., i = 1, 2, 3. t→+∞

1 t→+∞ t

Proof. From Lemma 4, we know that either lim Yi (t) = 0, a.s. or lim If lim Yi (t) = 0, then t→+∞

t

Yi (s)ds = a constant, i = 1, 2, 3.

0

M

t→+∞

Z

(26)

lim sup ln Xi (t)/t ≤ lim sup ln Yi (t)/t ≤ 0, a.s.. t→+∞

Z

t

t→+∞

lim

t→+∞

t→+∞

ED

1 Yi (s)ds = a constant, similar to the proof of (22), we have t 0 lim sup ln Xi (t)/t ≤ lim sup ln Yi (t)/t = 0, a.s..

If

t→+∞

lim ln Yi (t)/t = 0 a.s., hence

t→+∞

PT

Theorem 2. For model (2), under Hypotheses 1 and 2, we have: (I’) If Γ > 0, then γ1 > γ2 > γ3 and,

(i’) If γ1 < 1, then lim Xi (t) = 0, a.s., i = 1, 2, 3;

CE

t→+∞

AC

(ii’) If γ1 > 1 > γ2 , then lim Xi (t) = 0, a.s., i = 2, 3, and t→+∞

1 t→+∞ t lim

Z

t

X1 (s)ds =

0

a1 , a.s.; c11

(iii’) If γ2 > 1 > γ3 , then lim X3 (t) = 0, a.s. and t→+∞

1 t→+∞ t lim

Z

t

X1 (s)ds =

0

(iv’) If γ3 > 1, then 1 t→+∞ t lim

Z

Z

˜1 Θ1 − Θ , C33

1 t→+∞ t

Xi (s)ds =

Ci − C˜i , a.s., i = 1, 2, 3. C

t

0

(II’) If Γ < 0, then γ1 > π2 > π3 and, (v’) If γ1 < 1, then (i’) holds; (vi’) If γ1 > 1 > π2 , then (ii’) holds;

10

lim

0

t

X2 (s)ds =

˜2 Θ2 − Θ , a.s.; C33

(27)

ACCEPTED MANUSCRIPT

(vii’) If π2 > 1 > π3 , then lim X2 (t) = 0, a.s. and t→+∞

1 t→+∞ t lim

Z

t

X1 (s)ds =

0

˜1 Θ1 − Θ , C22

1 t→+∞ t lim

Z

t

X3 (s)ds =

0

˜3 Θ3 − Θ , a.s.; C22

(viii’) If π3 > 1, then (iv’) holds.

b1 b2 > 2 , β12 /2 β2 /2 therefore,

2 2 ˜ 2 = c11 β2 − c21 β1 > 0, Θ 2 2

CR IP T

Proof. We only prove (I’), the proof of (II’) is similar. Note that Θ2 = c11 b2 − c21 b1 > 0,

Θ2 c11 2b1 − = 2 (b1 β22 − b2 β12 ) > 0. ˜ ˜2 β12 Θ2 β1 Θ

On the other hand, we can compute that Θ2 Θ3 c11 Γ − = > 0, ˜2 ˜3 ˜ 2Θ ˜3 Θ Θ Θ

AN US

C3 c23 Γ Θ3 − = > 0. ˜3 ˜ 3 C˜3 Θ C˜3 Θ

Hence,

Θ2 C3 b1 > > . ˜2 β12 /2 Θ C˜3

M

Applying Itˆo’s formula to (2) results in Z t Z t Z t ln X1 (t) − ln X1 (0) = a1 t − c11 X1 (s)ds − c12 X2 (s − τ12 )ds − c13 X3 (s − τ13 )ds + β1 W1 (t), 0

t

0

X1 (s − τ21 )ds − c22

ED

ln X2 (t) − ln X2 (0) = a2 t − c21

Z

0

ln X3 (t) − ln X3 (0) = −a3 t + c31

Z

0

t

(28)

0

Z

X1 (s − τ31 )ds + c32

t

0

Z

0

X2 (s)ds − c23 t

Z

0

t

X3 (s − τ23 )ds + β2 W2 (t),

X2 (s − τ32 )ds − c33

Z

(29)

t

X3 (s)ds + β3 W3 (t).

(30)

0

PT

Dividing both sides of (28), (29) and (30) by t, one can obtain that

CE

Z t  Z 0 c12 1 X1 (t) ln =a1 − c11 hX1 (t)i − c12 hX2 (t)i − c13 hX3 (t)i + X2 (s)ds − X2 (s)ds t X1 (0) t t−τ12 −τ12 Z t  Z 0 β1 W1 (t) c13 X3 (s)ds − X3 (s)ds + , + t t t−τ13 −τ13

AC

Z t  Z 0 1 X2 (t) c21 ln =a2 − c21 hX1 (t)i − c22 hX2 (t)i − c23 hX3 (t)i + X1 (s)ds − X1 (s)ds t X2 (0) t t−τ21 −τ21 Z t  Z 0 c23 β2 W2 (t) + X3 (s)ds − X3 (s)ds + , t t t−τ23 −τ23 Z t  Z 0 1 X3 (t) c31 ln = − a3 + c31 hX1 (t)i + c32 hX2 (t)i − c33 hX3 (t)i − X1 (s)ds − X1 (s)ds t X3 (0) t t−τ31 −τ31  Z t Z 0 c32 β3 W3 (t) X2 (s)ds + − X2 (s)ds − . t t −τ32 t−τ32

(31)

(32)

(33)

To begin with, let us prove (i’). Notice that a1 < 0, a2 < 0, then the conditions in (a) of Lemma 4 are satisfied. Therefore, lim Xi (t) = 0, a.s., i = 1, 2, 3.

t→+∞

That is to say, all the populations go to extinction a.s.. 11

ACCEPTED MANUSCRIPT

Now we are in the position to show (ii’). Multiplying both sides of (31) and (32) by (−c21 ) and c11 , respectively, and then adding these two equations, we can observe that c11 X2 (t) c21 X1 (t) ln − ln t X2 (0) t X1 (0)

CR IP T

Z t  Z 0 ˜ 2 − C33 hX2 (t)i − C32 hX3 (t)i − c12 c21 X2 (s)ds − X2 (s)ds =Θ2 − Θ t t−τ12 −τ12   Z t Z t Z 0 Z 0 c11 c21 c13 c21 X3 (s)ds − X3 (s)ds + X2 (s)ds − X2 (s)ds − t t t−τ13 −τ13 t−τ21 −τ21 Z t  Z 0 −c21 β1 W1 (t) + c11 β2 W2 (t) + X3 (s)ds − X3 (s)ds + . t t−τ23 −τ23

Plugging (26), (35) into (34) leads to

(35)

AN US

In view of (6), (8) and (25), we can obtain that Z 1 t lim Xi (s)ds = 0, a.s., i = 1, 2, 3. t→+∞ t t−τ

(34)

c11 X2 (t) ˜ 2 + ε − C33 hX2 (t)i + −c21 β1 W1 (t) + c11 β2 W2 (t) ln ≤ Θ2 − Θ t X2 (0) t ˜ 2 +ε < 0, then it can be verified straightforwardly for sufficient large t, where ε is sufficiently small fulfilling Θ2 − Θ that

lim X2 (t) = 0, a.s..

M

t→+∞

ED

As a consequence, system (2) reduces to following two-species model:      dX1 (t) = X1 (t) b1 − c11 X1 (t) − c13 X3 (t − τ13 ) dt + β1 X1 (t)dW1 (t),       dX3 (t) = X3 (t) − b3 + c31 X1 (t − τ31 ) − c33 X3 (t) dt + β3 X3 (t)dW3 (t).

PT

Then with the similar proof of that in [35], it is not difficult to show that Z 1 t a1 lim X1 (s)ds = , lim X3 (t) = 0, a.s.. t→+∞ t 0 t→+∞ c11

AC

CE

Now let us move to (iii’). Let m and n satisfy the following system:   c22 m − c32 n = c12 , 

Thus,

m=

c23 m + c33 n = c13 .

C21 > 0, C11

n=−

After a calculation of (31) × (−1) + (32) × m + (33) × n, we get

C31 > 0. C11

Z t  Z 0 1 X1 (t) m X2 (t) n X3 (t) C1 − C˜1 C c12 ln = ln + ln + − hX1 (t)i + X2 (s)ds − X2 (s)ds t X1 (0) t X2 (0) t X3 (0) C11 C11 t t−τ12 −τ12 Z t  Z t  Z 0 Z 0 c13 c21 m X3 (s)ds − X3 (s)ds − X1 (s)ds − X1 (s)ds + t t t−τ13 −τ13 t−τ21 −τ21 (36) Z t  Z t  Z 0 Z 0 c23 m c31 n − X3 (s)ds − X3 (s)ds + X1 (s)ds − X1 (s)ds t t t−τ23 −τ23 t−τ31 −τ31 Z t  Z 0 c32 n β1 W1 (t) − mβ2 W2 (t) − nβ3 W3 (t) + X2 (s)ds − X2 (s)ds + . t t t−τ32 −τ32 12

ACCEPTED MANUSCRIPT

Substituting (26) and (35) into (36) gives 1 X1 (t) C1 − C˜1 C β1 W1 (t) − mβ2 W2 (t) − nβ3 W3 (t) ln ≤ +ε− hX1 (t)i + t X1 (0) C11 C11 t

(37)

for sufficiently large t. By C1 /C˜1 > 2b1 /β12 > γ2 > 1, Lemma 3 and the arbitrariness of ε, we have C1 − C˜1 , a.s.. C

Analogously, let m ˜ and n ˜ satisfy the following system:   c11 m ˜ − c31 n ˜ = c21 , 

Then, m ˜ =

c13 m ˜ + c33 n ˜ = c23 .

C21 > 0, C22

n ˜=−

C32 > 0. C22

(38)

CR IP T

hX1 (t)i∗ ≤

we get

AN US

Multiplying both sides of (31), (32) and (33) by m, ˜ (-1) and n ˜ , respectively, and then adding these three equations, hX2 (t)i∗ ≤ Plugging (35), (38) and (39) into (33) results in

C2 − C˜2 , a.s.. C

(40)

M

1 X3 (t) c33 (C3 − C˜3 ) β3 W3 (t) ln ≤ + ε − c33 hX3 (t)i + t X3 (0) C t

(39)

for sufficiently large t. Then according to C3 < C˜3 , Lemma 3 and the arbitrariness of ε, we can see that lim X3 (t) = 0, a.s..

ED

t→+∞

PT

Therefore, system (2) becomes the following two-species model:       dX1 (t) = X1 (t) b1 − c11 X1 (t) − c13 X3 (t − τ13 ) dt + β1 X1 (t)dW1 (t),      dX2 (t) = X2 (t) b2 − c21 X1 (t − τ21 ) − c22 X2 (t) dt + β2 X2 (t)dW2 (t).

CE

Then with the similar proof of that in [35], it follows that Z Z ˜1 ˜2 1 t Θ1 − Θ 1 t Θ2 − Θ lim X1 (s)ds = , lim X2 (s)ds = , a.s.. t→+∞ t 0 t→+∞ t 0 C33 C33

AC

(iv’). Note that C3 − C˜3 > 0. Then an application of Lemma 3 to (40) yields that hX3 (t)i∗ ≤

C3 − C˜3 , a.s.. C

Substituting (35), (39) and (41) into (31) elicits that 1 X1 (t) c12 (C2 − C˜2 ) c13 (C3 − C˜3 ) ln ≥a1 − ε − c11 hX1 (t)i − − t X1 (0) C C Z t  Z 0 c12 + X2 (s)ds − X2 (s)ds t t−τ12 −τ12 Z t  Z 0 c13 β1 W1 (t) + X3 (s)ds − X3 (s)ds + t t t−τ13 −τ13 ≥

c11 (C1 − C˜1 ) β1 W1 (t) − 2ε − c11 hX1 (t)i + C t 13

(41)

ACCEPTED MANUSCRIPT

for sufficiently large t. Hence, we can further get from Lemma 3 and the arbitrariness of ε that hX1 (t)i∗ ≥

C1 − C˜1 , a.s.. C

(42)

hX2 (t)i∗ ≥

C2 − C˜2 , a.s.. C

(43)

In the same way, we can show that

hX3 (t)i∗ ≥

CR IP T

Substituting (42) and (43) into (33), and then using (35), Lemma 3 and the arbitrariness of ε, we get C3 − C˜3 , a.s.. C

This, together with (38), (39), (41), (42) and (43), yields 1 lim t→+∞ t

Z

t

Xi (s)ds =

0

Ci − C˜1 , a.s. i = 1, 2, 3. C

AN US

This completes the proof. 2.2. Global attractivity

Theorem 3. Let Hypothesis 3 hold, then system (2) is globally attractive.

ED

M

Proof. Let hi be the cofactor of the i-th diagonal element of LC , where  c + c13 −c12 −c13  12  LC =  −c21 c21 + c23 −c23  −c31 −c32 c31 + c32



  . 

Making use of Kirchhoff’s Matrix Tree Theorem (see, e.g., [36]), one has hi > 0. Define  Z 3 X e cij hi ln Xi (t) − ln Xi (t) +

PT

V (t) =

3 X i=1

j=1,j6=i

t

t−τij

Reckoning the right differential dV (t) deduces that 3 X

CE dV (t) =

AC

i=1

+

ei (t))d(Xi (t) − X ei (t)) hi sgn(Xi (t) − X

3 3 X X

i=1 j=1,j6=i

≤

3  X i=1

According to Theorem 2.3 in [37], 3 3 X X

i=1 j=1,j6=i

As a consequence, E(V (t)) +

 ej (s) ds . Xj (s) − X

  ej (t) − Xj (t − τij ) − X ej (t − τij ) hi cij Xj (t) − X

 3 X ei (t) + e dt. − hi cii Xi (t) − X h c X (t) − X (t) i ij j j j=1,j6=i

3 3 X X e e hi cij Xj (t) − Xj (t) = hi cij Xi (t) − Xi (t) .

Z tX 3 0 i=1

i=1 j=1,j6=i

  3 X ei (s) ds ≤ V (0) < ∞. hi cii − |cij | E Xi (s) − X j=1,j6=i

14

ACCEPTED MANUSCRIPT

Therefore

ei (t) ∈ L1 [0, +∞). E Xi (t) − X

According to (2),

(44)

 Z t 2 E(X1 (t)) = X1 (0) + E(X1 (s))b1 − c11 E(X1 (s)) − c12 E(X1 (s)X2 (s − τ12 )) − c13 E(X1 (s)X3 (s − τ13 )) ds. 0

CR IP T

Clearly, E(X1 (t)) is differentiable. By virtue of (3), dE(X1 (t)) = E(X1 (t))b1 − c11 E(X1 (t))2 − c12 E(X1 (t)X2 (t − τ12 )) − c13 E(X1 (t)X3 (t − τ13 )) dt ≤ E(X1 (t))b1 ≤ b1 K1 ,

where K1 > 0 is a constant. Consequently, E(X1 (t)) is uniformly continuous. In the same way, E(X2 (t)) and E(X3 (t)) are also uniformly continuous functions. By virtue of (44) and Barbalat’s conclusion ([38]), we derive

2.3. Stability in distribution

AN US

the required assertion.

Theorem 4. If Hypotheses 1-3 hold, Γ > 0 (or Γ < 0) and γ3 > 1, then system (2) is stable in distribution. Proof. We only present the proof of Γ > 0, the proof of Γ < 0 is analogous. Let P (t, ξ, M ) stand for the probability of event X(t; ξ) ∈ M with initial value ξ(θ) ∈ U. In view of Lemma 3 and Chebyshev’s inequality,

M

one can observe that the family of p(t, ξ, ·) is tight.

Denote by P(U) all the probability measures on U. For any P1 , P2 ∈ P, define the following Kantorovich

ED

metric (see, e.g. [39, 40]):

Z dF (P1 , P2 ) = sup f ∈F

3 R+

PT

where

f (X)P1 (dX) −

F =



Z

3 R+

f (X)P2 (dX) ,

 f : U → R |f (φ) − f (ϕ)| ≤ ||φ − ϕ||, |f (·)| ≤ 1 .

AC

CE

For any f ∈ F and t, s > 0, one obtains Ef (X(t + s; ξ)) − Ef (X(t; ξ))    = E E f (X(t + s; ξ))|Fs − Ef (X(t; ξ)) Z = Ef (X(t; ζ))p(s, ξ, dζ) − Ef (X(t; ξ)) ≤

Z

3 R+

3 R+

Ef (X(t; ζ)) − Ef (X(t; ξ)) p(s, ξ, dζ).

According to Lemma 3, there is a T > 0 such that for t ≥ T , sup Ef (X(t; ζ)) − Ef (X(t; ξ)) ≤ ε. f ∈F

In other words, |Ef (X(t + s; ξ)) − Ef (X(t; ξ))| ≤ 3ε. Thanks to the arbitrariness of f , one can derive that sup Ef (X(t + s; ξ)) − Ef (X(t; ξ)) ≤ 3ε. f ∈F

15

ACCEPTED MANUSCRIPT

Hence for ∀ t ≥ T, s > 0,

  dF p(t + s, ξ, ·), p(t, ξ, ·) ≤ 3ε,

Consequently, for arbitrary ξ ∈ U, {p(t,ξ, ·) : t ≥ 0}  is Cauchy in P(U). Let ξ(θ) ≡ 0.5. Thereby there is a unique µ(·) ∈ P(U) satisfying lim dF p(t, ξ, ·), µ(·) = 0. Since Γ > 0 and γ3 > 1, then µ is with nowhere t→+∞

zero density. By Lemma 3,

  lim dF p(t, ϕ, ·), p(t, ξ, ·) = 0.

CR IP T

t→+∞

Hence for any ϕ ∈ U       lim dF p(t, ϕ, ·), µ(·) ≤ lim dF p(t, ϕ, ·), p(t, ξ, ·) + lim dF p(t, ξ, ·), µ(·) = 0. t→+∞

t→+∞

t→+∞

which is the desired assertion.

Proof of Theorem 1. First of all, let us prove (iv). In view of (iv’) in Theorem 2, we have (27).

AN US

Theorem 4 means that system (2) is stable in distribution. That is to say, there is a unique probability 3 , such that p(t, ξ, ·) converges weakly to µ3 . It then measure, denoted by µ3 , with nowhere zero density in R+

follows from Kolmogorov-Chapman equation that µ3 is invariant. Thanks to Corollary 3.4.3 in [41], one gets µ3 is strong mixing and therefore is ergodic ([41]). Thanks to (3.3.2) in [41] and (27), we obtain

t→+∞

Z

3 R+

xi µ3 (dx1 , dx2 , dx3 ) =

Ci − C˜i , i = 1, 2, 3, a.s.. C

M

lim hXi (t)i =

Now let us show (iii). In view of (iii’) in Theorem 2,

t→+∞

lim hX1 (t)i =

ED

lim X3 (t) = 0,

t→+∞

˜1 Θ1 − Θ , C33

lim hX2 (t)i =

t→+∞

˜2 Θ2 − Θ , a.s.. C33

CE

PT

Model (2) reduces to the following competitive model     ˆ 1 (t) = X ˆ 1 (t) b1 − c11 X ˆ 1 (t) − c11 X ˆ 2 (t − τ12 ) dt + β1 X ˆ 1 (t)dW1 (t),  dX      ˆ 2 (t) = X ˆ 2 (t) b2 − c21 X ˆ 1 (t − τ21 ) − c22 X ˆ 2 (t) dt + β2 X ˆ 2 (t)dW2 (t),  dX

(45)

ˆ i (θ) = Xi (θ), i = 1, 2. Similar to the proof of (iv), there is a UEID denoted by µ2 such with initial value X ˆ 1 (t), X ˆ 2 (t))T converges weakly to µ2 . Note that lim X3 (t) = 0 a.s., then that the transition probability of (X t→+∞

AC

ˆ 1 (t), X ˆ 2 (t))T . This completes the proof of (iii). The (X1 (t), X2 (t))T has the same asymptotic properties with (X proof of (ii) is analogous, and (i’) in Theorem 2 means (i).

3. Numerical Simulations Now let us introduce some numerical figures to illustrate the theoretical results. Choosing the parameters as follows:      dX1 (t) = X1 (t) 0.6 − 0.5X1 (t) − 0.2X2 (t − 1) − 0.06X3 (t − 2) dt + β1 X1 (t)dW1 (t),        dX2 (t) = X2 (t) 0.4 − 0.1X1 (t − 2) − 0.4X2 (t) − 0.04X3 (t − 1) dt + β2 X2 (t)dW2 (t),          dX3 (t) = X3 (t) − 0.01 + 0.2X1 (t − 3) + 0.05X2 (t − 2) − 0.4X3 (t) dt + β3 X3 (t)dW3 (t), 16

(46)

ACCEPTED MANUSCRIPT

with initial value X1 (θ) = 0.5 + 0.1 sin θ, X2 (θ) = 0.25 + 0.07 sin θ, X3 (θ) = 0.15 + 0.05 sin θ. It is easy to check that C = 0.0759, C1 = 0.0642, C2 = 0.0561, C3 = 0.0372, C23 = 0.019, C31 = −0.016, C32 = 0.014, C33 = 0.18, Θ1 = 0.16, Θ2 = 0.14, Θ3 = 0.115. (i) Set β12 /2 = 0.3, β22 /2 = 0.25 and β32 /2 = 0.05. Then by direct calculation, we have b1 /β12 = 1 > b2 /β22 = 0.8

CR IP T

and C˜1 = 0.0271 > 0, C˜2 = 0.0379 > 0, C˜3 = 0.0278 > 0, Γ = 0.0020 > 0. Furthermore, γ3 = C3 /C˜3 = 1.3405 > 1. Therefore all conditions in (iv) of Theorem 2 have been checked. Then X(t) is SDS and Z Z C1 − C˜1 C2 − C˜2 1 t 1 t X1 (s)ds = X2 (s)ds = lim = 0.4895, lim = 0.2398, t→+∞ t 0 t→+∞ t 0 C C Z 1 t C3 − C˜3 lim X3 (s)ds = = 0.1245. t→+∞ t 0 C

the numerical method mentioned in [43, 44].

AN US

See Fig.1. Fig.1(a) is obtained by applying the Milstein method ([42]) and Fig.1(b) is obtained by using (ii) Choose β12 /2 = 0.3, β22 /2 = 0.25 and β32 /2 = 0.2. Then we get b1 /β12 = 1 > b2 /β22 = 0.8 and C˜1 = 0.0247 > ˜ 2 = 1.4737 > 1 > γ3 = 0, C˜2 = 0.0358 > 0, C˜3 = 0.0548 > 0, Γ = 0.0230 > 0. Note that γ2 = Θ2 /Θ C3 /C˜3 = 0.6795. Consequently, all requirements in (iii) of Theorem 2 hold. According to Theorem 2, we

M

can easily derive that the species 3 goes to extinction and Z Z ˜1 ˜2 1 t Θ1 − Θ 1 t Θ2 − Θ lim X1 (s)ds = = 0.5, lim X2 (s)ds = = 0.25. t→+∞ t 0 t→+∞ t 0 C33 C33 See Fig.2.

ED

(iii) Set β12 /2 = 0.3, β22 /2 = 0.44 and β32 /2 = 0.2. Then we can see that b1 /β12 = 1 > b2 /β22 = 0.4545 and C˜1 = 0.0089 > 0, C˜2 = 0.0761 > 0, C˜3 = 0.0520 > 0, Γ = 0.0011 > 0. Then we can obtain that

PT

˜ 2 = 0.7368, which implies that all requirements in (ii) of Theorem 2 are γ1 = 2b1 /β12 = 2 > 1 > γ2 = Θ2 /Θ

CE

met. Making use of Theorem 2 yields that X2 and X3 go to extinction and Z a1 1 t X1 (s)ds = = 0.6. lim t→+∞ t 0 c11 See Fig.3.

(iv) Choose β12 /2 = 0.65, β22 /2 = 0.44 and β32 /2 = 0.2. Then we can observe that b1 /β12 = 0.4615 > b2 /β22 =

AC

0.4545 and C˜1 = 0.0656 > 0, C˜2 = 0.593 > 0, C˜3 = 0.0782 > 0, Γ = 0.0288 > 0. What is more, γ1 = 2b1 /β12 = 0.9231 < 1. To sum up, all conditions in (i) of Theorem 2 are satisfied. As a result, all the

species in model (2) go to extinction and the simulation results are shown in Fig.4.

4. Conclusions and Discussions This paper is concerned with the SDS of a stochastic delay population model with two competing preys and one predator. The classical SDS problem of this type of model is difficult because the traditional approaches can not be used. In this paper we have utilized an asymptotic approach to investigate this question, and established the sharp sufficient criteria for the SDS of each species. Our main result is Theorem 2, which reveals the effects of stochastic perturbations on the persistence and extinction of X1 , X2 , X3 . Here we only discuss the case Γ > 0. 17

ACCEPTED MANUSCRIPT

(a) For the prey species 1, Theorem 2 demonstrates that the extinction or not of this species only relies on γ1 . It is clear to know that d(γ1 )/d(β12 ) < 0, d(γ1 )/d(β22 ) = 0 and d(γ1 )/d(β32 ) = 0. Hence, with the increasing of noise intensity β12 , the prey species 1 tends to extinction, but β22 and β32 have no impact on the extinction or not of X1 . (b) For the prey species 2, Theorem 2 presents that the extinction or not of this species only relies on γ2 . By simple calculation, we can show that d(γ2 )/d(β12 ) > 0, d(γ2 )/d(β22 ) < 0 and d(γ2 )/d(β32 ) = 0. As a no impact.

CR IP T

result, the species 2 tends to extinction with the increasing of β22 or the decreasing of β12 , whereas β32 has

(c) For the predator species 3, we can see that the extinction or not of this species only relies on γ3 . Thanks to d(γ3 )/d(βi2 ) < 0, i = 1, 2, 3,

we can obtain that the species 3 tends to extinction with the increasing of βi2 , i = 1, 2, 3.

AN US

(d) The position of asymptotic invariant distribution indicates the persistence level of a population. The right the position, the better the persistence level. Our Theorem 2 shows that the stochastic perturbations can affect the persistent level of all species. Here we only consider the case γ3 > 0. It is easy to calculate that d(C˜1 )/d(β12 ) > 0, d(C˜1 )/d(β22 ) < 0 and d(C˜1 )/d(β32 ) < 0. Hence, the position of asymptotic invariant distribution of species 1 will tend to left with the increasing of β12 , and will tend to right with the increasing

M

of β22 and β32 . That is to say, with the increasing of β12 , β22 and β32 , the species 1 will be persistent worse, better and better, respectively. In the same way, we can analyze the species 2 and 3, and the conclusions

ED

are shown in the following Table 1.

β12 ↑

β22 ↑

β32 ↑

Species 1

←

→

→

Species 2

←

←

→

Species 3

←

←

←

CE

PT

Table 1: Effects of stochasticity on the position of asymptotic distribution when Γ > 0 and γ3 > 1

Some questions deserve further exploration. In the first place, it is attractive to study the delay population

AC

model with other perturbations, such as Markovian switching (see, e.g., [19]), L´evy jumps (see, e.g., [22]), or reaction-diffusion (see e.g., [45]). Another problem is to consider population models with different functional responses, such as Holling II-IV type schemes and Beddington-DeAngelis functional response. We leave these investigations for future work.

Acknowledgements The authors thank the editor and the anonymous referees for those valuable comments. This research was supported by the Natural Science Foundation of PR China [grant number 11301207], Project Funded by Chinese Postdoctoral Science Foundation [grant numbers 2015M571349, 2016T90236], “333 High-level Personnel Training Project” and Qing Lan Project of Jiangsu province.

18

ACCEPTED MANUSCRIPT

References [1] Y. Kuang, Delay Differential Equations with Applications in Population Dynamics, Academic Press, Boston, 1993. [2] H. Freedman, P. Waltman, Mathematical analysis of some three-species food-chain models, Math. Biosci. 33 (1977) 257-276.

Bull. Math. Biol. 45 (1983) 877-900.

CR IP T

[3] Y. Takeuchi, N. Adachi, Existence of bifurcation of stable equilibrium in two-prey, one-predator communities,

[4] V. Hutson, G. Vickers, A criterion for permanent co-existence of species, with an application to a two-prey one-predator system, Math. Biosci. 63 (1983) 253-269.

[5] W. Feng, Coexistence, stability, and limiting behavior in a one-predator-two-prey model, J. Math. Anal. Appl.

AN US

179 (1993) 592-609.

[6] S. Ahmad, I.M. Stamova, Almost necessary and sufficient conditions for survival of species, Nonlinear Anal. Real World Appl. 5 (2004) 219-229.

[7] R. M. May, Stability and Complexity in Model Ecosystems, Princeton Univ., 1973.

M

[8] S. Pasquali, The stochastic logistic equation: stationary solutions and their stability, Rend. Sem. Mat. Univ. Padova 106 (2001) 165-183.

ED

[9] R. Rudnicki, Long-time behaviour of a stochastic prey-predator model, Stoch. Process. Appl. 108 (2003) 93-107.

(2007) 108-119.

PT

[10] R. Rudnicki, K. Pich´ or, Influence of stochastic perturbation on prey-predator systems, Math. Biosci. 206

CE

[11] 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. [12] C. Ji, D. Jiang, Dynamics of a stochastic density dependent predator-prey system with Beddington-DeAngelis

AC

functional response, J. Math. Anal. Appl. 381 (2011) 441-453.

[13] X. Mao, Stationary distribution of stochastic population systems, Systems Control Lett. 60 (2011) 398-405. [14] D. Jiang, C. Ji, X. Li, D. O’Regand, Analysis of autonomous Lotka-Volterra competition systems with random perturbation, J. Math. Anal. Appl. 390 (2012) 582-595. [15] X. Zou, D. Fan, K. Wang, Stationary distribution and stochastic Hopf bifurcation for a predator-prey system with noises, Discrete Contin. Dyn. Syst. Ser. B 18 (2013) 1507-1519. [16] Y. Zhao, S. Yuan, J. Ma, Survival and stationary distribution analysis of a stochastic competitive model of three species in a polluted environment, Bull. Math. Biol. 77 (2015) 1285-1326.

19

ACCEPTED MANUSCRIPT

[17] D. Jiang, N. Shi, A note on nonautonomous logistic equation with random perturbation, J. Math. Anal. Appl. 303 (2005) 164-172. [18] L. Hung, Stochastic delay population systems, Appl. Anal. 88 (2009) 1303-1320. [19] Q. Luo, X. Mao, Stochastic population dynamics under regime switching II, J. Math. Anal. Appl. 355 (2009) 577-593.

CR IP T

[20] C. Zhu, G. Yin, On hybrid competitive Lotka-Volterra ecosystems, Nonlinear Anal. 71 (2009) e1370-e1379. [21] X. Li, X. Mao, Population dynamical behavior of non-autonomous Lotks-Volterra competitive system with random perturbation, Discrete Contin. Dyn. Syst. 24 (2009) 523-545.

[22] J. Bao, X. Mao, G. Yin, C. Yuan, Competitive Lotka-Volterra population dynamics with jumps, Nonlinear Anal. 74 (2011) 6601-6616.

AN US

[23] W. Li, X. Zhang, C. Zhang, A new method for exponential stability of coupled reaction diffusion systems with mixed delays: Combining Razumikhin method with graph theory, J. Franklin Inst. 352 (2015) 1169-1191. [24] C. Zhang, W. Li, K. Wang, Graph-theoretic method on exponential synchronization of stochastic coupled networks with Markovian switching, Nonlinear Anal. Hybrid Syst. 15 (2015) 37-51. [25] C. Zhang, W. Li, K. Wang, Graph-theoretic approach to stability of multi-group models with dispersal,

M

Discrete Contin. Dyn. Syst. Ser. B 20 (2015) 259-280.

[26] M.Liu, C.Bai, Dynamics of a stochastic one-prey two-predator model with L´evy jumps, Appl. Math. Comput.

ED

284 (2016) 308-321.

[27] M.Liu, M.Fan, Stability in distribution of a three-species stochastic cascade predator-prey system with time

PT

delays, IMA J. Appl. Math. 82 (2017) 396-423. [28] M.Liu, C.Bai, Optimal harvesting of a stochastic delay competitive model, Discrete Contin. Dyn. Syst. Ser.

CE

B 22 (2017) 1493-1508.

[29] Y.Zhu, M.Liu, Permanence and extinction in a stochastic service-resource mutualism model, Appl. Math.

AC

Lett. 69 (2017) 1-7.

[30] M. Liu, C. Bai, Y. Jin, Population dynamical behavior of a two-predator one-prey stochastic model with time delay, Discrete Contin. Dyn. Syst. 37 (2017) 2513-2538.

[31] R.Z. Hasminskii, Stochastic Stability of Differential Equations, Netherlands, 1980. [32] G. Hu, K. Wang, Stability in distribution of neutral stochastic functional differential equations with Markovian switching, J. Math. Anal. Appl. 385 (2012) 757-769. [33] N. Ikeda, S. Watanabe, A comparison theorem for solutions of stochastic differential equations and its applications, Osaka J. Math. 14 (1977) 619-633. [34] M.Liu, C.Bai, Analysis of a stochastic tri-trophic food-chain model with harvesting, J. Math. Biol. 73 (2016) 597-625. 20

ACCEPTED MANUSCRIPT

[35] Q. Liu, Q. Chen, Z. Liu, Analysis on stochastic delay Lotka-Volterra systems driven by L´evy noise, Appl. Math. Comput. 235 (2014) 261-271. [36] J. Moon, Counting labelled tress, Canadian Mathematical Congress, Montreal, 1970. [37] M.Y. Li, Z. Shuai, Global-stability problem for coupled systems of differential equations on networks, J. Differential Equations 248 (2010) 1-20.

CR IP T

[38] I. Barbˇalat, Syst`ems d’´equations diff´erentielles d’oscillations non lin´eaires, Rev. Math. Pures Appl. 4 (1959) 267-270.

[39] C.A. Cabrelli, U.M. Molter, The Kantorovich metric for probability measures on the circle, J. Comput. Appl. Math. 57 (1995) 345-361.

[40] Y. Deng, W. Du, The Kantorovich metric in computer science: a brief survey, Electron. Notes Theor.

AN US

Comput. Sci. 253 (2009) 73-82.

[41] D. Prato, J. Zabczyk, Ergodicity for Infinite Dimensional Systems, Cambridge University Press, Cambridge, 1996.

[42] D.J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM

M

Rev. 43(3) (2001) 525-546.

[43] Y. Cai, Y. Kang, M. Banerjee, W. Wang, A stochastic epidemic model incorporating media coverage.

ED

Commun. Math. Sci 14 (2016) 892-910.

[44] Y. Cai, Y. Kang, W. Wang, A stochastic SIRS epidemic model with nonlinear incidence rate. Appl. Math.

PT

Comput. 305 (2017) 221-240.

[45] C.Bai, Multiplicity of solutions for a class of nonlocal elliptic pperators systems, Bull. Korean Math. Soc.,

AC

CE

doi: 10.4134/BKMS.b150489.

21

ACCEPTED MANUSCRIPT

0.7

3 X1(t)

0.6

X2(t)

2.5

X3(t)

0.5

t−1∫t0X1(s)ds

Density of X3

2

t−1∫t0X2(s)ds

0.4

t−1∫t X (s)ds

1.5

0 3

Density of X2

0.3

Density of X1

0.5

0.1 0

0

500

1000 Time

1500

0

2000

0

0.1

CR IP T

1 0.2

0.2

0.3

0.4

0.5

0.6

0.7

Figure 1: Eqs. (46) with β12 /2 = 0.3, β22 /2 = 0.25 and β32 /2 = 0.05. (a): paths of X1 (t), X2 (t), X3 (t), hX1 (t)i, hX2 (t)i and hX3 (t)i;

AN US

(b): probability density functions of X1 (t), X2 (t) and X3 (t) at t = 2000.

0.7

X (t) 1

X2(t)

0.6

X (t) 3

t−1∫t X (s)ds

0.5

0 1

t−1∫t0X2(s)ds

M

0.4 0.3

ED

0.2 0.1

0

PT

0

500

1000 Time

1500

2000

AC

CE

2

1.8 1.6 1.4 1.2

Density of X2

1 Density of X1

0.8 0.6 0.4 0.2 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Figure 2: Eqs. (46) with β12 /2 = 0.3, β22 /2 = 0.25 and β32 /2 = 0.2. (a): paths of X1 (t), X2 (t), X3 (t), hX1 (t)i and hX2 (t)i; (b): probability density functions of X1 (t) and X2 (t) at t = 2000.

22

ACCEPTED MANUSCRIPT

0.8 0.7 0.6 0.5 X1(t)

0.4

X3(t)

t−1∫t0X1(s)ds

0.2 0.1 0

CR IP T

X2(t)

0.3

0

500

1000 Time

1500

1

AN US

0.9

2000

0.8

Density of X1

0.7 0.6 0.5 0.4

M

0.3 0.2 0.1 0.3

0.4

ED

0 0.2

0.5

0.6

0.7

0.8

0.9

1

Figure 3: Eqs. (46) with β12 /2 = 0.3, β22 /2 = 0.44 and β32 /2 = 0.2. (a): paths of X1 (t), X2 (t), X3 (t) and hX1 (t)i; (b): probability

AC

CE

PT

density function of X1 (t) at t = 2000.

0.7 X1(t) X2(t)

0.6

X3(t)

0.5 0.4 0.3 0.2 0.1 0

0

50

100

150

Time

Figure 4: The paths of X1 (t), X2 (t) and X3 (t) in Eqs. (46) with β12 /2 = 0.65, β22 /2 = 0.44 and β32 /2 = 0.2.

23