Dynamics of a Stochastic Cooperative Predator-Prey System with Beddington-DeAngelis Functional Response

Abstract

Stochastic cooperative predator-prey system with Beddington-DeAngelis functional response is studied. It presents an investigation of dynamic properties of the system. Our results show that there exists a unique positive solution to the system for any positive initial value, and the positive solution is stochastically bounded. Moreover, under some conditions, we analyze global asymptotic stability of the positive solutions. With small environmental noises, the stochastic system is getting more similar to the corresponding deterministic system. Neither of the species in the system will die out. Finally, simulations are carried out to conform to our result.

1. Introduction

As we all know, in mathematical biology, predator-prey system, competitive system, and cooperative system are the three rudimentary and important ecological systems. The dynamic relationship between species has long been and will continue to be a dominant theme in ecology due to its universal existence and importance. It is well-known that predator-prey systems are very important and extensive in the nature fields. One significant component of the predator-prey relationship is predator’s functional response, that is, the rate of prey consumption by an average predator. There are many significant functional responses in order to model various different situations. In fact, most of the functional responses are prey-dependent; however, some biologists have argued that in many cases, especially when predators have to search for food and therefore have to share or compete for food, the traditional predator-prey systems with prey-dependent functional response fail to model the interference among predators, the functional response should be predator-dependent. In [1, 2], Beddington and DeAngelis proposed the following predator-prey model with Beddington-DeAngelis functional response: 𝑑𝑥()−𝑑𝑡=𝑥𝑟−𝜃𝑥𝐸𝑥𝑦,𝑎+𝑏𝑥+𝑐𝑦𝑑𝑦𝑑𝑡=−𝑑𝑦+𝛽𝑥𝑦.𝑎+𝑏𝑥+𝑐𝑦 (1) Skalski and Gilliam [3] compared statistical evidence from 19 predator-prey systems, and then they claimed that three predator-dependent functional responses (Hassell-Varley, Beddington-DeAngelis, and Growley-Martin) can provide better description of predator feeding over a range of predator-prey abundances. And the Beddington-DeAngelis type functional response was even suitable in some cases.

But most of this work is restricted to predator-prey systems, little has been done for cooperative systems [4, 5]. May [6] suggested the following set of equations: 𝑑𝑥𝑑𝑡=𝑟₁𝑥𝑥1−𝑎₁+𝑏₁𝑦−𝑐₁𝑥,𝑑𝑦𝑑𝑡=𝑟₂𝑦𝑦1−𝑎₂+𝑏₂𝑥−𝑐₂𝑦, (2) to describe a pair of mutualists.

However, there is often the interaction among multiple species in nature, whose relationship is more complex than those in two species. Therefore, it is more realistic to consider the multiple-species predator-prey systems. In order to continue studying such models, in this paper, we consider a cooperative predator-prey system with Beddington-DeAngelis functional responses at first: 𝑑𝑥=𝑥𝑎₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧𝑎𝑑𝑡,𝑑𝑦=𝑦₂−𝑏₂ℎ𝑦−²𝑦𝑓₂+𝑔₂𝑧+𝑑₂𝑥1+𝛼₂𝑥+𝛽₂𝑦𝑎𝑑𝑡,𝑑𝑧=𝑧₃−𝑏₃ℎ𝑧−³𝑧𝑓₃+𝑔₃𝑦+𝑑₃𝑥1+𝛼₃𝑥+𝛽₃𝑧𝑑𝑡, (3) where species 𝑥 is the prey of 𝑦 and 𝑧, 𝑦, and 𝑧 are cooperative species. All the parameters in system (3) are positive constants.

In fact, population dynamics is inevitably affected by environmental white noise which is an important component in an ecosystem. But the model (3) is deterministic and does not incorporate the effect of environmental noise. May [6] also pointed out the fact that due to environmental fluctuation, the birth rates, carrying capacity, competition coefficients, and other parameters involved in system exhibit random fluctuation to a greater or a lesser extent. Therefore many scholars rewrote the deterministic models as stochastic ones subjected to stochastic noises, for studying the effect of environmental variability on the population dynamics [7–9].

The parameters in the real ecosystems are often subject to lots of environmental noises, since they relate to climate, geographical distribution, geological features, human disaster, human intervention, and other environmental factors. Therefore, the logistics and energy flow, in which they are determined by groups, are fluctuating. The oscillation in population biomass is directly manifested as birth and death rates of random perturbation. Currently, one of the main ways considered in the literature to model the effect of the environmental fluctuations in population dynamics is to assume that the most sensitive parameter is the intrinsic growth rate. Thus, in this paper we introduce some stochastic perturbation into the intrinsic growth rate. Therefore, the intrinsic growth rate can be written as an average growth rate adding some small random perturbed terms. In general, by the well-known central limit theorem, the small terms follow some normal distributions, so we can use standard Brownian motions to represent the environmental fluctuations.

In this paper, taking into account the effect of randomly fluctuating environment, we introduce stochastic perturbation into growth rates 𝑎₁, 𝑎₂, and 𝑎₃ to become 𝑎₁+𝜎₁𝐵̇₁(𝑡), 𝑎₂+𝜎₂𝐵̇₂(𝑡), and 𝑎₃+𝜎₃𝐵̇₃(𝑡) in system, where 𝜎²_𝑖 represents the intensity of the noise and 𝐵̇_𝑖(𝑡) is a standard white noise, namely, is 𝐵_𝑖(𝑡) a standard Brownian motion defined on a complete probability space (Ω,ℱ,𝒫). Then the stochastic system takes the following form: 𝑑𝑥=𝑥𝑎₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧𝑑𝑡+𝜎₁𝑥𝑑𝐵₁(𝑡),𝑑𝑦=𝑦𝑎₂−𝑏₂ℎ𝑦−²𝑦𝑓₂+𝑔₂𝑧+𝑑₂𝑥1+𝛼₂𝑥+𝛽₂𝑦𝑑𝑡+𝜎₂𝑦𝑑𝐵₂(𝑡),𝑑𝑧=𝑧𝑎₃−𝑏₃ℎ𝑧−³𝑧𝑓₃+𝑔₃𝑦+𝑑₃𝑥1+𝛼₃𝑥+𝛽₃𝑧𝑑𝑡+𝜎₃𝑧𝑑𝐵₃(𝑡). (4) Considering system (4), the initial conditions 𝑥(0)>0, 𝑦(0)>0, and 𝑧(0)>0 will be referred to.

2. Global Positive Solutions

Lemma 1. For any initial value (𝑥₀,𝑦₀,𝑧₀)∈𝑅³₊, where 𝑅³₊={(𝑥,𝑦,𝑧)∣𝑥>0,𝑦>0,𝑧>0}, system (4) has a unique positive local solution (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) for 𝑡∈[0,𝜏_𝑒] almost surely, where 𝜏_𝑒 is the explosion time.

Proof. Consider the following equations: 𝑑𝑢=𝑎₁−12𝜎²₁−𝑏₁𝑒^𝑢−𝑐₂𝑒^𝑣1+𝛼₂𝑒^𝑣+𝛽₂𝑒^𝑢−𝑐₃𝑒^𝑤1+𝛼₃𝑒^𝑣+𝛽₃𝑒^𝑤𝑑𝑡+𝜎₁𝑑𝐵₁(𝑡),𝑑𝑣=𝑎₂−12𝜎²₂−𝑏₂𝑒^𝑣−ℎ₂𝑒^𝑣𝑓₂+𝑔₂𝑒^𝑤+𝑑₂𝑒^𝑢1+𝛼₂𝑒^𝑣+𝛽₂𝑒^𝑢𝑑𝑡+𝜎₂𝑑𝐵₂(𝑡),𝑑𝑤=𝑎₃−12𝜎²₃−𝑏₃𝑒^𝑤−ℎ₃𝑒^𝑣𝑓₃+𝑔₃𝑒^𝑤+𝑑₃𝑒^𝑢1+𝛼₃𝑒^𝑤+𝛽₃𝑒^𝑢𝑑𝑡+𝜎₃𝑑𝐵₃(𝑡), (5) on 𝑡>0 with initial value 𝑢(0)=ln𝑥(0), 𝑣(0)=ln𝑦(0), and 𝑤(0)=ln𝑧(0). The coefficients of (5) satisfy the local Lipschitz condition, thus there is a unique local solution (𝑢(𝑡),𝑣(𝑡),𝑤(𝑡)) on [0,𝜏_𝑒). Then 𝑥(𝑡)=𝑒^𝑢(𝑡), 𝑦(𝑡)=𝑒^𝑣(𝑡), and 𝑧(𝑡)=𝑒^𝑤(𝑡) are the unique positive local solutions with initial value 𝑥₀>0, 𝑦₀>0, and 𝑧₀>0 by Itô’s formula.

Theorem 2. For any initial value (𝑥₀,𝑦₀,𝑧₀)∈𝑅³₊, there is a unique solution (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) of system (4) on 𝑡⩾0, and the solution will remain in 𝑅³₊ with probability 1.

