POPULATION DYNAMICAL BEHAVIOR OF A TWO-PREDATOR ONE-PREY STOCHASTIC MODEL WITH TIME DELAY

. In this paper, the convergence of the distributions of the solutions (CDS) of a stochastic two-predator one-prey model with time delay is consid- ered. Some traditional methods that are used to study the CDS of stochastic population models without delay can not be applied to investigate the CDS of stochastic population models with delay. In this paper, we use an asymptotic approach to study the problem. By taking advantage of this approach, we show that under some simple conditions, there exist three numbers ρ 1 > ρ 2 > ρ 3 , which are represented by the coeﬃcients of the model, closely related to the CDS of our model. We prove that if ρ 1 < 1, then lim t ) 0 almost surely, i = 1 , 2 , 3; If ρ i > 1 > ρ i +1 , i = 1 , 2, then lim t → + ∞ N j ( t ) = 0 almost surely, j = i + 1 ,..., 3, and the distributions of ( N 1 ( t ) ,...,N i ( t )) T converge to a unique ergodic invariant distribution (UEID); If ρ 3 > 1, then the distributions of ( N 1 ( t ) ,N 2 ( t ) ,N 3 ( t )) T converge to a UEID. We also discuss the eﬀects of stochastic noises on the CDS and introduce several numerical examples to illustrate the theoretical results.

