The Asymptotic Behavior of a Stochastic Predator-Prey System with Holling II Functional Response

Abstract

We discuss a stochastic predator-prey system with Holling II functional response. First, we show that this system has a unique positive solution as this is essential in any population dynamics model. Then, we deduce the conditions that there is a stationary distribution of the system, which implies that the system is permanent. At last, we give the conditions for the system that is going to be extinct.

1. Introduction

One of the most popular predator-prey model is the one with Michaelis-Menten type (or Holling Type II) functional response [1, 2]: 𝑥̇(𝑡)(𝑡)(𝑡)−(𝑡)=𝑥𝑎−𝑏𝑥𝛼𝑦(𝑡),𝑦̇(𝑡)(𝑡)(𝑡)1+𝛽𝑥=𝑦−𝑒+𝑘𝛼𝑥(𝑡),1+𝛽𝑥 (1.1) where 𝑥(𝑡) and 𝑦(𝑡) are the population densities of prey and predator at time 𝑡, respectively. The constants 𝑎, 𝑏/𝑎, 𝛼, 𝛽, 𝑒, and 𝑘 are positive constants that stand for prey intrinsic growth rate, carrying capacity, the maximum ingestion rate, half-saturation constant, predator death rate, and the conversion factor, respectively. This model exhibits the well-known but highly controversial “paradox of enrichment” observed by Hairston et al. [3] and by Rosenzweig [4] which is rarely reported in nature. It is very important to study the existence and asymptotical stability of equilibria and limit cycle for autonomous predator-prey systems with Holling II functional response. If 𝑘𝑎𝛼𝛽>𝑎𝑒𝛽²+𝑘𝑏𝛼+𝑏𝑒𝛽, then system (1.1) has a unique limit cycle which is stable. If 𝑎𝑘𝛼>𝑎𝑒𝛽+𝑏𝑒, then system (1.1) has a unique positive equilibrium: 𝑥^∗=𝑒𝑘𝛼−𝑒𝛽,𝑦^∗=𝑘𝛼𝑎𝑘𝛼−𝑎𝑒𝛽−𝑏𝑒𝑘𝛼−𝑒𝛽², (1.2) which is a stable node or focus (see [5]).

However, countless organisms live in seasonally or diurnally forced environments. Hence, authors considered models with periodic ecological parameters or perturbations. For example, Liu and Chen [6] introduced periodic constant impulsive immigration of predator into system (1.1) and gave conditions for the system to be extinct and permanence, respectively. Zhang and Chen [7] studied a Holling II functional response food chain model with impulsive perturbations. Zhang et al. [8] further considered system (1.1) with periodic constant impulsive immigration of predator and periodic variation in the intrinsic growth rate of the prey.

On the other hand, the white noise is always present, and we cannot omit the influence of the white noise to the system. May [9] pointed out that due to continuous fluctuation in the environment, the birth rates, death rates, carrying capacity, competition coefficients, and all other parameters involved with the model exhibit random fluctuation to a great lesser extent, and as a result the equilibrium population distribution never attains a steady value, but fluctuates randomly around some average value. Many authors studied the effect of the stochastic perturbation to the predator-prey system with different functional responses, such as [10–14]. Therefore, in this paper, we also introduce stochastic perturbation system (1.1) and obtain the following stochastic system: (𝑡)(𝑡)(𝑡)−(𝑡)𝑑𝑥=𝑥𝑎−𝑏𝑥𝛼𝑦1+𝛽𝑥(𝑡)𝑑𝑡+𝜎₁𝑥(𝑡)𝑑𝐵₁(𝑡),(𝑡)(𝑡)(𝑡)𝑑𝑦=𝑦−𝑒+𝑘𝛼𝑥1+𝛽𝑥(𝑡)𝑑𝑡+𝜎₂𝑦(𝑡)𝑑𝐵₂(𝑡), (1.3) where 𝐵₁(𝑡) and 𝐵₂(𝑡) are mutually independent Brownian motion with 𝐵₁(0)=𝐵₂(0)=0, and 𝜎²₁, 𝜎²₂ are intensities of the white noise.

The aim of this paper is to discuss the long time behavior of system (1.3). As the deterministic population models, we are also interested in the permanence and extinction of the system. The global stability of the positive equilibrium means that the system is permanence. But, for the stochastic system, there is no positive equilibrium. Hence, it is impossible that the solution of system (1.3) will tend to a fixed point. In this paper, we show that there is a stationary distribution of system (1.3) mainly according to the theory of Has’meminskii [15], if the white noise is small. While if the white noise is large, based on the techniques developed in [16, 17], we prove that the predator population will die out a.s. and the prey population will either extinct or its distribution converges to a probability measure. It does not happen that both the prey population and the predator population in system (1.3) will die out, which is brought by large white noise, such as weather, epidemic disease. From this point, we say the stochastic model is more realistic than the deterministic model.

The rest of this paper is organized as follows. In Section 2, we show that there is a unique nonnegative solution of system (1.3). In Section 3, we show that there is a stationary distribution under small white noise. While in Section 4, we consider the situation when the white noise is large. We prove that the system will be extinct. Finally, we give an appendix containing the stationary distribution theory used in Section 3.

2. Existence and Uniqueness of the Nonnegative Solution