Proof. According to Lemma 1, we only need to show that 𝜏_𝑒=∞. Let 𝑛₀>0 be sufficiently large for 𝑥₀, 𝑦₀, and 𝑧₀ lying within the interval [1/𝑛₀,𝑛₀]. For each integer 𝑛>𝑛₀, we define the stopping times 𝜏_𝑛=inf𝑡∈0,𝜏_𝑒(𝑡)∉1∶𝑥𝑛,𝑦(𝑡)∉1,𝑛𝑛(𝑡)∉1,𝑛or𝑧𝑛,𝑛. (6) Obviously, 𝜏_𝑛 is increasing as 𝑛→∞. Set 𝜏_∞=lim_𝑛→∞𝜏_𝑛; hence, 𝜏_∞⩽𝜏_𝑒 almost surely. Now, we only need to show that 𝜏_∞=∞. If this statement is false, there is a pair of constants 𝑇>0 and 𝜖∈(0,1) such that 𝒫{𝜏_∞⩽𝑇}>𝜖. Thus there exists an integer 𝑛₁>𝑛₀ such that 𝒫𝜏_𝑛⩽𝑇⩾𝜖,𝑛⩾𝑛₁. (7) Define a 𝐶³ function 𝑉: 𝑅³₊→𝑅₊ by 𝑉𝑥,𝑦,𝑧=(𝑥−1−ln𝑥)+𝑦−1−ln𝑦+(𝑧−1−ln𝑧). (8) The nonnegativity of this function can be seen from 𝑢−1−ln𝑢⩾0 and ∀𝑢>0. If (𝑥₀,𝑦₀,𝑧₀)∈𝑅³₊, we have 𝑑𝑉=(𝑥−1)𝑎₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧𝜎𝑑𝑡+²¹2+𝑎𝑑𝑡𝑦−1₂−𝑏₂ℎ𝑦−²𝑦𝑓₂+𝑔₂𝑧+𝑑₂𝑥1+𝛼₂𝑥+𝛽₂𝑦𝜎𝑑𝑡+²²2+()𝑎𝑑𝑡𝑧−1₃−𝑏₃ℎ𝑧−³𝑧𝑓₃+𝑔₃𝑦+𝑑₃𝑥1+𝛼₃𝑥+𝛽₃𝑧𝜎𝑑𝑡+²³2+()𝜎𝑑𝑡𝑥−1₁𝑑𝐵₁(𝑡)+𝑦−1𝜎₂𝑑𝐵₂(𝑡)+(𝑧−1)𝜎₃𝑑𝐵₃(𝑡), (9) then 𝐿𝑉=(𝑥−1)𝑎₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧+𝜎²₁2+𝑦−1𝑎₂−𝑏₂ℎ𝑦−²𝑦𝑓₂+𝑔₂𝑧+𝑑₂𝑥1+𝛼₂𝑥+𝛽₂𝑦+𝜎²₂2+(𝑧−1)𝑎₃−𝑏₃ℎ𝑧−³𝑧𝑓₃+𝑔₃𝑦+𝑑₃𝑥1+𝛼₃𝑥+𝛽₃𝑧+𝜎²₃2⩽(𝑥−1)𝑎₁−𝑏₁𝑥+𝑦−1𝑎₂−𝑏₂𝑑𝑦+²𝛼₂+(𝑧−1)𝑎₃−𝑏₃𝑑𝑧+³𝛼₃+𝜎²₁+𝜎²₂+𝜎²₃2⩽𝑥𝑎₁−𝑏₁𝑥+𝑦𝑎₂−𝑏₂𝑑𝑦+²𝛼₂+𝑧𝑎₃−𝑏₃𝑑𝑧+³𝛼₃+𝜎²₁+𝜎²₂+𝜎²₃2⩽𝐾, (10) where 𝐾 is a positive number. Substituting this inequality into (9), we see that 𝑑𝑉𝑥,𝑦,𝑧⩽𝐾𝑑𝑡+(𝑥−1)𝜎₁𝑑𝐵₁(𝑡)+𝑦−1𝜎₂𝑑𝐵₂(𝑡)+(𝑧−1)𝜎₃𝑑𝐵₃(𝑡). (11) Integrating both sides of the above inequality from 0 to 𝜏_𝑘∧𝑇 and then taking the expectations leads to 𝐸𝑉𝑥𝜏_𝑛∧𝑇,𝑦𝜏_𝑛∧𝑇,𝑧𝜏_𝑛𝑥∧𝑇⩽𝑉₀,𝑦₀,𝑧₀+𝐾𝑇. (12) Set Ω_𝑛={𝜏_𝑛⩽𝑇}, then we get 𝒫(Ω_𝑛)⩾𝜖 by inequality (7). Obviously, for every 𝜔∈Ω_𝑛, there are at least 𝑥(𝜏_𝑛,𝜔), 𝑦(𝜏_𝑛,𝜔), and 𝑧(𝜏_𝑛,𝜔) which equal either 𝑘 or 1/𝑘, then 𝑉(𝑥(𝜏_𝑛,𝜔),𝑦(𝜏_𝑛,𝜔),𝑧(𝜏_𝑛,𝜔)) is no less than ),1min(𝑛−1−ln𝑛𝑛1−1−ln𝑛. (13) It then follows from (9) that 𝑉𝑥₀,𝑦₀,𝑧₀1+𝐾𝑇⩾𝐸_Ω𝑛(𝜔)𝑉𝑥𝜏_𝑛,𝜔,𝑦𝜏_𝑛,𝜔,𝑧𝜏_𝑛),1,𝜔⩾𝜖min(𝑛−1−ln𝑛𝑛1−1−ln𝑛, (14) where 1_Ω𝑛(𝜔) is the indicator function of Ω_𝑛, letting 𝑛→∞, we have that ∞>𝑉𝑥₀,𝑦₀,𝑧₀+𝐾𝑇=∞ (15) This completes the proof.

3. Stochastic Boundedness

Definition 3. The solution (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) of system (4) is said to be stochastically ultimately bounded, if for any 𝜖∈(0,1), there is a positive constant 𝛿=𝛿(𝜖), such that for any initial value (𝑥₀,𝑦₀,𝑧₀)∈𝑅³₊, the solution (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) of system (4) has the property that limsup_𝑛→∞𝒫𝑥(𝑡)(𝑡)(𝑡=,𝑦,𝑧)𝑥²(𝑡)+𝑦²(𝑡)+𝑧²(𝑡)>𝛿<𝜖. (16)

Assumption A. For any initial value (𝑥₀,𝑦₀,𝑧₀)∈𝑅³₊, there exists 𝑝⩾1 such that 𝑥(0)<𝑎₁+(1/2)𝑝−1𝜎²₁𝑏₁,𝑦(0)<𝑎₂+𝑑₂/𝛼₂+(1/2)𝑝−1𝜎²₂𝑏₂,𝑧(0)<𝑎₃+𝑑₃/𝛼₃+(1/2)𝑝−1𝜎²₃𝑏₃. (17)

Lemma 4. Assume that Assumption A holds. Let (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) be a positive solution of (4) with any initial value (𝑥₀,𝑦₀,𝑧₀)∈𝑅³₊, for all 𝑝>1, then 𝐸𝑥^𝑝(𝑡)⩽𝐾₁𝑝,𝐸𝑦^𝑝(𝑡)⩽𝐾₂𝑝,𝐸𝑧^𝑝(𝑡)⩽𝐾₃𝑝, (18) where 𝐾₁𝑝𝑎∶=¹+(1/2)𝑝−1𝜎²₁𝑏₁^𝑝,𝐾₂𝑝𝑎∶=²+𝑑₂/𝛼₂+(1/2)𝑝−1𝜎²₂𝑏₂^𝑝,𝐾₃𝑝𝑎∶=³+𝑑₃/𝛼₃+(1/2)𝑝−1𝜎²₃𝑏₃^𝑝. (19)

Proof. Define the function 𝑉₁=𝑥^𝑝, for (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡))∈𝑅³₊ and 𝑝>0. By Itô’s formula we get 𝑑𝑉₁=𝑝𝑥₁^𝑝−11𝑑𝑥+2𝑝𝑝−1𝑥₁^𝑝−2(𝑑𝑥)²=𝑝𝑥^𝑝𝑎₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧𝑑𝑡+𝜎₁𝑑𝐵₁(𝑡+1)2𝑝𝑝−1𝜎²₁𝑑𝑡=𝑝𝑥^𝑝𝑎₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧+12𝑝𝑝−1𝜎²₁𝑑𝑡+𝜎₁𝑑𝐵₁(𝑡). (20) Integrating from 0 to 𝑡 and taking expectations yields 𝐸𝑥^𝑝(𝑡)−𝐸𝑥^𝑝(0=_𝑡0)𝑝𝐸𝑥^𝑝𝑎₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧+12𝑝𝑝−1𝜎²₁𝑑𝑠. (21) So, 𝑑𝐸𝑥^𝑝(𝑡)𝑥𝑑𝑡=𝑝𝐸^𝑝𝑎₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧+12𝑝𝑝−1𝜎²₁⩽𝑝𝑎₁+12𝑝−1𝜎²₁𝐸𝑥^𝑝(𝑡)−𝑏₁𝐸𝑥^𝑝+1(𝑡𝑎)⩽𝑝₁+12𝑝−1𝜎²₁𝐸𝑥^𝑝(𝑡)−𝑏₁𝐸𝑥^𝑝(𝑡)^{(𝑝+1)/𝑝}=𝑝𝐸𝑥^𝑝(𝑡)𝑎₁+12𝑝−1𝜎²₁−𝑏₁𝐸𝑥^𝑝(𝑡)^1/𝑝. (22) Let 𝑋(𝑡)=𝐸[𝑥^𝑝(𝑡)], then we have 𝑑𝑋(𝑡)(𝑡𝑎𝑑𝑡⩽𝑝𝑋)₁+12𝑝−1𝜎²₁−𝑏₁𝑋^1/𝑝(𝑡). (23) From (17), we know 0<𝑏₁𝑋^1/𝑝(0)=𝑏₁𝑥(0)<𝑎₁+12𝑝−1𝜎²₁, (24) which by the standard comparison argument shows that 𝐸𝑥^𝑝(𝑡)^1/𝑝=𝑋^1/𝑝⩽𝑎₁+(1/2)𝑝−1𝜎²₁𝑏₁, (25) that is, 𝐸𝑥^𝑝(𝑡⩽𝑎)¹+(1/2)𝑝−1𝜎²₁𝑏₁^𝑝. (26)
Define the function 𝑉₂=𝑦^𝑝, for (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡))∈𝑅³₊ and 𝑝>0. By Itô’s formula we get 𝑑𝑉₂=𝑝𝑦^𝑝−11𝑑𝑦+2𝑝𝑝−1𝑦^𝑝−2𝑑𝑦²=𝑝𝑦^𝑝𝑎₂−𝑏₂ℎ𝑦−²𝑦𝑓₂+𝑔₂𝑧+𝑑₂𝑥1+𝛼₂𝑥+𝛽₂𝑦𝑑𝑡+𝜎₂𝑑𝐵₂(𝑡)+12𝑝−1𝜎²₂𝑑𝑡=𝑝𝑦^𝑝𝑎₂−𝑏₂ℎ𝑦−²𝑦𝑓₂+𝑔₂𝑧+𝑑₂𝑥1+𝛼₂𝑥+𝛽₂𝑦+12𝑝−1𝜎²₂𝑑𝑡+𝜎₂𝑑𝐵₂(𝑡). (27) Integrating from 0 to 𝑡 and taking expectations yields 𝐸𝑦^𝑝(𝑡)−𝐸𝑦^𝑝(0=_𝑡0)𝑝𝐸𝑦^𝑝𝑎₂−𝑏₂ℎ𝑦−²𝑦𝑓₂+𝑔₂𝑧+𝑑₂𝑥1+𝛼₂𝑥+𝛽₂𝑦+12𝑝−1𝜎²₂𝑑𝑠. (28) So, 𝑑𝐸𝑦^𝑝(𝑡)𝑦𝑑𝑡=𝑝𝐸^𝑝𝑎₂−𝑏₂ℎ𝑦−²𝑦𝑓₂+𝑔₂𝑧+𝑑₂𝑥1+𝛼₂𝑥+𝛽₂𝑦+12𝑝−1𝜎²₂⩽𝑝𝑎₂+𝑑₂𝛼₂+12𝑝−1𝜎²₂𝐸𝑦^𝑝(𝑡)−𝑏₂𝐸𝑦^𝑝+1(𝑡𝑎)⩽𝑝₂+𝑑₂𝛼₂+12𝑝−1𝜎²₂𝐸𝑦^𝑝(𝑡)−𝑏₂𝐸𝑦^𝑝(𝑡)^{(𝑝+1)/𝑝}=𝑝𝐸𝑦^𝑝(𝑡)𝑎₂+𝑑₂𝛼₂+12𝑝−1𝜎²₂−𝑏₂𝐸𝑦^𝑝(𝑡)^1/𝑝. (29) Let 𝑌(𝑡)=𝐸[𝑦^𝑝(𝑡)], then we have 𝑑𝑌(𝑡)(𝑡𝑎𝑑𝑡⩽𝑝𝑌)₂+𝑑₂𝛼₂+12𝑝−1𝜎²₂−𝑏₂𝑌^1/𝑝(𝑡). (30) From (17), we know that 0<𝑏₂𝑌^1/𝑝(0)=𝑏₂𝑦(0)<𝑎₂+𝑑₂𝛼₂+12𝑝−1𝜎²₂, (31) which by the standard comparison argument shows that 𝐸𝑦^𝑝(𝑡)^1/𝑝=𝑌^1/𝑝(𝑡)⩽𝑎₂+𝑑₂/𝛼₂+(1/2)𝑝−1𝜎²₂𝑏₂, (32) that is, 𝐸𝑦^𝑝(𝑡⩽𝑎)²+𝑑₂/𝛼₂+(1/2)𝑝−1𝜎²₂𝑏₂^𝑝. (33) Similarly, we can show that 𝐸𝑧^𝑝(𝑡⩽𝑎)³+𝑑₃/𝛼₃+(1/2)𝑝−1𝜎²₃𝑏₃^𝑝. (34) This completes the proof.