N i (t) = 0 almost surely, i = 1, 2, 3; If ρ i > 1 > ρ i+1 , i = 1, 2, then lim t→+∞ N j (t) = 0 almost surely, j = i + 1, ..., 3, and the distributions of (N 1 (t), ..., N i (t)) T converge to a unique ergodic invariant distribution (UEID); If ρ 3 > 1, then the distributions of (N 1 (t), N 2 (t), N 3 (t)) T converge to a UEID. We also discuss the effects of stochastic noises on the CDS and introduce several numerical examples to illustrate the theoretical results. population models cannot tend to some positive fixed point. As a result, several authors began to consider the convergence of the distributions of the solutions (CDS) of stochastic population models (see, e.g., [9,17,33,41,44,45,50,52,53,55]). By solving the corresponding Fokker-Planck equation, Pasquali [41] investigated the CDS of stochastic logistic equations. Moreover, on the basis of the Markov semigroup theory, CDS of stochastic predator-prey systems were thought over in [9,44,45,52]. In addition, according to the Lyapunov function methods (see e.g. [12]), Ji, Jiang and Shi [17] explored the CDS of a stochastic predator-prey system with modified Leslie-Gower and Holling-type II schemes; Mao [33] studied the CDS of a stochastic Lotka-Volterra facultative mutualism system; Zhao et al. [53] considered the CDS of a stochastic competitive model in a polluted environment.
A great deal of researchers have devoted themselves to single-species or two species stochastic models without time delay. However, on the one hand, it has been recognized that single-species or two species ecological models can't describe the natural phenomena accurately and many vital behaviors can only be exhibited by systems with three or more species ( [13,40,43]). On the other hand, as is well known, time delay should not be neglected because the processes of the reproductions of species are not instantaneous and time delay can reflect natural phenomena more authentically ( [11,20,23]). For instance, when a group of sheep eat some grass, the number of their population will not increase at once owing to the fact that the processes of metabolism, growth and breed all need plenty of time which cannot be neglected. Therefore, it is important and interesting to investigate the CDS of stochastic three-species models with time delay.
As is well known, when it comes to the CDS of stochastic population models, using the explicit solution of the corresponding Fokker-Planck equation is one of the traditional methods (see e.g., [41]). However, for almost all stochastic population models with time delay, it is impossible to obtain the explicit solution of the corresponding Fokker-Planck equation. Besides that, another traditional method is to use the Markov semigroup theory, see e.g. [7,9,44,45,52]. However, this method cannot be applied to deal with delay population models owing to the fact that the Markov semigroup theory needs to have some standard measures. The phase space of ODEs is a subset of R n which has a standard measure-the Lebesgue measure. The phase space of DDEs is some space of functions where it is tough to decide what is a standard measure (communicated with Professor Ryszard Rudnick). The third traditional method to study the CDS of stochastic population models is to use the Lyapunov function method (see, e.g., [17,33,55]), but allow for the fact that the Lyapunov function method requires that the solution of a stochastic model must be a Markov process while the solution of a delay stochastic model is not a Markov process. And as a result, this approach cannot be used to study delay stochastic population models either.
In this paper, consider the fact that it is a generally natural phenomenon that several predators compete for a prey, we formulate a two-predator one-prey stochastic population model with time delay in Section 2. In Section 3, we analyze the CDS of the model by using an asymptotic approach ( [28]) and discuss the results and their implications from the viewpoint of biology. In Section 4, we present some numerical simulations to illustrate the theoretical findings. In Section 5, we give some concluding remarks, and show that the approach can be applied to study the CDS of other stochastic population models with/without delay. As an example, we establish the sufficient conditions for the CDS of a three-species stochastic delay TWO-PREDATOR ONE-PREY STOCHASTIC MODEL 2515 mutualism system. We also propose some interesting problems to be explored in the future.
2. Model derivations and main results. At the first place, consider the following Lotka-Volterra delay model with two competing predators and one prey, which has been widely investigated in [1,10,13,19,36,47]: where N i (t) is the population size of the i-th species at time t (i = 1, 2, 3). r 1 > 0 stands for the growth rate of specie 1, r 2 < 0 and r 3 < 0 are the death rates of species 2 and 3 respectively. a 11 > 0, a 22 > 0, a 33 > 0 mean the inside struggle of the specie i, i = 1, 2, 3. a 12 > 0 and a 13 > 0 represent the capture rates. a 21 and a 31 are negative constants representing the growth from food. a 23 > 0 and a 32 > 0 signify the competitions between species 2 and 3. τ ij ≥ 0 denotes the delay, φ i (θ) > 0 is a continuous function on [−τ , 0] (i, j = 1, 2, 3). Based on the above model, let us allow for the environmental perturbations in addition. In most cases, the environmental perturbations can be modeled by a colored noise ( [5,39]). Several researchers (see, e.g., [39,6]) have claimed that the colored noise could be approximated by a Gaussian white noise if it was not strongly correlated, and the approximation would be quite suitable.
Several approaches have been put forward to introduce a Gaussian white noise into deterministic population models. Here, enlightened by the approach used in [4,5,15,17,18,21,25,27,29,30,31,32,33,44,45,52,53,54,55], we consider the parameter perturbation and suppose that the environmental perturbations mainly affect the growth rates of the populations, with r i → r i + σ iẆi (t), where {W i (t)} t≥0 (i = 1, 2, 3) stand for standard independent Brownian motions which are defined on a complete probability space (Ω, F, {F} t≥0 , P ), σ 2 i is the intensity of the white noise. Then we obtain the following stochastic model with delay: with initial data (1). Model (2) is approximate "to age-structured populations, with populations growth taking place in discrete time steps"( [4,6]), and as a result, we utilize the Itô calculus rather than the Stratonovich calculus.
To make it more convenient, we introduce the following notations. Let R 3 stand for all the continue functions from [−τ , 0] to R 3 + . Denote α 1 = (a 11 , a 21 , a 31 ) T , α 2 = (a 12 , a 22 , a 32 ) T , α 1 = (a 13 , a 23 , a 33 ) T and define β = (r 1 , r 2 , Before we state our main results, let us introduce some assumptions. , which mean that model (2) has a positive equilibrium state if there is no stochastic perturbations. ∆ j > 0, j = 2, 3, which mean that both species 1 and species j can coexist if both stochastic perturbations and the other predator are absent.
The following assumption is a technical assumption to make the proof work.
Assumption 3. A 12 < 0, A 13 < 0, where A ij is the symbol of the complement minor of a ij in the determinant A, i, j = 1, 2, 3.
We also suppose that R > 0. By (30) below, we have ∆ 2 /∆ 2 > ∆ 3 /∆ 3 . Note that ∆ i /∆ i measures the persistence ability of species i when its competitor is absent, i = 2, 3. Hence R > 0 means that the persistence ability of species 2 is greater that that of species 3.
(ii) If ρ 1 > 1 > ρ 2 , then species 2 and 3 go to extinction a.s., and the distribution of N 1 (t) converges weakly to a unique invariant distribution ν 1 which is ergodic: (iii) If ρ 2 > 1 > ρ 3 , then species 3 goes to extinction, and the distribution of (N 1 (t), N 2 (t)) T converges weakly to a unique invariant distribution ν 2 which is ergodic: (iv) If ρ 3 > 1, then the distribution of (N 1 (t), N 2 (t), N 3 (t)) T converges weakly to a unique invariant distribution ν 3 which is ergodic: 2 . That is to say the intensity of the stochastic noise of the prey exceeds the resist ability of the prey, hence the prey goes to extinction.
2 ) , and 2 ) measure the "help" from the prey (i.e., the capture). That is to say, the intensity of the stochastic noise of the prey is small and it could resist, but the intensities of the stochastic noises of the predators are large and the "help" from the prey is invalid. Therefore the prey is persistent and the predators go to extinction.
The second term in the left side of (6) measures the competition form the species 2, and the second term in the right side measures the "help" from the prey. That is to say, for species 2, the intensity of the stochastic noise is small and it could be persistent with the "help" of the prey. However, for species 3, the intensity of the stochastic noise and the competition from species 2 are large, and the "help" from the prey is invalid.
That is to say, for species 3, the intensity of the stochastic noise and the competition from species 2 are small, it could be persistent with the "help" of prey.
3 ) < 0, which mean that when σ 2 1 or σ 2 3 increases, species 3 tends to extinction but the stochastic noise of species 2 is favorable for the persistence of species 3. These characteristics are reasonable from the biological point of view. When the noise intensity σ 2 1 increases, the species 1 tends to extinction. Species 3 will be restricted because the food of species 3 become less which gives rise to the abatement of species 3. As a result, the noise intensity σ 2 1 is unfavorable for the reproduction of species 3. On the other hand, the species 2 and 3 compete intensely with each other, and note that the noise intensity σ 2 2 can promote the extinction of N 2 , which means that species 3 will get more food. Thus the noise intensity σ 2 2 redounds the reproduction of species 3. The effects of σ 2 i on the extinction of species 2 can be obtained similar and hence are omitted. (b) From Theorem 2.1, one can find that ρ 1 decides the extinction or not of species 1. It is not difficult to see that 3 ) = 0, which signify that the increment of σ 2 1 leads to the extinction of N 1 while σ 2 2 and σ 2 3 have nothing to do with the extinction of species 1. That is to say, only by reducing the noise intensity σ 2 1 rather than σ 2 2 or σ 2 3 can we conserve the prey population under this circumstance, which is a meaningful and useful conclusion.
The proof of Theorem 2.1 is rather long, so we divide it into three parts.
Proof. To begin with, let us consider the following stochastic system with initial condition It is not difficult to see that the coefficients of model (7) satisfy the local Lipschitz condition, therefore system (7) has a unique local solution x(t) on [0, τ e ), where τ e represents the explosion time. According to Itô's formula, we can see that is the unique positive local solution to system (2). Now let us prove τ e = +∞. To this end, let us introduce the following auxiliary model: with initial data Taking advantage of the comparison theorem for stochastic equation [16] yields that In view of Theorem 4.2 in Jiang and Shi [18], the explicit solution of system (8) is It is easy to see that y 1 (t), y 2 (t) and y 3 (t) are existent on t ≥ 0, thereby τ e = +∞.
(II)If there exist three positive constants T, λ and λ 0 such that for all t ≥ T , then z * ≥ λ/λ 0 a.s.
In the same way, we can get that We know that Using (14) in (15) gives that Multiplying (15) and (16) by −a 21 and −a 11 respectively, and adding them, one can derive that An application of (14) gives that Consequently, utilizing (18), (19), (20) and Lemma 3.5 yields that In the similar way, we can get that when Now we are in the position to prove (c). If b 2 − a 21 b 1 /a 11 < 0, applying (18), (19), (20) and Lemma 3.5 to (16), one can see that lim t→+∞ y 2 (t) = 0, a.s.
The proofs of (d) and (e) are similar to that of (b) and (c) respectively, and hence are omitted. By Hence we can obtain the following result.
To get the required assertion, let us analyze the five outcomes in Lemma 3.6 one by one. First of all, in situation (a) we have lim t→+∞ y i (t) = 0, a.s., i = 1, 2, 3.
Substituting the above two inequalities into (31) and using Lemma 3.5, we can obtain that lim , a.s. (31), (32) and (33) by t, we can derive the following equations:

Now let us turn to (iii ). Dividing
Denote m, n as the solution of the following equations:
Consequently, model (2) reduces to the following predator-prey model: which has already been investigated in [30]. Then similar to the proof of Theorem 1 in [30], the following identities can be derived: , a.s. Now we are in the position to prove (iv ). By (42), since ρ 3 = A 3 /Ã 3 > 1, we know from the arbitrariness of ε and Lemma 3.5 that Denote p, q as the solution of the following equations: Then we have According to Lemma 3.8, for arbitrarily given ε > 0, there exists T 2 > 0 such that Multiplying (36), (37) and (38) by −p, 1 and −q respectively and adding them, we can obtain that for t > T 2 Using (27) in (45) gives that Note that A 2 /Ã 2 > A 3 /Ã 3 > 1. According to the arbitrariness of ε and Lemma 3.5, we have It follows that for any sufficiently small ε, there exist T 3 and T 4 such that Substituting (47) into (36) results in that for sufficiently large t, According to the arbitrariness of ε and Lemma 3.5, we have By a 31 < 0, a 32 > 0, one can observe that for every sufficiently small ε, there exist Using (49) in (38), we can get that when t is large enough, According to the arbitrariness of ε and Lemma 3.5, we can obtain that Combining (44) with (50), one can observe that In the similar way, using (44), (48) in (37) and then combining them with (46) yield that Subsequently, substituting (51), (52) into (36) and combining them with (48), we can attain that This completes the proof.