To investigate the dynamical behavior, the first concern is the global existence of the solutions. Hence in this section we show that the solution of system (1.3) is global and nonnegative. It is not difficult to check the uniqueness and global existence of solutions if the coefficients of the equation satisfy the linear growth condition and local Lipschitz condition (cf. [18]). However, the coefficients of system (1.3) do not satisfy the linear growth condition, but locally Lipschitz continuous, so the solution of system (1.3) may explode at a finite time. In this section, by changing variables, we first show that system (1.3) has a local solution, then show that this solution is global.

Theorem 2.1. For any initial value (𝑥(0),𝑦(0))∈𝑅²₊, there is a unique solution (𝑥(𝑡),𝑦(𝑡)) of system (1.3) on 𝑡≥0, and the solution will remain in 𝑅²₊ with probability 1.

Proof. First, consider the following system, by changing variables, 𝑥(𝑡)=𝑒^𝑢(𝑡), 𝑦(𝑡)=𝑒^𝑣(𝑡), (𝑡)=𝜎𝑑𝑢𝑎−²¹2−𝑏𝑒^𝑢(𝑡)−𝛼𝑒^𝑣(𝑡)1+𝛽𝑒^𝑢(𝑡)𝑑𝑡+𝜎₁𝑑𝐵₁(𝑡),(𝑡)=𝜎𝑑𝑣−𝑒−²²2+𝑘𝛼𝑒^𝑢(𝑡)1+𝛽𝑒^𝑢(𝑡)𝑑𝑡+𝜎₂𝑑𝐵₂(𝑡). (2.1) It is clear that the coefficients of system (2.1) are locally Lipschitz continuous for the given initial value (log𝑥(0),log𝑦(0))∈𝑅² there is a unique local solution (𝑢(𝑡),𝑣(𝑡)) on 𝑡∈[0,𝜏_𝑒), where 𝜏_𝑒 is the explosion time (see [18]). Hence, by Itô formula, we know (𝑒^𝑢(𝑡),𝑒^𝑣(𝑡)),𝑡∈[0,𝜏_𝑒) is a unique positive local solution of system (1.3). To show that this solution is global, we need to show that 𝜏_𝑒=∞ a.s. Let 𝑚₀≥1 be sufficiently large so that 𝑥(0),𝑦(0) all lie within the interval [1/𝑚₀,𝑚₀]. For each integer 𝑚≥𝑚₀, define the stopping time: 𝜏_𝑚=inf𝑡∈0,𝜏_𝑒𝑥(𝑡)(𝑡≤1∶min,𝑦)𝑚ormax𝑥(𝑡),𝑦(𝑡)≥𝑚, (2.2) Where, throughout this paper, we set inf∅=∞ (as usual ∅ denotes the empty set). Clearly, 𝜏_𝑚 is increasing as 𝑚→∞. Set 𝜏_∞=lim_𝑚→∞𝜏_𝑚, whence 𝜏_∞≤𝜏_𝑒 a.s. If we can show that 𝜏_∞=∞ a.s., then 𝜏_𝑒=∞ and (𝑥(𝑡),𝑦(𝑡))∈𝑅²₊ a.s. for all 𝑡≥0. In other words, to complete the proof all we need to show is that 𝜏_∞=∞ a.s. If this statement is false, then there is a pair of constants 𝑇>0 and 𝜖∈(0,1) such that 𝑃𝜏_∞≤𝑇>𝜖. (2.3) Hence there is an integer 𝑚₁≥𝑚₀ such that 𝑃𝜏_𝑚≤𝑇≥𝜖∀𝑚≥𝑚₁. (2.4) Define a 𝐶²-function 𝑉∶𝑅²₊→𝑅₊ by 𝑉=𝑥𝑥,𝑦𝑥−𝑐−𝑐log𝑐+1𝑘𝑦−1−log𝑦, (2.5) where 𝑐 is a positive constant to be determined later. The nonnegativity of this function can be seen from 𝑢−1−log𝑢≥0,forall𝑢>0. Using Itô's formula, we get 𝑑𝑉∶=𝐿𝑉𝑑𝑡+𝜎₁(𝑥−𝑐)𝑑𝐵₁(𝑡)+𝜎₂𝑘𝑦−1𝑑𝐵₂(𝑡), (2.6) where 𝐿𝑉=(𝑥−𝑐)𝑎−𝑏𝑥−𝛼𝑦+1+𝛽𝑥𝑐𝜎²₁2+1𝑘𝑦−1−𝑒+𝑘𝛼𝑦+𝜎1+𝛽𝑥2𝑘²²=−𝑎𝑐+𝑐𝜎²₁2+𝑒𝑘+𝜎²₂+()𝑒2𝑘𝑎+𝑏𝑐𝑥−𝑘𝑦−𝑏𝑥²+𝛼𝑐𝑦1+𝛽𝑥≤−𝑎𝑐+𝑐𝜎²₁2+𝑒𝑘+𝜎²₂+()2𝑘𝑎+𝑏𝑐𝑥−𝑏𝑥²−𝑒𝑘−𝛼𝑐𝑦. (2.7) Choose 𝑐=𝑒/𝛼𝑘 such that 𝑒/𝑘−𝛼𝑐=0, then 𝐿𝑉≤−𝑎𝑐+𝑐𝜎²₁2+𝑒𝑘+𝜎²₂+()2𝑘𝑎+𝑏𝑐𝑥−𝑏𝑥²≤𝐾, (2.8) where 𝐾 is a positive constant. Therefore 𝑑𝑉𝑥^𝜏𝑚₀^∧𝑇(𝑡),𝑦(𝑡)≤^𝜏𝑚₀^∧𝑇+𝐾𝑑𝑡_{𝜏𝑚0∧𝑇}𝜎₁(𝑥(𝑠)−𝑐)𝑑𝐵₁(𝑠)+𝜎₂𝑘𝑦(𝑠)−1𝑑𝐵₂(𝑠), (2.9) which implies that, 𝐸𝑉𝑥𝜏_𝑚∧𝑇,𝑦𝜏_𝑚∧𝑇≤𝑉𝑥(0),𝑦(0)+𝐸^𝜏𝑚₀^∧𝑇𝐾𝑑𝑡≤𝑉𝑥(0),𝑦(0)+𝐾𝑇. (2.10) Set Ω_𝑚={𝜏_𝑚≤𝑇} for 𝑚≥𝑚₁, then by (2.4), we know that 𝑃(Ω_𝑚)≥𝜖. Note that for every 𝜔∈Ω_𝑚, there is at least one of 𝑥(𝜏_𝑚,𝜔) and 𝑦(𝜏_𝑚,𝜔) equals either 𝑚 or 1/𝑚, then 𝑉𝑥𝜏_𝑚,𝑦𝜏_𝑚≥𝑚𝑚−𝑐−𝑐log𝑐∧1𝑚(∧1−𝑐+𝑐log𝑐𝑚)𝑘∧1𝑚−1−log𝑚𝑘1𝑚−1+log𝑚. (2.11) It then follows from (2.4) and (2.10) that 𝑉𝑥(0),𝑦(0)+𝐾𝑇≥𝐸1_Ω𝑚(𝜔)𝑉𝑥𝜏_𝑚,𝑦𝜏_𝑚𝑚≥𝜖𝑚−𝑐−𝑐log𝑐∧1𝑚()∧1−𝑐+𝑐log𝑐𝑚𝑘∧1𝑚−1−log𝑚𝑘1𝑚−1+log𝑚, (2.12) where 1_Ω𝑚(𝜔) is the indicator function of Ω_𝑚. Letting 𝑚→∞ leads to the contradiction that ∞>𝑉(𝑥(0),𝑦(0))+𝐾𝑇=∞. So we must therefore have 𝜏_∞=∞ a.s.