Theorem 5. Assume that Assumption A holds, the solutions of system (4) with initial value (𝑥₀,𝑦₀,𝑧₀)∈𝑅³₊ are stochastically ultimately bounded.

Proof. If (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡))∈𝑅³, its norm here is denoted by |𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)|=(𝑥²(𝑡)+𝑦²(𝑡)+𝑧²(𝑡))^1/2, then 𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)^𝑝⩽3^𝑝/2|𝑥(𝑡)|^𝑝+𝑦(𝑡)^𝑝+|𝑧(𝑡)|^𝑝 (35) by Lemma 4, 𝐸[|𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)|^𝑝]⩽𝐾(𝑝), 𝑡∈(0,+∞). 𝐾(𝑝) is dependent on (𝑥₀,𝑦₀,𝑧₀) and is defined by 𝐾(𝑝)=3^𝑝/2(𝐾₁(𝑝)+𝐾₂(𝑝)+𝐾₃(𝑝)). By virtue of Chebyshev inequality, we can easily obtain that the solution (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) of system (4) is stochastically ultimately bounded.

4. Stochastic Permanence

Definition 6 (see [10]). The solution (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) of system (4) is said to be stochastically permanent, if for any 𝜖∈(0,1), there exist a pair of positive constants 𝛿=𝛿(𝜖) and 𝜒=𝜒(𝜖) such that for any initial value (𝑥₀,𝑦₀,𝑧₀)∈𝑅³₊, the solution (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) of system (4) has the properties that liminf_𝑡⇀∞𝒫𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)⩾𝛿⩾1−𝜖,liminf_𝑡⇀∞𝒫𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)⩽𝜒⩾1−𝜖. (36)

Assumption B. One has
1/2(max{𝜎₁,𝜎₂,𝜎₃})²<min{𝑎₁,𝑎₂−(𝑐₂/𝛼₂),𝑎₃−(𝑐₃/𝛼₃)}.

Theorem 7. Under Assumption B, for any initial value (𝑥₀,𝑦₀,𝑧₀)∈𝑅³₊, the solution (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) of system (4) satisfies that limsup_𝑡→∞𝐸1𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)^𝜃⩽𝐻, (37) where 𝜃 is an arbitrary positive constant satisfying 𝜃+12max𝜎₁,𝜎₂,𝜎₃²<min𝑎₁,𝑎₂−𝑐₂𝛼₂,𝑎₃−𝑐₃𝛼₃, (38) and 𝑘 is an arbitrary positive constant satisfying 𝜃min𝑎₁,𝑎₂−𝑐₂𝛼₂,𝑎₃−𝑐₃𝛼₃−𝜃(𝜃+1)2max𝜎₁,𝜎₂,𝜎₃²−𝑘>0. (39)