3.2.
Proof of global attractivity. Let k i be the cofactor of the i-th diagonal element of LÂ. Then according to Kirchhoff's Matrix Tree Theorem (see, e.g., [37]), we obtain that k i > 0, i = 1, 2, 3. Define Calculating the right differential d + V (t), we can attain that where (c ij ) 3×3 =Â. By Theorem 2.3 in Li and Shuai [22], Then we have Namely that . Now we are in the position to prove that E|N i (t)−M i (t)| is uniformly continuous. In fact, it follows from (2) that That is to say, E(N i (t)), i = 1, 2, 3, are continuously differentiable functions with respect of t. On the other hand, in view of (54), and, where L 1 > 0 is a constant. Hence E(N i (t)), i = 1, 2, 3, are uniformly continuous functions. According to the well known Barbalat's conclusion ( [3]), we obtain lim t→+∞ E|N i (t) − M i (t)| = 0, i = 1, 2, 3.

Proof of stability in distribution.
Lemma 3.12. Model (2) is stable in distribution if Assumption 2 holds.
This completes the proof. Now we are in the position to prove Theorem 2.1.
Proof of Theorem 2.1. To begin with, let us show (iv). By (iv ) in Lemma 3.9, all the species are stable in time average: On the other hand, note that model (2) is stable in distribution, then the transition probability p(t, φ, ·) of N (t) converges weakly to a unique probability measure which is denoted by ν 3 . By Kolmogorov-Chapman equation, ν 3 is a invariant measure. According to Corollary 3.4.3 in [42], we can obtain that ν 3 is strong mixing. Thanks to Theorem 3.2.6 in [42], one can see that ν 3 is ergodic. By virtue of (3.3.2) in [42] and (55), one gets We are in the position to prove (iii). By (iii ) in Lemma 3.9, Consider the following predator-prey model: with initial dataÑ i (θ) = N i (θ), i = 1, 2. Since lim t→+∞ N 3 (t) = 0 a.s., the CDS of (N 1 (t), N 2 (t)) T is the same with that of (Ñ 1 (t),Ñ 2 (t)) T . Similar to the proof of Lemma 3.12, model (56) is stable in distribution. Then similar to the proof of (iv), we can show that the transition probability of (Ñ 1 (t),Ñ 2 (t)) T converges weakly to a unique invariant ergodic measure which is denoted by ν 2 , and The proof of (ii) is similar to that of (iii) and hence is omitted. (i) follows from (i ) in Lemma 3.9, and this completes the proof.
Traditional approaches to investigate the CDS of stochastic population models are to solve the corresponding Fokker-Planck equations ( [41]), or to use the Markov semigroup theory ( [44,45,52]), or to use the Lyapunov function method ( [17,33,55]). However, for most stochastic delay models, these three approaches can not be used. In this paper, we use an asymptotic method ( [28]) to study the CDS problems, which can avoid the difficulties of these three traditional approaches. The algorithm of this method can be described as follows: consider the CDS problem of the following stochastic Lotka-Volterra model with time delay: where N i (θ) = φ i (θ), θ ∈ [−τ , 0].
Step 1. Establish sufficient conditions for stability in time average of all species in model (57). For most two-species models, these sufficient conditions have been obtained in the literature (see e.g., [26,30,48]). However, for general n− species model (57), the sufficient conditions are difficult to establish. One possible way to obtain the sufficient conditions is to apply Itô's formula to ln N i and then using Lemma 3.5 and the comparison theorem for stochastic equation to estimate lim inf t→+∞ N i (t) and lim sup t→+∞ N i (t) .
Step 2. Prove that model (57) is globally attractive. There are two steps to show this. First, to give some conditions under which model (57) satisfies (53) (Liu [24] (Lemma 2.2) established some sufficient conditions (not necessary) for this, the author showed that if a ii > 4 3 j=1,j =i |a ij |, then model (57) satisfies (53)). The second step is to applied Itô's formula to and k i is the cofactor of the i-th diagonal element of the corresponding matrix (see the proof of Lemma 3.11), and then use (53) to complete the proof of global attractivity.
Step 3. Prove that model (57) is stability in distribution under the conditions in Step 2. The proof is a slight modification of that in Lemma 3.12 (i.e., R 3 + changes to R n + in the proof).