3. Permanence

There is no equilibrium of system (1.3). Hence we cannot show the permanence of the system by proving the stability of the positive equilibrium as the deterministic system. In this section we show that there is a stationary distribution of system (1.3).

Remark 3.1. Theorem 2.1 shows that there exists a unique positive solution (𝑥(𝑡),𝑦(𝑡)) of system (1.3) with any initial value (𝑥(0),𝑦(0))∈𝑅²₊. From the proof of Theorem 2.1, we obtain that 𝐿𝑉≤𝐾. Define 𝑉=𝑉+𝐾, then 𝐿𝑉≤𝑉, and it is clear that 𝑉_𝑅=inf_{(𝑥,𝑦)∈𝑅}²_+⧵𝐷𝑘𝑉(𝑥,𝑦)→∞ as 𝑘→∞, where 𝐷_𝑘=(1/𝑘,𝑘)×(1/𝑘,𝑘). Hence by Remark 2 of Theorem 4.1 of Has’meminskii, 1980, page 86 in [15], we obtain that the solution (𝑥(𝑡),𝑦(𝑡)) is a homogeneous Markov process in 𝑅²₊.

Theorem 3.2. If 𝑎𝑒𝛽+𝑏𝑒<𝑎𝑘𝛼<𝑎𝑒𝛽+𝑏𝑘𝛼/𝛽 and 𝜎₁>0,𝜎₂>0 such that 𝜎²₂<𝑘𝛼𝑥^∗/(1+𝛽𝑥^∗) and 12+𝑙₂𝑥^∗𝑥^∗𝜎²₁+1+𝛽𝑥^∗2+𝑙₂𝑦^∗𝑘𝑦^∗𝜎²₂𝑘1<min2𝑏−𝛼𝛽𝑦^∗1+𝛽𝑥^∗𝑥^∗2,𝑙₂2𝛼𝑥^∗𝑘1+𝛽𝑥^∗−𝜎²₂𝑘²𝑦^∗2, (3.1) where (𝑥^∗,𝑦^∗) is the positive equilibrium of system (1.1) and 𝑙₂ is defined as in the proof. Then system (1.3) has a stationary ergodic solution.