Proof. Define 𝑉(𝑥,𝑦,𝑧)=𝑥+𝑦+𝑧 for (𝑥,𝑦,𝑧)∈𝑅³₊, then 𝑎𝑑𝑉𝑥,𝑦,𝑧=𝑥₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧𝑎𝑑𝑡+𝑦₂−𝑏₂ℎ𝑦−²𝑦𝑓₂+𝑔₂𝑧+𝑑₂𝑥1+𝛼₂𝑥+𝛽₂𝑦𝑎𝑑𝑡+𝑧₃−𝑏₃ℎ𝑧−³𝑧𝑓₃+𝑔₃𝑦+𝑑₃𝑥1+𝛼₃𝑥+𝛽₃𝑧𝑑𝑡+𝜎₁𝑥𝑑𝐵₁(𝑡)+𝜎₂𝑦𝑑𝐵₂(𝑡)+𝜎₃𝑧𝑑𝐵₃(𝑡). (40) Define 𝑈(𝑥,𝑦,𝑧)=(1/(𝑉(𝑥,𝑦,𝑧))) for (𝑥,𝑦,𝑧)∈𝑅³₊, by Itô’s formula, we get 𝑑𝑈=−𝑈²𝑥𝑎₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧+𝑦𝑎₂−𝑏₂ℎ𝑦−²𝑦𝑓₂+𝑔₂𝑧+𝑑₂𝑥1+𝛼₂𝑥+𝛽₂𝑦+𝑧𝑎₃−𝑏₃ℎ𝑧−³𝑧𝑓₃+𝑔₃𝑦+𝑑₃𝑥1+𝛼₃𝑥+𝛽₃𝑧𝑑𝑡+𝑈³𝜎²₁𝑥²+𝜎²₂𝑦²+𝜎²₃𝑧²𝑑𝑡−𝑈²𝜎₁𝑥𝑑𝐵₁(𝑡)+𝜎₂𝑦𝑑𝐵₂(𝑡)+𝜎₃𝑧𝑑𝐵₃(𝑡)=𝐿𝑈𝑑𝑡−𝑈²𝜎₁𝑥𝑑𝐵₁(𝑡)+𝜎₂𝑦𝑑𝐵₂(𝑡)+𝜎₃𝑧𝑑𝐵₃(𝑡), (41) where 𝐿𝑈=−𝑈²𝑥𝑎₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧+𝑦𝑎₂−𝑏₂ℎ𝑦−²𝑦𝑓₂+𝑔₂𝑧+𝑑₂𝑥1+𝛼₂𝑥+𝛽₂𝑦+𝑧𝑎₃−𝑏₃ℎ𝑧−³𝑧𝑓₃+𝑔₃𝑦+𝑑₃𝑥1+𝛼₃𝑥+𝛽₃𝑧+𝑈³𝜎²₁𝑥²+𝜎²₂𝑦²+𝜎²₃𝑧². (42)
Under Assumption B, choosing a positive constant 𝜃 such that it satisfies (38). By Itô’s formula, we get 𝑑(1+𝑈)^𝜃=𝜃(1+𝑈)^𝜃−1𝜃()𝐿𝑈+𝜃−12𝑈⁴(1+𝑈)^𝜃−2×𝜎²₁𝑥²+𝜎²₂𝑦²+𝜎²₃𝑧²𝑑𝑡−𝜃𝑈²(1+𝑈)^𝜃−1𝜎₁𝑥𝑑𝐵₁(𝑡)+𝜎₂𝑦𝑑𝐵₂(𝑡)+𝜎₃𝑧𝑑𝐵₃(𝑡), (43) where 𝐿(1+𝑈)^𝜃=𝜃(1+𝑈)^𝜃−1𝜃()𝐿𝑈+𝜃−12𝑈⁴×(1+𝑈)^𝜃−2𝜎²₁𝑥²+𝜎²₂𝑦²+𝜎²₃𝑧², (44) then choosing a positive constant 𝑘 such that it obeys (39), by Itô’s formula, 𝑑𝑒^𝑘𝑡(1+𝑈)^𝜃=𝑒^𝑘𝑡𝑑(1+𝑈)^𝜃+𝑘𝑒^𝑘𝑡(1+𝑈)^𝜃𝑑𝑡, (45) where 𝐿𝑒^𝑘𝑡(1+𝑈)^𝜃=𝑒^𝑘𝑡𝐿(1+𝑈)^𝜃+𝑘𝑒^𝑘𝑡(1+𝑈)^𝜃=𝑒^𝑘𝑡(1+𝑈)^𝜃−2×𝑘(1+𝑈)²+𝜃(𝜃+1)2𝑈⁴𝜎²₁𝑥²+𝜎²₂𝑦²+𝜎²₃𝑧²+𝜃𝑈³𝜎²₁𝑥²+𝜎²₂𝑦²+𝜎²₃𝑧²−𝜃(1+𝑈)𝑈²×𝑥𝑎₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧+𝑦𝑎₂−𝑏₂ℎ𝑦−²𝑦𝑓₂+𝑔₂𝑧+𝑑₂𝑥1+𝛼₂𝑥+𝛽₂𝑦+𝑧𝑎₃−𝑏₃ℎ𝑧−³𝑧𝑓₃+𝑔₃𝑦+𝑑₃𝑥1+𝛼₃𝑥+𝛽₃𝑧⩽𝑒^𝑘𝑡(1+𝑈)^𝜃−2×𝑘(1+𝑈)²+𝜃(𝜃+1)2𝑈⁴𝜎²₁𝑥²+𝜎²₂𝑦²+𝜎²₃𝑧²+𝜃𝑈³𝜎²₁𝑥²+𝜎²₂𝑦²+𝜎²₃𝑧²−𝜃𝑈²𝑎₁𝑥+𝑎₂−𝑐₂𝛼₂𝑦+𝑎₃−𝑐₃𝛼₃𝑧−𝑏₁𝑥²+𝑏₂+ℎ₂𝑓₂𝑦+𝑏₃+ℎ₃𝑓₃𝑧−𝜃𝑈³𝑎₁𝑥+𝑎₂−𝑐₂𝛼₂𝑦+𝑎₃−𝑐₃𝛼₃𝑧−𝑏₁𝑥²+𝑏₂+ℎ₂𝑓₂𝑦+𝑏₃+ℎ₃𝑓₃𝑧⩽𝑒^𝑘𝑡(1+𝑈)^𝜃−2×𝑘+𝜃max𝑏₁,𝑏₂+ℎ₂𝑓₂,𝑏₃+ℎ₃𝑓₃+2𝑘−𝜃min𝑎₁,𝑎₂−𝑐₂𝛼₂,𝑎₃−𝑐₃𝛼₃+𝜃max𝑏₁,𝑏₂+ℎ₂𝑓₂,𝑏₃+ℎ₃𝑓₃+𝜃max𝜎₁,𝜎₂,𝜎₃²𝑈+𝑘−𝜃min𝑎₁,𝑎₂−𝑐₂𝛼₂,𝑎₃−𝑐₃𝛼₃+𝜃(𝜃+1)2max𝜎₁,𝜎₂,𝜎₃²𝑈², (46) where 𝑈³𝜎²₁𝑥²+𝜎²₂𝑦²+𝜎²₃𝑧²⩽max𝜎₁,𝜎₂,𝜎₃²𝜃()𝑈,𝜃+12𝑈⁴𝜎²₁𝑥²+𝜎²₂𝑦²+𝜎²₃𝑧²⩽𝜃(𝜃+1)2max𝜎₁,𝜎₂,𝜎₃²𝑈². (47) Hence, it implies that there exists a positive constant 𝐾 such that 𝐿𝑒^𝑘𝑡(1+𝑈)^𝜃⩽𝐾𝑒^𝑘𝑡. (48) Then we have 𝑑𝑒^𝑘𝑡(1+𝑈)^𝜃⩽𝐾𝑒^𝑘𝑡𝑑𝑡−𝑒^𝑘𝑡𝜃𝑈²(1+𝑈)^𝜃−1×𝜎₁𝑥𝑑𝐵₁(𝑡)+𝜎₂𝑦𝑑𝐵₂(𝑡)+𝜎₃𝑧𝑑𝐵₃(𝑡). (49) Integrating both sides of the above inequality from 0 to 𝑡 and then taking the expectations leads to 𝐸𝑒^𝑘𝑡(1+𝑈(𝑡))^𝜃⩽(1+𝑈(0))^𝜃+𝐾𝑘𝑒^𝑘𝑡=(1+𝑈(0))^𝜃+𝐻₁𝑒^𝑘𝑡, (50) where 𝐻₁=𝐾/𝑘. So limsup_𝑡→∞𝐸𝑈(𝑡)^𝜃⩽limsup_𝑡→∞𝐸(1+𝑈(𝑡))^𝜃⩽𝐻₁. (51) Since that (𝑥+𝑦+𝑧)^𝜃⩽3^𝜃(𝑥²+𝑦²+𝑧²)^𝜃/2=3^𝜃|𝑥,𝑦,𝑧|^𝜃, where (𝑥,𝑦,𝑧)∈𝑅³₊, obviously limsup_𝑡→∞𝐸1𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)^𝜃⩽3^𝜃limsup_𝑡→∞𝐸𝑈(𝑡)^𝜃⩽3^𝜃𝐻₁=∶𝐻, (52) as required.

Theorem 8. Under Assumption B, system (4) is stochastically permanent.

Proof. By Theorem 5, we know that limsup_𝑡→∞𝐸𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)^𝑝⩽𝐾𝑝. (53) Now, for any 𝜖>0, let 𝜒=(𝐾(𝑝)/𝜖)^1/𝑝. Then by Chebyshev’s inequality, we can obtain the conclusion easily.

5. Global Asymptotic Stability