Proof. Since 𝑎𝑘𝛼>𝑎𝑒𝛽+𝑏𝑒, then there is a positive equilibrium (𝑥^∗,𝑦^∗) of system (1.1), and 𝑎=𝑏𝑥^∗+𝛼𝑦^∗1+𝛽𝑥^∗,𝑒=𝑘𝛼𝑥^∗1+𝛽𝑥^∗. (3.2) Let 𝑉₁𝑥,𝑦=𝑥−𝑥^∗−𝑥^∗𝑥log𝑥^∗+𝑙₁𝑦−𝑦^∗−𝑦^∗𝑦log𝑦^∗, (3.3) where 𝑙₁ is a positive constant to be determined later. Let 𝐿 be the generating operator of system (1.3). Then 𝐿𝑉₁=𝑥−𝑥^∗𝑎−𝑏𝑥−𝛼𝑦+𝑥1+𝛽𝑥^∗𝜎²₁2+𝑙₁𝑦−𝑦^∗−𝑒+𝑘𝛼𝑥+𝑙1+𝛽𝑥¹𝑦^∗𝜎²₂2=𝑥−𝑥^∗−𝑏𝑥−𝑥^∗−𝛼1+𝛽𝑥𝑦−𝑦^∗+𝛼𝛽𝑦^∗1+𝛽𝑥^∗1+𝛽𝑥𝑥−𝑥^∗+𝑥^∗𝜎²₁2+𝑙₁𝑦−𝑦^∗𝑘𝛼𝑥−𝑥^∗1+𝛽𝑥^∗+𝑙1+𝛽𝑥¹𝑦^∗𝜎²₂2≤−𝑏−𝛼𝛽𝑦^∗1+𝛽𝑥^∗𝑥−𝑥^∗2−𝛼𝑙1+𝛽𝑥1−¹𝑘1+𝛽𝑥^∗𝑥−𝑥^∗𝑦−𝑦^∗+𝑥^∗𝜎²₁2+𝑙₁𝑦^∗𝜎²₂2. (3.4) Choose 𝑙₁=(1+𝛽𝑥^∗)/𝑘 such that 1−𝑙₁𝑘/(1+𝛽𝑥^∗)=0 and yields 𝐿𝑉₁≤−𝑏−𝛼𝛽𝑦^∗1+𝛽𝑥^∗𝑥−𝑥^∗2+𝑥^∗𝜎²₁2+1+𝛽𝑥^∗𝑦^∗𝜎²₂.2𝑘 (3.5) Let 𝑉₂=1𝑥,𝑦2𝑥−𝑥^∗+1𝑘𝑦−𝑦^∗2. (3.6) Note that 𝑑𝑥−𝑥^∗+1𝑘𝑦−𝑦^∗=𝑎𝑥−𝑏𝑥²−𝑒𝑘𝑦𝑑𝑡+𝜎₁𝑥𝑑𝐵₁(𝑡)+𝜎₂𝑘𝑦𝑑𝐵₂(𝑡)=−𝑏𝑥𝑥−𝑥^∗2𝑦+𝛼^∗𝑥−𝑥^∗−𝑥^∗𝑦−𝑦^∗1+𝛽𝑥^∗𝑑𝑡+𝜎₁𝑥𝑑𝐵₁(𝑡)+𝜎₂𝑘𝑦𝑑𝐵₂(𝑡), (3.7) then 𝐿𝑉₂=𝑥−𝑥^∗+1𝑘𝑦−𝑦^∗−𝑏𝑥𝑥−𝑥^∗2𝑦+𝛼^∗𝑥−𝑥^∗−𝑥^∗𝑦−𝑦^∗1+𝛽𝑥^∗+𝜎²₁2𝑥²+𝜎²₂2𝑘²𝑦²=−𝑏𝑥𝑥−𝑥^∗2+𝛼𝑦^∗1+𝛽𝑥^∗+𝑏𝑦^∗𝑘𝑥−𝑥^∗2−𝑏𝑘𝑥−𝑥^∗2𝑦−𝛼𝑥^∗𝑘1+𝛽𝑥^∗𝑦−𝑦^∗2+𝛼𝑦^∗𝑘1+𝛽𝑥^∗−𝛼𝑥^∗1+𝛽𝑥^∗−𝑏𝑥^∗𝑘𝑥−𝑥^∗𝑦−𝑦^∗+𝜎²₁2𝑥²+𝜎²₂2𝑘²𝑦²≤𝛼𝑦^∗1+𝛽𝑥^∗+𝑏𝑦^∗𝑘+𝜎²₁𝑥−𝑥^∗2−𝛼𝑥^∗𝑘1+𝛽𝑥^∗−𝜎²₂𝑘²𝑦−𝑦^∗2+𝛼𝑦^∗𝑘1+𝛽𝑥^∗−𝛼𝑥^∗1+𝛽𝑥^∗−𝑏𝑥^∗𝑘𝑥−𝑥^∗𝑦−𝑦^∗+𝜎²₁𝑥^∗2+𝜎²₂𝑘²𝑦^∗2, (3.8) where 𝐿 is also the generating operator of system (1.3). Note that 𝛼𝑦^∗𝑘1+𝛽𝑥^∗−𝛼𝑥^∗1+𝛽𝑥^∗−𝑏𝑥^∗𝑘𝑥−𝑥^∗𝑦−𝑦^∗≤𝛼𝑦^∗/𝑘1+𝛽𝑥^∗−𝛼𝑥^∗/1+𝛽𝑥^∗−𝑏𝑥^∗/𝑘²2𝛼𝑥^∗/𝑘1+𝛽𝑥^∗−𝜎²₂/𝑘²𝑥−𝑥^∗2+12𝛼𝑥^∗𝑘1+𝛽𝑥^∗−𝜎²₂𝑘²𝑦−𝑦^∗2∶=𝛿𝑥−𝑥^∗2+12𝛼𝑥^∗𝑘1+𝛽𝑥^∗−𝜎²₂𝑘²𝑦−𝑦^∗2, (3.9) then 𝐿𝑉₂≤𝛼𝑦^∗1+𝛽𝑥^∗+𝑏𝑦^∗𝑘+𝜎²₁+𝛿𝑥−𝑥^∗2−12𝛼𝑥^∗𝑘1+𝛽𝑥^∗−𝜎²₂𝑘²𝑦−𝑦^∗2+𝜎²₁𝑥^∗2+𝜎²₂𝑘²𝑦^∗2. (3.10) Now define 𝑉𝑥,𝑦=𝑉₁𝑥,𝑦+𝑙₂𝑉₂𝑥,𝑦, (3.11) where 𝑙₂ is a positive constant to be determined later. Then 𝐿𝑉≤−𝑏−𝛼𝛽𝑦^∗1+𝛽𝑥^∗−𝑙₂𝛼𝑦^∗1+𝛽𝑥^∗+𝑏𝑦^∗𝑘+𝜎²₁+𝛿𝑥−𝑥^∗2−𝑙₂2𝛼𝑥^∗𝑘1+𝛽𝑥^∗−𝜎²₂𝑘²𝑦−𝑦^∗2+12+𝑙₂𝑥^∗𝑥^∗𝜎²₁+1+𝛽𝑥^∗2+𝑙₂𝑦^∗𝑘𝑦^∗𝜎²₂𝑘. (3.12) Choose 𝑙₂>0 such that (𝑏−𝛼𝛽𝑦^∗/(1+𝛽𝑥^∗)−𝑙₂(𝛼𝑦^∗/(1+𝛽𝑥^∗)+𝑏𝑦^∗/𝑘+𝜎²₁+𝛿))=(1/2)(𝑏−𝛼𝛽𝑦^∗/(1+𝛽𝑥^∗)), then it follows from (3.12) that 1𝐿𝑉≤−2𝑏−𝛼𝛽𝑦^∗1+𝛽𝑥^∗𝑥−𝑥^∗2−𝑙₂2𝛼𝑥^∗𝑘1+𝛽𝑥^∗−𝜎²₂𝑘²𝑦−𝑦^∗2+12+𝑙₂𝑥^∗𝑥^∗𝜎²₁+1+𝛽𝑥^∗2+𝑙₂𝑦^∗𝑘𝑦^∗𝜎²₂𝑘. (3.13) Note that 12+𝑙₂𝑥^∗𝑥^∗𝜎²₁+1+𝛽𝑥^∗2+𝑙₂𝑦^∗𝑘𝑦^∗𝜎²₂𝑘1<min2𝑏−𝛼𝛽𝑦^∗1+𝛽𝑥^∗𝑥^∗2,𝑙₂2𝛼𝑥^∗𝑘1+𝛽𝑥^∗−𝜎²₂𝑘²𝑦^∗2, (3.14) then the ellipsoid −12𝑏−𝛼𝛽𝑦^∗1+𝛽𝑥^∗𝑥−𝑥^∗2−𝑙₂2𝛼𝑥^∗𝑘1+𝛽𝑥^∗−𝜎²₂𝑘²𝑦−𝑦^∗2+12+𝑙₂𝑥^∗𝑥^∗𝜎²₁+1+𝛽𝑥^∗2+𝑙₂𝑦^∗𝑘𝑦^∗𝜎²₂𝑘=0 (3.15) lies entirely in 𝑅²₊. We can take 𝑈 to be a neighborhood of the ellipsoid with 𝑈⊆𝐸_𝑙=𝑅²₊, so that for (𝑥,𝑦)∈𝑈⧵𝐸_𝑙, 𝐿𝑉≤−𝐶 (𝐶 is a positive constant), which implies condition (B.2) in Lemma A.1 is satisfied. Hence the solution (𝑥(𝑡),𝑦(𝑡)) is recurrent in the domain 𝑈, which together with Lemma A.3 and Remark 3.1 implies that (𝑥(𝑡),𝑦(𝑡)) is recurrent in any bounded domain 𝐷⊂𝑅²₊. Besides, for all 𝐷, there is an 𝑀=min{𝜎²₁𝑥²,𝜎²₂𝑦²,(𝑥,𝑦)∈𝐷}>0 such that ²_𝑖,𝑗=1𝜆_𝑖𝑗𝜉_𝑖𝜉_𝑗=𝜎²₁𝑥²𝜉²₁+𝜎²₂𝑦²𝜉²₂≥𝑀𝜉²all𝑥∈𝐷,𝜉∈𝑅², (3.16) which implies that condition (B.1) is also satisfied. Therefore, system (1.3) has a stable a stationary distribution 𝜇(⋅) and it is ergodic.