Definition 9. Let (𝑥₁(𝑡),𝑦₁(𝑡),𝑧₁(𝑡)) be a positive solution of system (4). If we say that (𝑥₁(𝑡),𝑦₁(𝑡),𝑧₁(𝑡)) is globally asymptotically stable in expectation, it means that any other solution (𝑥₂(𝑡),𝑦₂(𝑡),𝑧₂(𝑡)) of system (4) has 𝑡⩾0 and that we have initial value (𝑥₀,𝑦₀,𝑧₀)∈𝑅³₊. That is 𝒫lim_𝑡⇀∞𝐸𝑥₁(𝑡),𝑦₁(𝑡),𝑧₁(𝑡)−𝑥₂(𝑡),𝑦₂(𝑡),𝑧₂(𝑡)=0=1. (54)

Lemma 10 (see [11]). Suppose that an n-dimensional stochastic process 𝑋(𝑡) on 𝑡⩾0 satisfies the condition 𝐸|𝑋(𝑡)−𝑋(𝑠)|^𝛼⩾𝑐|𝑡−𝑠|^1+𝛽,0⩽𝑠,𝑡<∞, (55) for some positive constants 𝛼, 𝛽, and 𝑐. Then there exists a continuous modification 𝑋(𝑡) of 𝑋(𝑡) which has the property that for every 𝜗∈(0,𝛽/𝛼) there is a positive random variable ℎ(𝜔) such that 𝒫𝜔∶sup_{||(𝜔)0<𝑡−𝑠<ℎ,0⩽𝑠,𝑡<∞}𝑋(𝑡,𝜔)−𝑋(𝑠,𝜔)|𝑡−𝑠|^𝜗⩽21−2^−𝜗=1. (56) In other words, almost every sample path of 𝑋(𝑡) is locally but uniformly Hölder continuous with exponent 𝜗.

Lemma 11. Let (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) be a solution of system (4) on 𝑡⩾0 with initial value (𝑥₀,𝑦₀,𝑧₀)∈𝑅³₊, then almost every sample path of (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) is uniformly continuous on 𝑡⩾0.

Proof. From system (4), we have the following stochastic integral equation 𝑥(𝑡)(0)+=𝑥_𝑡0𝑥(𝑠)×𝑎₁−𝑏₁𝑥(𝑠)−𝑐₂𝑦(𝑠)1+𝛼₂𝑥(𝑠)+𝛽₂𝑦(𝑠)−𝑐₃𝑧(𝑠)1+𝛼₃𝑥(𝑠)+𝛽₃𝑧(𝑠)+𝑑𝑠_𝑡0𝜎₁𝑥(𝑠)𝑑𝐵₁(𝑠). (57) Let  𝑓(𝑠)=𝑥(𝑠)(𝑎₁−𝑏₁𝑥(𝑠)−(𝑐₂𝑦(𝑠)/(1+𝛼₂𝑥(𝑠)+𝛽₂𝑦(𝑠))) − (𝑐₃𝑧(𝑠)/(1+𝛼₃𝑥(𝑠)+𝛽₃𝑧(𝑠)))), 𝑔(𝑠) = 𝜎₁𝑥(𝑠), notice that 𝐸𝑓(𝑡)^𝑝=𝐸𝑥𝑎₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧^𝑝⩽12𝐸|𝑥|^2𝑝+12𝐸×𝑎₁−𝑏₁𝑐𝑥−²𝑦1+𝛼₂𝑥+𝛽₂𝑦−𝑐₃𝑧1+𝛼₃𝑥+𝛽₃𝑧^2𝑝⩽12𝐸|𝑥|^2𝑝+12𝐸𝑎₁+𝑏₁𝑥+𝑐₂𝑦+𝑐₃𝑧^2𝑝⩽12𝐸|𝑥|^2𝑝+4^2𝑝−12×𝑎₁^2𝑝+𝑏₁^2𝑝𝐸|𝑥|^2𝑝+𝑐₂𝐸𝑦^2𝑝+𝑐₃𝐸|𝑧|^2𝑝⩽12𝐾₁+42𝑝^2𝑝−12×𝑎₁^2𝑝+𝑏₁^2𝑝𝐾₁2𝑝+𝑐₂𝐾₂2𝑝+𝑐₃𝐾₃2𝑝=∶𝐾₄𝑝. (58) On the other hand, by the moment inequality (see [12]) for stochastic integrals, we have that for 0⩽𝑡₁<𝑡₂<∞ and 𝑝>2, 𝐸^𝑡2_𝑡1𝑔(𝑠)𝑑𝐵₁(𝑠)^𝑝⩽𝑝𝑝−12^𝑝/2𝑡₂−𝑡₁^{(𝑝−2)/2𝑡2}_𝑡1𝐸𝑔(𝑠)^𝑝⩽𝑝𝑑𝑠𝑝−12^𝑝/2𝑡₂−𝑡₁^{(𝑝−2)/2𝑡2}_𝑡1𝜎^𝑝₁𝑘₁𝑝=𝑝𝑑𝑠𝑝−12^𝑝/2𝑡₂−𝑡₁^𝑝/2𝐺₁𝑝, (59) where 𝐺₁(𝑝)=𝜎^𝑝₁𝑘₁(𝑝). Let 0⩽𝑡₁<𝑡₂<∞, 𝑡₂−𝑡₁⩽1, and (1/𝑝)+(1/𝑞)=1, we obtain 𝐸𝑥𝑡₂−𝑥𝑡₁^𝑝=𝐸^𝑡2_𝑡1𝑓(𝑠)𝑑𝑠+^𝑡2_𝑡1𝑔(𝑠)𝑑𝐵₁(𝑠)^𝑝⩽2^𝑝−1𝐸^𝑡2_𝑡1𝑓(𝑠)𝑑𝑠^𝑝+2^𝑝−1𝐸^𝑡2_𝑡1𝑔(𝑠)𝑑𝐵₁(𝑠)^𝑝⩽2^{𝑝−1𝑡2}_𝑡11^𝑞𝑑𝑠^𝑝/𝑞𝐸^𝑡2_𝑡1𝑓(𝑠)^𝑝𝑑𝑠+2^𝑝−1𝑝𝑝−12^𝑝/2𝑡₂−𝑡₁^𝑝/2𝐺₁𝑝⩽2^𝑝−1𝑡₂−𝑡₁^𝑝𝐾₄𝑝+2^𝑝−1𝑝𝑝−12^𝑝/2𝑡₂−𝑡₁^𝑝/2𝐺₁𝑝⩽2^𝑝−1𝑡₂−𝑡₁^𝑝/2𝑡₂−𝑡₁^𝑝/2+𝑝𝑝−12^𝑝/2𝐾₅𝑝⩽2^𝑝−1𝑝1+𝑝−12^𝑝/2𝐾₅𝑝𝑡₂−𝑡₁^𝑝/2, (60) where 𝐾₅(𝑝)∶=max{𝐾₄(𝑝),𝐺₁(𝑝)}. Then, we have that almost every sample path of 𝑥(𝑡) is locally but uniformly Hölder continuous with exponent 𝜗 for every 𝜗∈(0,(𝑝−2)/2𝑝) and therefore almost every sample path of 𝑥(𝑡) is uniformly continuous on 𝑡⩾0. Similarly, we can show that almost every sample path of 𝑦(𝑡) and 𝑧(𝑡) is uniformly continuous on 𝑡⩾0.

Lemma 12 (see [13]). Let 𝑓(𝑡) be a nonnegative function defined on [0,∞) such that 𝑓(𝑡) is integrable on [0,∞) and is uniformly continuous on 𝑡⩾0. Then lim_𝑡→∞𝑓(𝑡)=0.

Theorem 13. If 𝐴∶=𝑏₁−𝑐₂𝛼₂𝛽₂−𝑐₃𝛼₃𝛽₃−𝑑₂−𝑑₃>0,𝐵∶=𝑏₂−𝑐₂−2ℎ₂𝑓₂−𝑑₂𝛽₂𝛼₂−ℎ₃𝑔₃𝑎₂+𝑑₂/𝛼₂+𝜎²₂𝑏₂𝑓²₃>0,𝐶∶=𝑏₃−𝑐₃−2ℎ₃𝑓₃−𝑑₃𝛽₃𝛼₃−ℎ₂𝑔₂𝑎₃+𝑑₃/𝛼₃+𝜎²₃𝑏₃𝑓²₂>0, (61) then system (4) is globally asymptotically stable.

Proof. Define 𝑉(𝑡)=ln𝑥₁(𝑡)−ln𝑥₂(𝑡)+ln𝑦₁(𝑡)−ln𝑦₂(𝑡+)ln𝑧₁(𝑡)−ln𝑧₂(𝑡), (62) then 𝑉(𝑡) is a continuous positive function on 𝑡⩾0. A direct calculation of the right differential 𝑑⁺𝑉(𝑡) of 𝑉(𝑡), and then applying Itô’s formula, we have 𝑑⁺𝑉(𝑡)=sgn𝑥₁(𝑡)−𝑥₂(𝑡)×𝑑𝑥₁(𝑡)𝑥₁(𝑡)−𝑑𝑥₁(𝑡)²2𝑥²₁(𝑡)−𝑑𝑥₂(𝑡)𝑥₂(𝑡)−𝑑𝑥₂(𝑡)²2𝑥²₂(𝑡)+sgn𝑦₁(𝑡)−𝑦₂(𝑡)×𝑑𝑦₁(𝑡)𝑦₁(𝑡)−𝑑𝑦₁(𝑡)²2𝑦²₁(𝑡)−𝑑𝑦₂(𝑡)𝑦₂(𝑡)−𝑑𝑦₂(𝑡)²2𝑦²₂(𝑡)+sgn𝑧₁(𝑡)−𝑧₂(𝑡)×𝑑𝑧₁(𝑡)𝑧₁(𝑡)−𝑑𝑧₁(𝑡)²2𝑧²₁(𝑡)−𝑑𝑧₂(𝑡)𝑧₂(𝑡)−𝑑𝑧₂(𝑡)²2𝑧²₂(𝑡)=sgn𝑥₁(𝑡)−𝑥₂(𝑡)×−𝑏₁𝑥₁(𝑡)−𝑥₂(𝑡−𝑐₂)𝑦₁(𝑡)1+𝛼₂𝑥₁(𝑡)+𝛽₂𝑦₁(𝑡)−𝑐₂𝑦₂(𝑡)1+𝛼₂𝑥₂(𝑡)+𝛽₂𝑦₂(𝑡)−𝑐₃𝑧₁(𝑡)1+𝛼₃𝑥₁(𝑡)+𝛽₃𝑧₁(𝑡)−𝑐₃𝑧₂(𝑡)1+𝛼₃𝑥₂(𝑡)+𝛽₃𝑧₂(𝑡)𝑦𝑑𝑡+sgn₁(𝑡)−𝑦₂(𝑡)×−𝑏₂𝑦₁(𝑡)−𝑦₂(𝑡)−ℎ₂𝑦₁(𝑡)𝑓₂+𝑔₂𝑧₁(𝑡)−ℎ₂𝑦₂(𝑡)𝑓₂+𝑔₂𝑧₂(𝑡)+𝑑₂𝑥₁(𝑡)1+𝛼₂𝑥₁(𝑡)+𝛽₂𝑦₁(𝑡)−𝑑₂𝑥₂(𝑡)1+𝛼₂𝑥₂(𝑡)+𝛽₂𝑦₂(𝑡)𝑧𝑑𝑡+sgn₁(𝑡)−𝑧₂(𝑡)×−𝑏₃𝑧₁(𝑡)−𝑧₂(𝑡)−ℎ₃𝑧₁(𝑡)𝑓₃+𝑔₃𝑦₁(𝑡)−ℎ₃𝑧₂(𝑡)𝑓₃+𝑔₃𝑦₂(𝑡)+𝑑₃𝑥₁(𝑡)1+𝛼₃𝑥₁(𝑡)+𝛽₃𝑧₁(𝑡)−𝑑₃𝑥₂(𝑡)1+𝛼₂𝑥₂(𝑡)+𝛽₂𝑧₂(𝑡)𝑑𝑡. (63) Integrating from 0 to 𝑡 and taking expectations yields 𝐸[𝑉(𝑡)]−𝐸[𝑉(0)]=𝐸^𝑡₀sgn𝑥₁(𝑠)−𝑥₂(𝑠)×−𝑏₁𝑥₁(𝑠)−𝑥₂(𝑠)−𝑐₂𝑦₁(𝑠)1+𝛼₂𝑥₁(𝑠)+𝛽₂𝑦₁(𝑠)−𝑐₂𝑦₂(𝑠)1+𝛼₂𝑥₂(𝑠)+𝛽₂𝑦₂(𝑠)−𝑐₃𝑧₁(𝑠)1+𝛼₃𝑥₁(𝑠)+𝛽₃𝑧₁(𝑠)−𝑐₃𝑧₂(𝑠)1+𝛼₃𝑥₂(𝑠)+𝛽₃𝑧₂(𝑠)+sgn𝑦₁(𝑠)−𝑦₂(𝑠)×−𝑏₂𝑦₁(𝑠)−𝑦₂(𝑠)−ℎ₂𝑦₁(𝑠)𝑓₂+𝑔₂𝑧₁(𝑠)−ℎ₂𝑦₂(𝑠)𝑓₂+𝑔₂𝑧₂(𝑠)+𝑑₂𝑥₁(𝑠)1+𝛼₂𝑥₁(𝑠)+𝛽₂𝑦₁(𝑠)−𝑑₂𝑥₂(𝑠)1+𝛼₂𝑥₂(𝑠)+𝛽₂𝑦₂(𝑠)+sgn𝑧₁(𝑠)−𝑧₂(𝑠)×−𝑏₃𝑧₁(𝑠)−𝑧₂(𝑠)−ℎ₃𝑧₁(𝑠)𝑓₃+𝑔₃𝑦₁(𝑠)−ℎ₃𝑧₂(𝑠)𝑓₃+𝑔₃𝑦₂(𝑠)+𝑑₃𝑥₁(𝑠)1+𝛼₃𝑥₁(𝑠)+𝛽₃𝑧₁(𝑠)−𝑑₃𝑥₂(𝑠)1+𝛼₂𝑥₂(𝑠)+𝛽₂𝑧₂(𝑠)𝑑𝑠. (64) So 𝑑𝐸[𝑉(𝑡)]𝑥𝑑𝑡=𝐸sgn₁(𝑡)−𝑥₂(𝑡)×−𝑏₁𝑥₁(𝑡)−𝑥₂(𝑡)−𝑐₂𝑦₁(𝑡)1+𝛼₂𝑥₁(𝑡)+𝛽₂𝑦₁(𝑡)−𝑐₂𝑦₂(𝑡)1+𝛼₂𝑥₂(𝑡)+𝛽₂𝑦₂(𝑡)−𝑐₃𝑧₁(𝑡)1+𝛼₃𝑥₁(𝑡)+𝛽₃𝑧₁(𝑡)−𝑐₃𝑧₂(𝑡)1+𝛼₃𝑥₂(𝑡)+𝛽₃𝑧₂(𝑡)+sgn𝑦₁(𝑡)−𝑦₂(𝑡×)−𝑏₂𝑦₁(𝑡)−𝑦₂(𝑡)−ℎ₂𝑦₁(𝑡)𝑓₂+𝑔₂𝑧₁(𝑡)−ℎ₂𝑦₂(𝑡)𝑓₂+𝑔₂𝑧₂(𝑡)+𝑑₂𝑥₁(𝑡)1+𝛼₂𝑥₁(𝑡)+𝛽₂𝑦₁(𝑡)−𝑑₂𝑥₂(𝑡)1+𝛼₂𝑥₂(𝑡)+𝛽₂𝑦₂(𝑡)+sgn𝑧₁(𝑡)−𝑧₂(𝑡)×−𝑏₃𝑧₁(𝑡)−𝑧₂(𝑡)−ℎ₃𝑧₁(𝑡)𝑓₃+𝑔₃𝑦₁(𝑡)−ℎ₃𝑧₂(𝑡)𝑓₃+𝑔₃𝑦₂(𝑡)+𝑑₃𝑥₁(𝑡)1+𝛼₃𝑥₁(𝑡)+𝛽₃𝑧₁(𝑡)−𝑑₃𝑥₂(𝑡)1+𝛼₂𝑥₂(𝑡)+𝛽₂𝑧₂(𝑡)⩽−𝑏₁𝐸𝑥₁(𝑡)−𝑥₂(𝑡)−𝑐₂𝐸sgn𝑥₁(𝑡)−𝑥₂(𝑡)×𝑦₁(𝑡)1+𝛼₂𝑥₁(𝑡)+𝛽₂𝑦₁(𝑡)−𝑦₂(𝑡)1+𝛼₂𝑥₂(𝑡)+𝛽₂𝑦₂(𝑡)−𝑐₃𝐸sgn𝑥₁(𝑡)−𝑥₂(𝑡×𝑧₁)(𝑡)1+𝛼₃𝑥₁(𝑡)+𝛽₃𝑧₁(𝑡)−𝑧₂(𝑡)1+𝛼₃𝑥₂(𝑡)+𝛽₃𝑧₂(𝑡)−𝑏₂𝐸𝑦₁(𝑡)−𝑦₂(𝑡)−ℎ₂𝐸sgn𝑦₁(𝑡)−𝑦₂(𝑡)×𝑦₁(𝑡)𝑓₂+𝑔₂𝑧₁(𝑡)−𝑦₂(𝑡)𝑓₂+𝑔₂𝑧₂(𝑡)+𝑑₂𝐸sgn𝑦₁(𝑡)−𝑦₂(𝑡)×𝑥₁(𝑡)1+𝛼₂𝑥₁(𝑡)+𝛽₂𝑦₁(𝑡)−𝑥₂(𝑡)1+𝛼₂𝑥₂(𝑡)+𝛽₂𝑦₂(𝑡)−𝑏₃𝐸𝑧₁(𝑡)−𝑧₂(𝑡)−ℎ₃𝐸sgn𝑧₁(𝑡)−𝑧₂(𝑡)×𝑧₁(𝑡)𝑓₃+𝑔₃𝑦₁(𝑡)−𝑧₂(𝑡)𝑓₃+𝑔₃𝑦₂(𝑡)+𝑑₃𝐸sgn𝑧₁(𝑡)−𝑧₂(𝑡)×𝑥₁(𝑡)1+𝛼₃𝑥₁(𝑡)+𝛽₃𝑧₁(𝑡)−𝑥₂(𝑡)1+𝛼₂𝑥₂(𝑡)+𝛽₂𝑧₂(𝑡)⩽−𝑏₁𝐸𝑥₁(𝑡)−𝑥₂(𝑡)+𝑐₂𝐸𝑦₁(𝑡)1+𝛼₂𝑥₁(𝑡)+𝛽₂𝑦₁(𝑡)−𝑦₂(𝑡)1+𝛼₂𝑥₂(𝑡)+𝛽₂𝑦₂(𝑡)+𝑐₃𝐸𝑧₁(𝑡)1+𝛼₃𝑥₁(𝑡)+𝛽₃𝑧₁(𝑡)−𝑧₂(𝑡)1+𝛼₃𝑥₂(𝑡)+𝛽₃𝑧₂(𝑡)−𝑏₂𝐸𝑦₁(𝑡)−𝑦₂(𝑡)+ℎ₂𝐸𝑦₁(𝑡)𝑓₂+𝑔₂𝑧₁(𝑡)−𝑦₂(𝑡)𝑓₂+𝑔₂𝑧₂(𝑡)+𝑑₂𝐸𝑥₁(𝑡)1+𝛼₂𝑥₁(𝑡)+𝛽₂𝑦₁(𝑡)−𝑥₂(𝑡)1+𝛼₂𝑥₂(𝑡)+𝛽₂𝑦₂(𝑡)−𝑏₃𝐸𝑧₁(𝑡)−𝑧₂(𝑡)+ℎ₃𝐸𝑧₁(𝑡)𝑓₃+𝑔₃𝑦₁(𝑡)−𝑧₂(𝑡)𝑓₃+𝑔₃𝑦₂(𝑡)+𝑑₃𝐸𝑥₁(𝑡)1+𝛼₃𝑥₁(𝑡)+𝛽₃𝑧₁(𝑡)−𝑥₂(𝑡)1+𝛼₂𝑥₂(𝑡)+𝛽₂𝑧₂(𝑡)⩽−𝑏₁𝐸𝑥₁(𝑡)−𝑥₂(𝑡)+𝑐₂𝐸𝑦₁(𝑡)−𝑦₂(𝑡+𝑐₂)𝛼₂𝛽₂𝐸𝑥₁(𝑡)−𝑥₂(𝑡)+𝑐₃𝐸𝑧₁(𝑡)−𝑧₂(𝑡+𝑐₃)𝛼₃𝛽₃𝐸𝑥₁(𝑡)−𝑥₂(𝑡)−𝑏₂𝐸𝑦₁(𝑡)−𝑦₂(𝑡+ℎ₂)𝑓₂𝐸𝑦₁(𝑡)−𝑦₂(𝑡+ℎ₂)𝑔₂𝑓²₂𝐸𝑦₁(𝑡)𝐸𝑧₁(𝑡)−𝑧₂(𝑡+ℎ₂)𝑓₂𝐸𝑦₁(𝑡)−𝑦₂(𝑡)+𝑑₂𝐸𝑥₁(𝑡)−𝑥₂(𝑡+𝑑)²𝛽₂𝛼₂𝐸𝑦₁(𝑡)−𝑦₂(𝑡)−𝑏₃𝐸𝑧₁(𝑡)−𝑧₂(𝑡+ℎ)³𝑓₃𝐸𝑧₁(𝑡)−𝑧₂(𝑡+ℎ₃)𝑔₃𝑓²₃𝐸𝑧₁(𝑡)𝐸𝑦₁(𝑡)−𝑦₂(𝑡+ℎ₃)𝑓₃𝐸𝑧₁(𝑡)−𝑧₂(𝑡)+𝑑₃𝐸𝑥₁(𝑡)−𝑥₂(𝑡+𝑑₃)𝛽₃𝛼₃𝐸𝑧₁(𝑡)−𝑧₂(𝑡⩽)−𝑏₁+𝑐₂𝛼₂𝛽₂+𝑐₃𝛼₃𝛽₃+𝑑₂+𝑑₃×𝐸𝑥₁(𝑡)−𝑥₂(𝑡+)−𝑏₂+𝑐₂+2ℎ₂𝑓₂+𝑑₂𝛽₂𝛼₂+ℎ₃𝑔₃𝑓²₃𝐸𝑧³₁(𝑡)^1/3×𝐸𝑦₁(𝑡)−𝑦₂(𝑡+)−𝑏₃+𝑐₃+2ℎ₃𝑓₃+𝑑₃𝛽₃𝛼₃+ℎ₂𝑔₂𝑓²₂𝐸𝑦³₁(𝑡)^1/3×𝐸𝑧₁(𝑡)−𝑧₂(𝑡). (65) By Lemma 4𝐸𝑥³₁(𝑡)^1/3⩽𝑎₁+𝜎²₁𝑏₁,𝐸𝑦³₁(𝑡)^1/3⩽𝑎₂+𝑑₂/𝛼₂+𝜎²₂𝑏₂,𝐸𝑧³₁(𝑡)^1/3⩽𝑎₃+𝑑₃/𝛼₃+𝜎²₃𝑏₃. (66) Thus 𝑑𝐸[𝑉(𝑡)]⩽𝑑𝑡−𝑏₁+𝑐₂𝛼₂𝛽₂+𝑐₃𝛼₃𝛽₃+𝑑₂+𝑑₃𝐸𝑥₁(𝑡)−𝑥₂(𝑡+)−𝑏₂+𝑐₂+2ℎ₂𝑓₂+𝑑₂𝛽₂𝛼₂+ℎ₃𝑔₃𝑎₂+𝑑₂/𝛼₂+𝜎²₂𝑏₂𝑓²₃×𝐸𝑦₁(𝑡)−𝑦₂(𝑡+)−𝑏₃+𝑐₃+2ℎ₃𝑓₃+𝑑₃𝛽₃𝛼₃+ℎ₂𝑔₂𝑎₃+𝑑₃/𝛼₃+𝜎²₃𝑏₃𝑓²₂×𝐸𝑧₁(𝑡)−𝑧₂(𝑡𝑥)=𝐴𝐸₁(𝑡)−𝑥₂(𝑡)+𝐵𝐸𝑦₁(𝑡)−𝑦₂(𝑡𝑧)+𝐶𝐸₁(𝑡)−𝑧₂(𝑡). (67) Integrating both sides leads to 𝐸[𝑉(𝑡[𝑉(0+)]⩽𝐸)]_𝑡0𝐴𝐸𝑥₁(𝑠)−𝑥₂(𝑠𝑦)+𝐵𝐸₁(𝑠)−𝑦₂(𝑠)+𝐶𝐸𝑧₁(𝑠)−𝑧₂(𝑠)𝑑𝑠. (68) Therefore [𝑉(𝑡+0⩽𝐸)]_𝑡0𝐴𝐸𝑥₁(𝑠)−𝑥₂(𝑠)+𝐵𝐸𝑦₁(𝑠)−𝑦₂(𝑠𝑧)+𝐶𝐸₁(𝑠)−𝑧₂(𝑠[𝑉(0)𝑑𝑠⩽𝐸)]<∞. (69) Then we get 𝐸𝑥₁(𝑡),𝑦₁(𝑡),𝑧₁(𝑡)−𝑥₂(𝑡),𝑦₂(𝑡),𝑧₂(𝑡)=𝐸𝑥₁(𝑡)−𝑥₂(𝑡)²+𝑦₁(𝑡)−𝑦₂(𝑡)²+𝑧₁(𝑡)−𝑧₂(𝑡)^21/2⩽𝐸𝑥₁(𝑡)−𝑥₂(𝑡)+𝐸𝑦₁(𝑡)−𝑦₂(𝑡𝑧)+𝐸₁(𝑡)−𝑧₂(𝑡)∈𝐿¹[0,∞). (70) Therefore from Lemmas 11 and 12, we have lim_𝑡→+∞𝐸𝑥₁(𝑡),𝑦₁(𝑡),𝑧₁(𝑡)−𝑥₂(𝑡),𝑦₂(𝑡),𝑧₂(𝑡)=0. (71) This completes the proof of Theorem 13.