Note that 𝑑𝑥^𝑝=𝑝𝑥^𝑝𝑎−𝑏𝑥−𝛼𝑦1+𝛽𝑥𝑑𝑡+𝑝𝜎₁𝑥^𝑝𝑑𝐵₁(𝑡)+12𝑝𝑝−1𝜎²₁𝑥^𝑝𝑑𝑡≤𝑝𝑥^𝑝𝑝𝑎+2𝜎²₁−𝑏𝑥𝑑𝑡+𝑝𝜎₁𝑥^𝑝𝑑𝐵₁(𝑡), (3.17) then 𝑑𝐸𝑥^𝑝𝑑𝑡≤𝑝𝑎+𝑝𝜎²₁2𝐸𝑥^𝑝−𝑏𝐸𝑥^𝑝+1≤𝑝𝑎+𝑝𝜎²₁2𝐸𝑥^𝑝−𝑏𝐸[𝑥^𝑝]^{(𝑝+1)/𝑝}. (3.18) Hence by comparison theorem, we get limsup_𝑡→∞𝐸𝑥^𝑝(𝑡)≤𝑎+𝑝𝜎²₁/2𝑏^𝑝, (3.19) by which together with the continuity of 𝐸[𝑥^𝑝(𝑡)], we have that there exists a positive constant 𝐾=𝐾(𝑝) such that 𝐸𝑥^𝑝(𝑡)≤𝐾𝑝. (3.20)

By Doob's martingale inequality, together with the (3.20), for 𝛿>0, we have 𝑃𝜔∶sup_{(𝑛−1)𝛿≤𝑡≤𝑛𝛿}𝑥(𝑡)𝑡≤𝐸𝑥>𝛿^𝑝(𝑛𝛿)(𝑛𝛿)^𝑝𝛿≤𝐾𝑝𝑛^𝑝𝛿^𝑝+1,𝑝>1. (3.21) In view of the well-known Borel-Cantelli lemma, we see that for almost all 𝜔∈Ω, sup_{(𝑛−1)𝛿≤𝑡≤𝑛𝛿}𝑥(𝑡)𝑡≤𝛿 (3.22) holds for all but finitely many 𝑛. Hence there exists an 𝑛₀(𝜔), for all 𝜔∈Ω excluding a 𝑃-null set, for which (3.22) holds whenever 𝑛≥𝑛₀. Consequently, letting 𝛿→0, we have, for almost all 𝜔, lim_𝑡→∞𝑥(𝑡)𝑡=0. (3.23)

By the ergodic property, for any given constant 𝑚>0, we have lim_𝑡→∞1𝑡^𝑡₀𝑥^𝑝(𝑠)∧𝑚𝑑𝑠=_𝑅2+𝑧^𝑝₁∧𝑚𝜇𝑑𝑧₁,𝑑𝑧₂,a.s. (3.24) On the other hand, by dominated convergence theorem and (3.20), we get 𝐸lim_𝑡→∞1𝑡^𝑡₀𝑥^𝑝(𝑠)∧𝑚𝑑𝑠=lim_𝑡→∞1𝑡^𝑡₀𝐸𝑥^𝑝(𝑠)∧𝑚𝑑𝑠≤𝐾𝑝, (3.25) which together with (3.24) implies _𝑅2+𝑧^𝑝₁∧𝑚𝜇𝑑𝑧₁,𝑑𝑧₂≤𝐾𝑝. (3.26) Letting 𝑚→∞, we get _𝑅2+𝑧^𝑝₁𝜇𝑑𝑧₁,𝑑𝑧₂≤𝐾𝑝. (3.27) That is to say, the function 𝑓₁(𝑧)=𝑧^𝑝₁ is integrable with respect to the measure 𝜇. Therefore, by ergodicity property again, we get lim_𝑡→∞1𝑡^𝑡₀𝑥^𝑝(𝑠)𝑑𝑠=_𝑅2+𝑧^𝑝₁𝜇𝑑𝑧₁,𝑑𝑧₂,a.s. (3.28) Besides, 𝑥(𝑡)−𝑥(0)𝑡=𝑎𝑡^𝑡₀𝑥(𝑠)𝑏𝑑𝑠−𝑡^𝑡₀𝑥²(𝑠)𝛼𝑑𝑠−𝑡^𝑡₀𝑥(𝑠)𝑦(𝑠)(𝑠)𝜎1+𝛽𝑥𝑑𝑠+¹𝑡^𝑡₀𝑥(𝑠)𝑑𝐵₁(𝑠). (3.29) Let 𝑀₁(𝑡)=∫^𝑡₀𝑥(𝑠)𝑑𝐵₁(𝑠) which is a martingale with 𝑀₁(0)=0 and limsup_𝑡→∞𝑀₁,𝑀_1𝑡𝑡=lim_𝑡→∞1𝑡^𝑡₀𝑥²₁(𝑠)𝑑𝑠=_𝑅2+𝑧²₁𝜇𝑑𝑧₁,𝑑𝑧₂<∞, (3.30) then by the strong law of large numbers, we get lim_𝑡→∞𝑀₁(𝑡)𝑡=lim_𝑡→∞1𝑡^𝑡₀𝑥₁(𝑠)𝑑𝐵₁(𝑠)=0, (3.31) which together with (3.23) and (3.28) implies (3.29) that lim_𝑡→∞1𝑡^𝑡₀𝑥(𝑠)𝑦(𝑠)(𝑠)𝑎1+𝛽𝑥𝑑𝑠=𝛼_𝑅²₊𝑧₁𝜇𝑑𝑧₁,𝑑𝑧₂−𝑏𝛼_𝑅²₊𝑧²₁𝜇𝑑𝑧₁,𝑑𝑧₂. (3.32)