6. Extinction

Theorem 14. For any initial value (𝑥₀,𝑦₀,𝑧₀)∈𝑅³₊, the solution (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) of system (4) obeyslimsup_𝑡→∞ln𝑥(𝑡)𝑡⩽𝑎₁−𝜎²₁2,limsup_𝑡→∞ln𝑦(𝑡)𝑡⩽𝑎₂+𝑑₂𝛼₂−𝜎²₂2,limsup_𝑡→∞ln𝑧(𝑡)𝑡⩽𝑎₃+𝑑₃𝛼₃−𝜎²₃2. (72)

Proof. Define Lyapunov functions ln𝑥, ln𝑦, and ln𝑧, respectively, hence by Itô’s formula, we get ln𝑥(𝑡)=ln𝑥₀+𝑎₁−𝜎²₁2𝑡+^𝑡₀−𝑏₁𝑥(𝑠)−𝑐₂𝑦(𝑠)1+𝛼₂𝑥(𝑠)+𝛽₂𝑦(𝑠)−𝑐₃𝑧(𝑠)1+𝛼₃𝑥(𝑠)+𝛽₃𝑧(𝑠)(𝑡)𝑑𝑠,ln𝑦=ln𝑦₀+𝑎₂−𝜎²₂2𝑡+^𝑡₀−𝑏₂𝑦(𝑠)−ℎ₂𝑦(𝑠)𝑓₂+𝑔₂𝑧(𝑠)+𝑑𝑠_𝑡0𝑑₂𝑥(𝑠)1+𝛼₂𝑥(𝑠)+𝛽₂𝑦(𝑠)(𝑡)𝑑𝑠,ln𝑧=ln𝑧₀+𝑎₃−𝜎²₃2𝑡+^𝑡₀−𝑏₃𝑧(𝑠)−ℎ₃𝑧(𝑠)𝑓₃+𝑔₃𝑦(𝑠)+𝑑𝑠_𝑡0𝑑₃𝑥(𝑠)1+𝛼₃𝑥(𝑠)+𝛽₃𝑧(𝑠)𝑑𝑠. (73) Then we have ln𝑥(𝑡)⩽ln𝑥₀+𝑎₁−𝜎²₁2𝑡+𝜎₁𝐵₁(𝑡),ln𝑦(𝑡)⩽ln𝑦₀+𝑎₂+𝑑₂𝛼₂−𝜎²₂2𝑡+𝜎₂𝐵₂(𝑡),ln𝑧(𝑡)⩽ln𝑧₀+𝑎₃+𝑑₃𝛼₃−𝜎²₃2𝑡+𝜎₃𝐵₃(𝑡). (74) Dividing 𝑡 on the both sides and letting 𝑡→∞, we can derive limsup_𝑡→∞ln𝑥(𝑡)𝑡⩽𝑎₁−𝜎²₁2,limsup_𝑡→∞ln𝑦(𝑡)𝑡⩽𝑎₂+𝑑₂𝛼₂−𝜎²₂2,limsup_𝑡→∞ln𝑧(𝑡)𝑡⩽𝑎₃+𝑑₃𝛼₃−𝜎²₃2, (75)(76)(77) as required.