Hence from these arguments, we get the following result.

Theorem 3.3. Assume the same conditions as in Theorem 3.2. Then one has lim_𝑡→∞1𝑡^𝑡₀𝑥^𝑝(𝑠)𝑑𝑠=_𝑅2+𝑧^𝑝₁𝜇𝑑𝑧₁,𝑑𝑧₂,a.s.,lim_𝑡→∞1𝑡^𝑡₀𝑥(𝑠)𝑦(𝑠)(𝑠)𝑎1+𝛽𝑥𝑑𝑠=𝛼_𝑅²₊𝑧₁𝜇𝑑𝑧₁,𝑑𝑧₂−𝑏𝛼_𝑅²₊𝑧²₁𝜇𝑑𝑧₁,𝑑𝑧₂. (3.33)

4. Extinction

In this section, we show the situation when system (1.3) will be extinct.

Case 1. (𝑎<𝜎²₁/2).
Obviously, 𝑑𝑥≤𝑥(𝑎−𝑏𝑥)𝑑𝑡+𝜎₁𝑥𝑑𝐵₁(𝑡), (4.1) then when 𝑎<𝜎²₁/2, lim_𝑡→∞𝑥(𝑡)=0,a.s. (4.2) That is to say, for all 0<𝜖₁<𝑒+𝜎²₂/2, there exist 𝑇₁=𝑇₁(𝜔) and a set Ω_𝜖1 such that 𝑃(Ω_𝜖1)>1−𝜖₁ and 𝑘𝛼𝑥(𝑡)≤𝜖₁ for 𝑡≥𝑇₁ and 𝜔∈Ω_𝜖1. Then −𝑒𝑦(𝑡)𝑑𝑡+𝜎₂𝑦(𝑡)𝑑𝐵₂(𝑡)≤𝑑𝑦(𝑡)≤𝑦(𝑡)−𝑒+𝜖₁𝑑𝑡+𝜎₂𝑦(𝑡)𝑑𝐵₂(𝑡), (4.3) and so lim_𝑡→∞𝑦(𝑡)=0,a.s. (4.4)

Case 2. (𝑎>𝜎²₁/2,𝑒+𝜎²₂/2>𝑘𝛼/𝛽).
Note that 𝑑𝑦(𝑡)≤𝑦(𝑡)−𝑒+𝑘𝛼𝛽𝑑𝑡+𝜎₂𝑦(𝑡)𝑑𝐵₂(𝑡), (4.5) then if 𝑒+𝜎²₂/2>𝑘𝛼/𝛽, we have lim_𝑡→∞𝑦(𝑡)=0,a.s.. (4.6) In this situation, for all 0<𝜖₂<𝑎−𝜎²₁/2, there exist 𝑇₂=𝑇₂(𝜔) and a set Ω_𝜖2 such that 𝑃(Ω_𝜖2)>1−𝜖₂ and 𝛼𝑦(𝑡)≤𝜖₂ for 𝑡≥𝑇₂ and 𝜔∈Ω_𝜖2. Then 𝑥(𝑡)𝑎−𝑏𝑥(𝑡)−𝜖₂𝑑𝑡+𝜎₁𝑥(𝑡)𝑑𝐵₁(𝑡)≤𝑑𝑥(𝑡)≤𝑥(𝑡)(𝑎−𝑏𝑥(𝑡))𝑑𝑡+𝜎₁𝑥(𝑡)𝑑𝐵₁(𝑡),𝜎𝑎−²¹2−𝑏𝑥(𝑡)−𝜖₂𝑑𝑡+𝜎₁𝑑𝐵₁(𝑡)(𝑡)≤𝜎≤𝑑𝑥𝑎−²¹2−𝑏𝑥(𝑡)𝑑𝑡+𝜎₁𝑑𝐵₁(𝑡). (4.7)(4.8) Consider the following equation: (𝑡)=𝜎𝑑Φ𝑎−²¹2−𝑏𝑒^Φ(𝑡)𝑑𝑡+𝜎₁𝑑𝐵₁(𝑡). (4.9) If 𝑎>𝜎²₁/2, (4.9) has the density 𝑔_∗(𝜁) such that 12𝜎²₁𝑔^′_∗(𝜁)=𝑎₁−𝜎²₁2−𝑏₁₁𝑒^𝜁𝑔_∗(𝜁). (4.10) Therefore from (4.8) and the arbitrary of 𝜖₂, we get that the distribution of log𝑥(𝑡) converges weakly to the probability measure with density 𝑔_∗. Thus, from (4.10), we obtain that the distribution of 𝑥(𝑡) converges weakly to the probability measure with density 𝑓^∗(𝜁)=𝐶₀𝜁^{2(𝑎−𝜎21/2)/𝜎21−1}𝑒^{−2𝑏𝜁/𝜎21}, where 𝐶₀=(2𝑏/𝜎²₁)^{2(𝑎−𝜎21/2)/𝜎21}/Γ(2(𝑎−𝜎²₁/2)/𝜎²₁). Besides, from the ergodic theorem and (4.10), it follows that lim_𝑡→∞1𝑡^𝑡₀𝑥(𝑠)𝑑𝑠=^∞_−∞𝑒^𝜁𝑔_∗(𝜁)𝑑𝜁=^∞_−∞𝑎−𝜎²₁/2𝑏𝑔_∗(𝜁)𝑑𝜁=𝑎−𝜎²₁/2𝑏a.s. (4.11)
Therefore, by the above arguments, we obtain the following.

Theorem 4.1. Let (𝑥(𝑡),𝑦(𝑡)) be the solution of system (1.3) with the initial value (𝑥(0),𝑦(0))∈𝑅²₊. Then,(i)if 𝑎<𝜎²₁/2, then lim_𝑡→∞𝑥(𝑡)=0,lim_𝑡→∞𝑦(𝑡)=0,a.s., (4.12)(ii)if 𝑎>𝜎²₁/2,𝑒+𝜎²₂/2>𝑘𝛼/𝛽, then the distribution of 𝑥(𝑡) converges weakly to the probability measure with density 𝑓^∗(𝜁)=𝐶₀𝜁^{2(𝑎−𝜎21/2)/𝜎21−1}𝑒^{−2𝑏𝜁/𝜎21}, where 𝐶₀=(2𝑏/𝜎²₁)^{2(𝑎−𝜎21/2)/𝜎21}/Γ(2(𝑎−𝜎²₁/2)/𝜎²₁), and lim_𝑡→∞1𝑡^𝑡₀𝑥(𝑠)𝑑𝑠=𝑎−𝜎²₁/2𝑏,lim_𝑡→∞𝑦(𝑡)=0,a.s. (4.13)

Appendix

For the completeness of the paper, in this section, we list some theories about stationary distribution (see [15]).

Let 𝑋(𝑡) be a homogeneous Markov process in 𝐸_𝑙 (𝐸_𝑙 denotes euclidean 𝑙-space) described by 𝑑𝑋(𝑡)=𝑏(𝑋)𝑑𝑡+^𝑘_𝑟=1𝑔_𝑟(𝑋)𝑑𝐵_𝑟(𝑡). (A.1) The diffusion matrix is 𝐴(𝑥)=(𝑎_𝑖𝑗(𝑥)),𝑎_𝑖𝑗(𝑥)=∑^𝑘_𝑟=1𝑔^𝑖_𝑟(𝑥)𝑔^𝑗_𝑟(𝑥).

Assumption B. There exists a bounded domain 𝑈⊂𝐸_𝑙 with regular boundary Γ, having the following properties. (B.1) In the domain 𝑈 and some neighbourhood thereof, the smallest eigenvalue of the diffusion matrix 𝐴(𝑥) is bounded away from zero. (B.2) If 𝑥∈𝐸_𝑙⧵𝑈, the mean time 𝜏 at which a path issuing from 𝑥 reaches the set 𝑈 is finite, and sup_𝑥∈𝐾𝐸_𝑥𝜏<∞ for every compact subset 𝐾⊂𝐸_𝑙.

Lemma A.1 (see [15]). If (B) holds, then the Markov process 𝑋(𝑡) has a stationary distribution 𝜇(⋅). Let 𝑓(⋅) be a function integrable with respect to the measure 𝜇. Then 𝑃_𝑥{lim_𝑇→∞1/𝑇∫^𝑇₀𝑓(𝑋(𝑡))𝑑𝑡=∫_𝐸𝑙𝑓(𝑥)𝜇(𝑑𝑥)}=1 for all 𝑥∈𝐸_𝑙.

Remark A.2. The proof is given in [15]. Exactly, the existence of stationary distribution with density is referred to Theorem 4.1, Page 119, and Lemma 9.4, Page 138, in [4]. The weak convergence and the ergodicity is obtained in Theorem 5.1, Page 121, and Theorem 7.1, Page 130, in [4].
To validate (B.1), it suffices to prove that 𝐹 is uniformly elliptical in any bounded domain 𝐷, where 𝐹𝑢=𝑏(𝑥)⋅𝑢_𝑥+(1/2)𝑡𝑟(𝐴(𝑥)𝑢_𝑥𝑥); that is, there is a positive number 𝑀 such that ∑^𝑘_𝑖,𝑗=1𝑎_𝑖𝑗(𝑥)𝜉_𝑖𝜉_𝑗≥𝑀|𝜉|²,𝑥∈𝐷,𝜉∈𝑅^𝑘 (see Chapter 3, Page 103 of [19] and Rayleigh's principle in [20, Chapter 6, Page 349]). To verify (B.2), it is sufficient to show that there exists some neighborhood 𝑈 and a nonnegative 𝐶²-function such that and for any 𝐸_𝑙⧵𝑈,𝐿𝑉 is negative (for details refer to [21, Page 1163]).

Lemma A.3. Let 𝑋(𝑡) be a regular temporally homogeneous Markov process in 𝐸_𝑙. If 𝑋(𝑡) is recurrent relative to some bounded domain 𝑈, then it is recurrent relative to any nonempty domain in 𝐸_𝑙.

Acknowledgments

The work was supported by the Ministry of Education of China (no. 109051), the Ph.D. Programs Foundation of Ministry of China (no. 200918), NSFC of China (no. 10971021), and the Graduate Innovative Research Project of NENU (no. 09SSXT117), Youth Fund of Jiangsu Province (no. BK2012208), and the Tian Yuan Special Funds of the National Natural Science Foundation of China (no. 11226205).