So we can obtain that if 𝑎₁−(𝜎²₁/2)<0, 𝑎₂+(𝑑₂/𝛼₂)−(𝜎²₂/2)<0, and 𝑎₃+(𝑑₃/𝛼₃)−(𝜎²₃/2)<0 hold, then for any initial value (𝑥₀,𝑦₀,𝑧₀)∈𝑅³₊, the solution (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) of system (4) will be extinct exponentially with probability one.

7. Numerical Simulations

At last, we numerically simulate the solution (𝑥(𝑡),𝑦(𝑡),𝑧(𝑡)) of system (4) to substantiate the analytical findings. By the method mentioned in [14], we consider the discretization equation: 𝑥_𝑘+1=𝑥_𝑘+𝑥_𝑘×𝑎₁−𝑏₁𝑥_𝑘−𝑐₂𝑦_𝑘1+𝛼₂𝑥_𝑘+𝛽₂𝑦_𝑘−𝑐₃𝑧_𝑘1+𝛼₃𝑥_𝑘+𝛽₃𝑧_𝑘Δ𝑡+𝜎₁𝑥_𝑘𝜉_𝑘√1Δ𝑡+2𝜎²₁𝑥_𝑘𝜉²_𝑘𝑦−1Δ𝑡,_𝑘+1=𝑦_𝑘+𝑦_𝑘×𝑎₂−𝑏₂𝑦_𝑘−ℎ₂𝑦_𝑘𝑓₂+𝑔₂𝑧_𝑘+𝑑₂𝑥_𝑘1+𝛼₂𝑥_𝑘+𝛽₂𝑦_𝑘Δ𝑡+𝜎₂𝑦_𝑘𝜂_𝑘√1Δ𝑡+2𝜎²₂𝑦_𝑘𝜂²_𝑘𝑧−1Δ𝑡,_𝑘+1=𝑧_𝑘+𝑧_𝑘×𝑎₃−𝑏₃𝑧_𝑘−ℎ₃𝑧_𝑘𝑓₃+𝑔₃𝑦_𝑘+𝑑₃𝑥_𝑘1+𝛼₃𝑥_𝑘+𝛽₃𝑧_𝑘Δ𝑡+𝜎₃𝑧_𝑘𝜁_𝑘√1Δ𝑡+2𝜎²₃𝑧_𝑘𝜁²_𝑘−1Δ𝑡, (78) where 𝜉_𝑘, 𝜂_𝑘, and 𝜁_𝑘 are Gaussian random variables that follow 𝑁(0,1).

For example, in Figure 1, we choose initial value (𝑥₀,𝑦₀,𝑧₀)=(0.5,0.5,0.5) and parameters satisfying conditions of Theorem 7; the system is stochastically permanent. In Figure 2, (a) is the solution to the deterministic system (4), (b) is the solution to the stochastic system (4), and (c) is the phase diagram of the stochastic system (4). Comparing (a) with (b), we find that with small environmental noise, the stochastic system is getting more similar to the deterministic. In Figure 3, (a) is the solution to the deterministic system (4) and (b) is the solution to the stochastic system (4). Comparing (a) with (b), we find that the sufficiently large environmental noise make the stochastic system extinct. The phase diagram of (𝑥,𝑦,𝑧) is displayed in (c).

Figure 1 

Solutions of system (4) for (𝑥0,𝑦0,𝑧0) = (0.5,0.5,0.5),𝑎1 = 2;𝑎2=1;𝑎3 = 1.5;𝑏1=3;𝑏2=2;𝑏3=2;𝑐2 = 0.2;𝑐3=0.1;𝑑2=1;𝑑3 = 1.5;𝛼2=1;𝛽2=0.5;𝛼3=0.6;𝛽3=0.8;ℎ2=0.5;ℎ3=1;𝑓2=1;𝑓3=0.3;𝑔2 = 0.5;𝑔3=1;𝜎1=0.05;𝜎2 = 0.05;𝜎3=0.05.

(a)

(b)

Figure 2 

Solutions of system (4) for (𝑥0,𝑦0,𝑧0) = (0.6,0.5,0.5),𝑎1=2;𝑎2=1;𝑎3 = 1.5;𝑏1=2.5;𝑏2 = 2;𝑏3=1.5;𝑐2=0.2;𝑐3 = 0.1;𝑑2=1.5;𝑑3=1;𝛼2 = 1;𝛽2=0.5;𝛼3=1;𝛽3=0.8;ℎ2=0.5;ℎ3=1;𝑓2=1;𝑓3=0.3;𝑔2=1;𝑔3=0.5;𝜎1=0.05;𝜎2 = 0.05;𝜎3=0.05.

(a)

(b)

(c)

Figure 3 

Solutions of system (4) for (𝑥0,𝑦0,𝑧0)=(1.5,1,1),𝑎1 = 1.1;𝑎2=1.2;𝑎3=1.5;𝑏1=0.8;𝑏2=1.1;𝑏3=1.2;𝑐2=0.02;𝑐3 = 0.01;𝑑2=1.2;𝑑3=1;𝛼2=0.7;𝛽2=0.5;𝛼3 = 0.8;𝛽3=0.5;ℎ2=0.8;ℎ3=0.5;𝑓2=1;𝑓3=0.3;𝑔2=0.8;𝑔3=0.5;𝜎1 = 1.8;𝜎2=1.8;𝜎3=1.8.

8. Discussion

This paper has been devoted to dynamics of a stochastic cooperative predator-prey system with Beddington-DeAngelis functional. Firstly we show that, although the coefficients in the model neither satisfy the linear growth condition nor local Lipschitz continuity, the stochastic model has a globally positive solution. Then we know that the positive solution is stochastically bounded. Moreover, under some conditions, we analyze global asymptotic stability of the positive solutions. We can find that the stochastic model will preserve this nice property provided that the noise is sufficiently small. Some meaningful questions deserve further investigation. One way we can consider colored noise in the models owing to sudden environmental changes caused by seasons or other reasons. Moreover, it is worth to rebuild our model with some parameters, not the ones studied above, which are also subject to stochastic excitation. The ecosystem we considered has some limitations, for example, both predators (𝑦 and 𝑧) with the common prey (𝑥) and we only established a model with a mutual cooperation between both predators. It is worth to discuss mechanics of the population contained predators with a mutual competition in other ecosystems.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This research was supported by the National Natural Science Fund of China (NNSF: 11102132).