Bull Math Biol DOI 10.1007/s11538-015-0086-4 ORIGINAL ARTICLE

Survival and Stationary Distribution Analysis of a Stochastic Competitive Model of Three Species in a Polluted Environment Yu Zhao1,2 · Sanling Yuan1 · Junling Ma3

Received: 29 November 2014 / Accepted: 8 May 2015 © Society for Mathematical Biology 2015

Abstract In this paper, we develop and study a stochastic model for the competition of three species with a generalized dose–response function in a polluted environment. We first carry out the survival analysis and obtain sufficient conditions for the extinction, non-persistence, weak persistence in the mean, strong persistence in the mean and stochastic permanence. The threshold between weak persistence in the mean and extinction is established for each species. Then, using Hasminskii’s methods and a Lyapunov function, we derive sufficient conditions for the existence of stationary distribution for each population. Numerical simulations are carried out to support our theoretical results, and some biological significance is presented. Keywords Stochastic competition model · Polluted environment · Itô’s formula · Survival analysis · Stationary distribution Mathematics Subject Classification

B

34K50 · 92B05 · 60J28

Sanling Yuan [email protected] Junling Ma [email protected]

1

College of Science, University of Shanghai for Science and Technology, Shanghai 200093, China

2

School of Mathematics and Computer Science, Ningxia Normal University, Ningxia 756000, Guyuan, China

3

Department of Mathematics and Statistics, University of Victoria, Victoria, BC V8W 3R4, Canada

123

Y. Zhao et al.

1 Introduction In recent decades, rapid industrial development and other human activities introduced a mass of toxicants into ecosystems. These toxicants have seriously threatened the survival of the exposed populations such as human, fishes and planktons. For example, the chemical wastes and sewage leakages from factories contain complex and highly toxic organic compositions in high density and are difficult to degrade, causing species extinctions as toxins destroy their natural habitats and community structures. Therefore, it is important to assess environmental toxicity and to develop optimal policies to reduce environmental and economic losses, which require qualitative and quantitative estimates of the survival risk of the species in a polluted environment. Using mathematical models to estimate the effects of a toxin on a population has recently attracted considerable interests (Mezquita et al. 1999; Yang and Liu 2010; Hallam and De Luna 1984; Wang 2010; Chen and Chen 1993). Species commonly compete with other species for, e.g., foods and habitats. It is thus important to consider the toxins’ effect on the survival of competing species. Since Hallam et al. (1983) first proposed a single-species model in a polluted environment, many remarkable studies about the persistence of two or more species in toxic environments have been carried out. For example, Liu and Ma (1991) studied a Lotka– Volterra model of two species in a polluted environment and obtained the threshold for persistence and extinction of each species. This model was then extended to a competitive Lotka–Volterra model of three species with a toxin by Dong and Wang (1997), which takes the form ⎧ 3 dxi (t) ⎪ ⎨ dt = xi (t)[ri0 − ri1 C0 (t) − j=1 aij x j (t)], i = 1, 2, 3, dC0 (t) = kC E (t) − (g + m)C0 (t), dt ⎪ ⎩ dC E (t) = −hC (t) + u(t), E dt

(1)

with initial conditions xi (0) = xi0 > 0, i = 1, 2, 3 and C0 (0) = 0, C E (0) = 0, where xi (t)(i = 1, 2, 3) is the density of the ith species at time t, C0 (t) is the toxin concentration in the organism at time t, and C E (t) is the toxin concentration in the environment at time t. The intrinsic growth rate of the ith species in the absence of toxicant, denoted as ri0 , is assumed to be constant. The dose–response rate of species i to the organismal toxicant concentration, denoted as ri1 , and the action of species j upon the growth rate of species i, denoted as aij for i, j = 1, 2, 3, are assumed to be positive constants. In particular, aii represents the intraspecific competition coefficient of species i. The toxin uptake rate per unit biomass, namely k, the rate of toxin loss in the environment h, the net ingestion rate g and the depuration rate m are also assumed to be positive constants. The environmental toxin input rate is u(t). By the method of integral average, Dong and Wang (1997) obtained the survival thresholds for each species. Other related research can also be found in , e.g., He and Wang (2009), Wang and Ma (1994), Duan et al. (2004), Ma and Hallam (1987), Pan et al. (2000), Zhao

123

Survival and Stationary Distribution Analysis of a. . .

et al. (2009), Hallam et al. (1979), Hallam and Ma (1986) and Dubeya and Hussainb (2006). These studies are all based on deterministic models. However, the growth of populations is inevitably affected by environmental noise. In fact, it has been noted that the deterministic models in ecology are based on the assumption that, when population sizes approach infinity, random environmental fluctuations are small and can be ignored (Melbourne and Hastings 2008; Bandyopadhyay and Chattopadhyay 2005). But in reality, environmental stochasticity may involve the variations of factors such as necessary nutrients, climate, acidity and habitats, which may affect the growth rate. May (1973) has claimed that owing to random noises, the growth rates in population systems should be stochastic, and as a result, the solution of the system will not tend to a steady positive point, but fluctuate around some average values. Hence, stochastic models are generally more realistic than their deterministic counterparts, and it is necessary to incorporate environmental stochasticity into model (1). Notice that the parameter ri0 in model (1) represents the intrinsic growth rate of the ith species. To model the random changes in the growth rate produced by the randomly varying environment, we use the following mean-reverting Ornstein– Uhlenbeck process (Duffie 1996; Wu et al. 2008): dri0 (t) = αi (rie − ri0 (t))dt + βi dB i (t), i = 1, 2, 3,

(2)

where αi , rie , βi are all positive constants. αi is the reversion rate, rie is the mean reversion level or long-run equilibrium of growth rate ri0 , βi is the intensity of volatility, and B i (t) is a standard Brownian motion. Combining model (1) with Eq. (2), we derive the following stochastic system: ⎧  ⎪ dxi (t) = xi (t)[ri0 (t) − ri1 H (C0 (t)) − 3j=1 aij x j (t)]dt, i = 1, 2, 3 ⎪ ⎪ ⎨ dri0 (t) = αi (rie − ri0 (t))dt + βi dB i (t), i = 1, 2, 3, dC0 (t) ⎪ ⎪ dt = kC E (t) − (g + m)C0 (t), ⎪ ⎩ dC E (t) = −hC E (t) + u(t), dt

(3)

where H (C0 ), which is a non-decreasing function with respect to C0 with H (0) = 0, is the dose–response function (see Liu and Wang 2011; Butler 1979; Filov et al. 1979 for the reason why we use this generalized-type function), and it is assumed that |H (C0 (t))| ∈ L 1 [0, ∞). Survival analysis is one of main current approaches to deal with stochastic population model, which reveals the persistence or extinction of populations in a random environment. On the other hand, in order to study the statistic characteristic of the long-term behaviors of the sample trajectories, it is useful to further analyze the stationary distribution, which is an important measure of stochastic population model from the ecological view point. Recently, the stationary distribution analysis has attracted considerable interests. For example, Mao (2011) proved the existence of a unique stationary distribution and developed a useful method to compute the mean and variance of the stationary distribution. Danga et al. (2014) considered a population dynamics of two species described by Kolmogorov systems of competitive type under telegraph

123

Y. Zhao et al.

noise and established the existence of an invariant (or a stationary) probability measure. One may also refer to Gard (1992), Li (2013), Rudnicki and Pichor (2007), Hasminskii (1980), Mandal and Banerjee (2012), Ji et al. (2011), Jiang et al. (2012), Liu and Wang (2012) and Rudnicki (2003) for more related studies. In this paper, we will devote our main attention to the analysis of model (3). We first carry out the survival analysis for the system. Then, we prove the existence of its stationary distribution. Notice that model (3) is a non-autonomous system. As far as we know, few works have been performed for the existence of the stationary distribution of a non-autonomous model in the literature. The organization of this paper is as follows. The model and some basic preliminaries are presented in the next section. In Sect. 3, we perform the survival analysis for the model and obtain the sufficient conditions for extinction, non-persistence, weak persistence in the mean, strong persistence in the mean and stochastic permanence. The threshold between weak persistence in the mean and extinction will be established for each population. Using of Hasminskii’s methods and constructing Lyapunov function, in Sect. 4, sufficient conditions ensuring the existence of stationary distribution of the model are obtained. To illustrate our theoretical results, some numerical simulations are carried out in Sect. 5. Finally, a discussion on the biological significance of the obtained results is presented in the last section.

2 Some Basic Preliminaries From now on, unless otherwise specified, let (, F , {Ft }t≥0 , P) be a complete probability space with a filtration {Ft }t≥0 satisfying the usual conditions (i.e., it is right continuous and F0 contains all P-null sets). Notice that Eq. (2) has the following explicit form solution: ri0 (t) = rie + (ri0 (0) − rie )e

−αi t

 + βi

t

e−αi (t−s) dB i (s), i = 1, 2, 3.

(4)

0

t β2

Since the term βi 0 e−αi (t−s) dB i (s) follows the normal distribution N 0, 2αi i 1 −  e−2αi t , it is a.e. equal to its version

βi

(1 − e−2αi t ) dBi (t) , 2αi dt

where Bi (t) is a standard Brownian motion. Thus, (4) can be almost surely written as the following form: ri0 (t) = rie + (ri0 (0) − rie )e−αi t + σi (t)

dBi (t) , i = 1, 2, 3, dt

(5)

 i , i = 1, 2, 3. Submitting (5) into the where σi (t) = σi (1 − e−2αi t ) and σi = √β2α i first three equations in model (3), we can get the following stochastic model:

123

Survival and Stationary Distribution Analysis of a. . .

  ⎧ 3 −αi t − r H (C (t)) − ⎪ dx (t) = x (t) r + (r (0) − r )e a x (t) dt i i ie i0 ie i1 0 ij j ⎪ j=1 ⎪ ⎨ + xi (t)σi (t)dBi (t), i = 1, 2, 3, (6) dC0 (t) ⎪ = kC ⎪ dt E (t) − (g + m)C 0 (t), ⎪ ⎩ dC E (t) = −hC E (t) + u(t), dt with initial data xi (0) = xi0 > 0, i = 1, 2, 3 and 0 ≤ C0 (0) ≤ 1, 0 ≤ C E (0) ≤ 1,

(7)

Note that model (6) is a five-dimensional model and C0 (t) and C E (t) can be explicitly solved (as functions of u(t)) from the last two equations. So we need only to consider the following three-dimensional model:  3 dxi (t) = xi (t) rie + (ri0 (0) − rie )e−αi t − ri1 H (C0 (t)) −

j=1

+ xi (t)σi (t)dBi (t), i = 1, 2, 3,

 aij x j (t) dt (8)

with initial data xi (0) = xi0 > 0,

i = 1, 2, 3.

(9)

For the sake of convenience and simplicity, we introduce the following notations:       σ2  a11 r1e − 21 r11   a11 a12 a13      2 A =  a21 a22 a23  , R =  a21 r2e − σ2 r21  , 2    a31 a32 a33  2  a31 r3e − σ3 r31  2       σ2  r1e − 21 a12 a13   r11 a12 a13      2 A1 =  r2e − σ2 a22 a23  ,  A1 =  r21 a22 a23  , 2    r31 a32 a33  2  r3e − σ3 a32 a33  2       σ2  a11 r1e − 21 a13   a11 r11 a13      2 A2 =  a21 r21 a23  , A2 =  a21 r2e − σ2 a23  ,  2    a31 r31 a33  2  a31 r3e − σ3 a33  2      σ2   a11 a12 r1e − 21   a11 a12 r11      2 A3 =  a21 a22 r21  . A3 =  a21 a22 r2e − σ2  ,  2   a31 a32 r31  2   a31 a32 r3e − σ3  2 Now that model (8) describes the situation that three species compete with each other in a polluted environments, it is very important to find the conditions under which each species will be persistent or go to extinction. Denote

123

Y. Zhao et al.

 f (t) =

1 t

 0

t

f (s)ds, f ∗ = lim sup f (t), f ∗ = lim inf f (t). t→∞

t→∞

The following definitions about extinction, non-persistence, weakly persistence, strongly persistence in the mean and stochastically permanent can be found in Wang (2010). Definition 2.1 (i) The population x(t) is said to go to extinction if limt→+∞ x(t) = 0. (ii) The population x(t) is said to be stochastic non-persistence in the mean if x∗ = 0, a.s. (iii) The population x(t) is said to be stochastic weak persistence in the mean if x∗ > 0, a.s. (iv) The population x(t) is said to be stochastic strong persistence in the mean if x∗ > 0, a.s. Definition 2.2 Model (8) is said to be stochastically permanent if for any  ∈ (0, 1), there exists a pair of positive constant γ = γ () and χ = χ () such that for any initial data (9), the solution x(t) of model (8) has the property that lim inf P {|x(t)| ≤ χ } ≥ 1 − , t→∞

lim inf P {|x(t)| ≥ γ } ≥ 1 − . t→∞

It is easy to have the following basic results for model (6) or (8) (see Wang 2010 or Liu et al. 2011 for details). Lemma 2.1 For model (6), if 0 < k < g + m and supt≥0 u(t) ≤ h, then 0 ≤ C0 (t) < 1, 0 ≤ C E (t) < 1 for all t ∈ R+ . Lemma 2.2 For any given initial value (9), there is a unique solution x(t) = (x1 (t), x2 (t), x3 (t)) of model (8) on t > 0 and the solution will remain in R3+ with probability one, namely x(t) ∈ R3+ for all t ≥ 0 almost surely. From now on, we impose 0 < k < g + m and supt≥0 u(t) ≤ h on model (3) or (8). The following two lemmas, whose proofs can be found in Liu et al. (2011), will play an important role in the survival analysis of model (8). Lemma 2.3 For model (8), we have lim supt→∞ [ ln xti (t) ] ≤ 0, i = 1, 2, 3. Lemma 2.4 Suppose x(t) ∈ C[ × R3+ , R0+ ], where R0+ = {a|a > 0, l ∈ R}. (i) If there are positive constants λ0 , λ and T ≥ 0, such that  ln x(t) ≤ λt − λ0

0

t

x(s)ds +

n   i=1

t

σi (s)dBi (s),

0

for t > T , where σi is a bounded function, 1 ≤ i ≤ n, then x∗ ≤

123

λ λ0

a.s.

Survival and Stationary Distribution Analysis of a. . .

(ii) If there are positive constants λ0 , λ and T ≥ 0, such that  ln x(t) ≥ λt − λ0

t

x(s)ds +

0

n   i=1

t

σi (s)dBi (s),

0

for t > T , where σi is a bounded function, 1 ≤ i ≤ n, then x∗ ≥

λ λ0

a.s.

3 Survival Analysis of the Model In the previous section, Lemma 2.2 shows that the solution of model (8) with a positive initial data will remain positive. This nice property provides us with a great opportunity to discuss in more detail how the solution varies in R3+ . We first give the following fundamental assumptions. Assumption 3.1 rie >

σi2 2 .

Assumption 3.2 Ai > 0, i = 1, 2, 3. Assumption 3.3 R, R21 > 0; R11 , R31 ≥ 0. Assumption 3.4 (1) (Ai ) j j > 0, i, j = 1, 2, 3; (2) A12 , A21 , A23 , A32 > 0, A13 , A31 < 0; i > 0, i = 1, 2, 3, (  A2 )33 > 0. Here Φij denotes the complement minor of (3) A the element φij in the determinant Φ = (φij )3×3 . Assumption 3.1 implies that, under the effect of environment noise, each species can live alone in a non-polluted environment provided the intensity of the noise is small enough such that rie > σ2

σi2 2 . Under Assumptions 3.1 and 3.2, we can easily prove

from rie > 2i and Ai > 0, i = 1, 2, 3 that A > 0; this implies that when there is no environment perturbations (i.e., σi = 0), model (8) has a globally stable positive equilibrium in a non-polluted environment. Assumption 3.3 implies that r1e − r11

σ12 2

r2e − ≥ r21

σ22 2

r3e − ≥ r31

σ32 2

r1e − and r11

σ12 2

r3e − > r31

σ32 2

,

(10)

which in turn implies that the order of antitoxic and antinoise ability of population (from strongest to weakest) is x1 , x2 , x3 . Other conditions in the above assumptions will be needed in the following survival analysis. Now, we are in a position to establish the following threshold theorems. Denote r1e − μ1 = r11

σ12 2

, μ2 =

(A2 )33 A3 , μ3 = .   ( A2 )33 A3

We can easily prove that μ1 ≥ μ2 > μ3 .

(11)

123

Y. Zhao et al.

In fact, simple computations show that

μ1 − μ2 = =

r1e − r11

σ12 2

(A2 ) − 33 =  A2 33

a11

 r1e −



 r21 − r2e −

 A2 33 r11 

σ12 2

σ22 2



 r11

a11 R31

 A2 33 r11 

and that μ2 − μ3 =

A2 )33 A3 − A3 (  (A2 )33 A3 (A2 )33  A33 R − = = .      ( A2 )33 ( A2 )33 A3 ( A2 )33  A3 A3

Noting that Aii = (Ai )ii > 0, it follows from the assumptions above that (11) holds. Theorem 3.1 If H (C0 )∗ > μi , then xi will go to extinction a.s. (i = 1, 2, 3). Proof For species x1 , define a function V1 = ln x1 (t), applying Itô’s formula to the first equation of model (8) and integrating both sides from 0 to t, we obtain   σ12 1 x1 (t) ln = r1e − − r11 H (C0 (t)) − a11 x1 (t) − a12 x2 (t) t x1 (0) 2  1 t − a13 x3 (t) + f 1 (t) + σ1 (s)dB1 (s), t 0

(12)

where f 1 (t) =

 σ12 e−2α1 t (r10 (0) − r1e )

+ 1 − e−α1 t . 2 tα1

Obviously, lim f 1 (t) = 0. Noting that t→∞

t 0

σ1 (s)dB1 (s) = σ1

t  0

(1 − e−2α1 s )d

B1 (s) is a local martingale, it follows from the strong law of large number of local martingale (Theorem 3.4, Mao 2007) that 1 t→∞ t



t

lim

σ1 (s)dB1 (s) = 0.

0

Taking upper limit on both sides of (12) and making use of H (C0 )∗ > μ1 yields    σ12 x1 (t) ≤ r1e − lim sup ln − r11 H (C0 )∗ < 0, x1 (0) 2 t→∞ 

which implies limt→∞ x1 (t) = 0 a.s. For species x2 , define a function V2 = ln x2 (t), applying Itô’s formula to the second equation of model (8) and integrating both sides from 0 to t, we can get

123

Survival and Stationary Distribution Analysis of a. . . 1 t

ln

x2 (t) x2 (0)

 = r2e −

σ22 2



− r21 H (C0 (t)) − a21 x1 (t) − a22 x2 (t) t − a23 x3 (t) + f 2 (t) + 1t 0 σ2 (t)dB2 (s),

(13)

where f 2 (t) =

 σ22 e−2α2 t (r20 (0) − r2e )

+ 1 − e−α2 t , 2 tα2

Obviously, lim f 2 (t) = 0. Computing (13) × a11 − (12) × a21 yields t→∞

a11 x2 (t) a21 x1 (t) ln − ln = (A2 )33 − (  A2 )33 H (C0 (t)) − A33 x2 (t) t x2 (0) t x1 (0) − A32 x3 (t) + a11 f 2 (t) − a21 f 1 (t)   a11 t a21 t + σ2 (s)dB2 (s) − σ1 (s)dB1 (s). t 0 t 0 (14) Analyzing as that for (12), it follows from Lemma 2.3 that if H (C0 )∗ > μ2 , then lim sup t→∞

a

11

t

ln

x2 (t)  ≤ (A2 )33 − (  A2 )33 H (C0 )∗ < 0, x2 (0)

(15)

which implies limt→∞ x2 (t) = 0 a.s. For species x3 , define a function V3 = ln x3 (t), applying Itô’s formula to the third equation of model (8) and integrating both sides from 0 to t yields   σ32 1 x3 (t) ln = r3e − − r31 H (C0 (t)) − a31 x1 (t) − a32 x2 (t) t x3 (0) 2  1 t − a33 x3 (t) + f 3 (t) + σ3 (s)dB3 (s), t 0

(16)

where f 3 (t) =

σ32 e−2α3 t (r30 (0) − r3e ) + (1 − e−α3 t ), 2 tα3

Obviously, lim f 3 (t) = 0. Multiplying both sides of (12), (13) and (16), respect→∞ tively, by m , n and −1, and then adding these three equalities, we get that

123

Y. Zhao et al.

m  x1 (t)  n x2 (t) A3 −  1 x3 (t) A3 H (C0 (t)) A ln = ln + ln + − x3 (t) t x3 (0) t x1 (0) t x2 (0) A33 A33 t A13 0 σ1 (s)dB1 (s) +m  f 1 (t) +  n f 2 (t) − f 3 (t) + A33 t t t σ (s)dB (s) A23 0 σ2 (s)dB2 (s) 3 3 − 0 , (17) + A33 t t 

a11 m  + a21 n = a31 , > which implies m  = − AA13 33  + a22 n = a32 , a12 m > 0. Analyzing as that for (14), it follows from Lemma 2.3 that if 0, n = AA23 33 H (C0 )∗ > μ3 , then where m , n satisfy the equations

lim sup t→∞

1 t

ln

3 H (C0 )∗ A3 − A x3 (t)  ≤ < 0, x3 (0) A33

from which we know that limt→∞ x3 (t) = 0 a.s.

(18) 

Theorem 3.2 If H (C0 )∗ = μi , then xi will be nonpersistent in the mean a.s. (i=1,2,3). Proof For species x1 , noting that limt→∞ f 1 (t) = 0, we have that for any  > 0, there exists a T > 0 such that H (C0 (t)) ≥ H (C0 )∗ − 2 and f 1 (t) ≤ r112  for all t > T . Therefore, we have from (12) that if H (C0 )∗ = μ1 , then    σ12 x1 (t) ln ≤ r1e − t − r11 H (C0 (t))∗ t + r11 t − a11 x1 (t)t x1 (0) 2  t + σ1 (s)dB1 (s) 0  t  t ≤ r11 t − a11 x1 (s)ds + σ1 (s)dB1 (s). 

0

(19)

0

According to Lemma 2.4, we must have that x1 ∗ ≤ ra1111 . Since  is arbitrary, we have that x1 ∗ = 0. For species x2 , from Lemma 2.3 and limt→∞ f i (t) = 0(i = 1, 2), we have that for any  > 0, there exists a T > 0 such that a21 ln xx101 ≤ 3 , a11 f 2 (t) − a21 f 1 (t) ≤ 3  and H (C0 (t)) ≥ H (C0 )∗ − 3(  for all t > T . Substituting these inequalities A2 )33 into (14) and making use of Lemma 2.4, we can have that if H (C0 )∗ = μ2 , then x2 ∗ ≤ /(  A2 )33 . Since  is arbitrary, we have that x2 ∗ = 0. For x3 , using Lemmas 2.3 and 2.4, we can similarly prove from (17) that if  H (C0 )∗ = μ3 , then x3 ∗ = 0.

123

Survival and Stationary Distribution Analysis of a. . .

Theorem 3.3 If H (C0 )∗ < μi , then xi will be weakly persistent in the mean a.s. (i=1,2,3). Proof We first consider species x1 . If H (C0 )∗ < μ1 , it then follows from Lemma 2.3 and (12) that a11 x1 ∗ + a12 x2 ∗ + a13 x3 ∗ ≥ r1e −

σ12 − r11 H (C0 )∗ > 0. 2

(20)

We claim that x1 ∗ > 0 a.s. In fact, for any ω ∈ {x1 (ω)∗ = 0}, we must have that x2 (ω)∗ > 0 or x3 (ω)∗ > 0, which has three cases in all. Case (i) x2 (ω)∗ > 0 and x3 (ω)∗ = 0. From (A1 )33 , (A2 )33 > 0 and R31 ≥ 0 we have that r1e − a12 < a22 r2e −

σ12 2 σ22 2

r1e − a11 r11 < and ≤ a21 r21 r2e −

σ12 2 σ22 2

,

11 which imply that rr21 . ∈ 0, aa11 21

 (a) If rr11 ∈ 0, aa12 , then (  A1 )33 ≤ 0. Computing (12) ×a22 −(13) × a12 leads to 21 22 a22 t

1 t

ln

x1 (t,ω) x1 (0)

−

a12 t

ln

x2 (t,ω) x2 (0)

= (A1 )33 − (  A1 )33 H (C0 (t)) − A33 x1 (t, ω) + A31x3 (t, ω) + a22 f 1 (t) − a12 f 2 (t) t t + at22 0 σ1 (s)dB1 (s) − at12 0 σ2 (s)dB2 (s).

Taking upper limit on both sides of the above equality, it follows from lim supt→∞  ∗ ∗ ln xx11(t,ω) (0) ≤ 0, x 1 (ω) = 0 and x 3 (ω) = 0 that a12 lim sup t→∞

 ln x (t, w)  2 ≤ −(A1 )33 + (  A1 )33 H (C0 )∗ < 0. t

Thus, limt→∞ x2 (t, ω) = 0, which contradicts with x2 (ω)∗ > 0. Therefore, ∈ / (0, aa12 ]. 22

 r11 A1 33 > 0. Computing (13) ×r11 − (12) ×r21 yields , a11 ), then  (b) If r21 ∈ ( aa12 22 a21 r11 r21

r11 x2 (t, ω) r21 x1 (t, ω) ln − ln = −R31 + (  A2 )33 x1 (t, ω) − (  A1 )33 x2 (t, ω) t x2 (0) t x1 (0) − ( A1 )32 x3 (t, ω) − r21 f 1 (t) + r11 f 2 (t) t t σ2 (s)dB2 (s) 0 σ1 (s)dB1 (s) + r11 0 . − r21 t t (21) Noting that lim supt→∞

1 t

ln

x1 (t,ω)  x1 (0)

≤ 0, x1 (ω)∗ = 0 and x3 (ω)∗ = 0 imply

that for any ε > 0, there exists a T such that when t > T, r21 t ln

x1 (t,ω) x1 (0)

< ε/4, and

123

Y. Zhao et al.

  max |(  A2 )33 x1 (t, ω)|, |(  A1 )32 x3 (t, ω)|, |r11 f 2 (t) − r21 f 1 (t)| < ε/4, it then follows from (21) that when t > T , ln

 t x2 (t, ω) −1 −1  (−R31 + ε)t − r11 ( A1 )33 x2 (s, ω)ds ≤ r11 x2 (0) 0  t  t −1 − r11 r21 σ1 (s)dB1 (s) + σ2 (s)dB2 (s). 0

0

According to Lemma 2.4, we obtain that x2 (ω)∗ ≤ (−R31 + ε)/(  A1 )33 . Then, it follows from the arbitrariness of ε and −R31 ≤ 0 that x2 (ω)∗ = 0, which ∈ / aa12 , a11 . In other words, x1 ∗ > 0 contradicts with x2 (ω)∗ > 0. Therefore, rr11 21 22 a21 a.s. Case (ii) x2 (ω)∗ = 0 and x3 (ω)∗ > 0. From (A3 )22 > 0, (A1 )22 > 0, R21 > 0, we have that r1e − a13 < a33 r3e − which implies that (a) If

r11 r31

r11 r31

σ12 2 σ32 2

r1e − a11 r11 < and < a31 r31 r3e −

σ12 2 σ32 2

,

∈ (0, aa11 ). 31

∈ (0, aa13 ], then (  A1 )22 ≤ 0. Computing (12) ×a33 − (13) ×a13 yields 33

a33 x1 (t, ω) a13 x3 (t, ω) ln − ln = (A1 )22 −(  A1 )22 H (C0 (t)) − A22 x1 (t, ω) t x1 (0) t x3 (0) − A21 x2 (t, ω) + a33 f 1 (t) − a13 f 3 (t)  a33 t + σ1 (s)dB1 (s) t 0  a13 t − σ3 (s)dB3 (s). (22) t 0 It follows from lim supt→∞ that a13 lim sup t→∞

r11 r31

1 t

ln

x1 (t,ω)  x1 (0)

≤ 0, x1 (ω)∗ = 0 and x2 (ω)∗ = 0

 ln x (t, ω)  3 ≤ −(A1 )22 + (  A1 )22 H (C0 )∗ < 0. t

That is to say, limt→∞ x3 (t, ω) = 0, which contradicts with x3 (ω)∗ > 0. So ∈ / (0, aa13 ]. 33

123

Survival and Stationary Distribution Analysis of a. . .

(b) If

r11 r31

∈



a13 a11 a33 , a31

, then (  A1 )22 > 0. Computing (16) ×r11 − (12) ×r31 leads to,

r11 x3 (t, ω) r31 x1 (t, ω) A3 )22 x1 (t, ω)+(  A3 )21 x2 (t, ω) ln − ln = −R21 + (  t x3 (0) t x1 (0) − ( A1 )22 x3 (t, ω) + r11 f 3 (t) − r31 f 1 (t)  r31 t − σ1 (s)dB1 (s) t 0  r11 t + σ3 (s)dB3 (s). (23) t 0   ≤ 0, x1 (ω)∗ = 0 and x2 (ω)∗ = 0 imply Noting that lim supt→∞ 1t ln xx11(t,ω) (0) that for any ε > 0, there exists a T such that when t > T , r31 x1 (t, ω) ln < ε/4, t x1 (0) and   A3 )21 x2 (t, ω)|, |r11 f 3 (t) − r31 f 1 (t)| < ε/4, max |(  A3 )22 x1 (t, ω)|, |(  it then follows from (23) that when t > T ,  t x3 (t, ω) −1 −1  ≤ r11 (−R21 + ε)t − r11 ( A1 )22 ln x3 (s, ω)ds x3 (0) 0  t  t −1 − r11 r31 σ1 (s)dB1 (s) + σ3 (s)dB3 (s). 0

0

According to Lemma 2.4, we obtain that A1 )22 . x3 (ω)∗ ≤ (−R21 + ε)/(  Then, it follows from the arbitrariness of ε and −R21 < 0 that x3 (ω)∗ = 0, which ∈ / aa13 , a11 . In other words, x1 ∗ > 0 contradicts with x3 (ω)∗ > 0. Therefore, rr11 31 33 a31 a.s. Case (iii) x2 (ω)∗ > 0 and x3 (ω)∗ > 0. Noting from (11) and H (C0 )∗ < μ1 , we need consider, respectively, the following three situations. (a) If μ2 < H (C0 )∗ < μ1 , then it follows from (15) that limt→∞ x2 (t, ω) = 0, / (μ2 , μ1 ). which contradicts with x2 (ω)∗ > 0. Therefore, H (C0 )∗ ∈ (b) If μ3 < H (C0 )∗ ≤ μ2 , then it follows from (18) that limt→∞ x3 (t, ω) = 0, / (μ3 , μ2 ]. which contradicts with x3 (ω)∗ > 0. Therefore, H (C0 )∗ ∈ (c) If H (C0 )∗ ≤ μ3 , we first estimate x2 , x3 . Multiplying both sides of (12), (13) and (16), respectively, by m, ¯ −1 and n, ¯ then adding these three equalities, we get that

123

Y. Zhao et al.

m¯ x1 (t, ω) n¯ x3 (t, ω) A2 −  1 x2 (t, ω) A2 H (C0 (t)) ln = ln + ln + t x2 (0) t x1 (0) t x3 (0) A22 A − x2 (t, ω) + m¯ f 1 (t) − f 2 (t) + n¯ f 3 (t) A22 t t t A12 0 σ1 (s)dB1 (s) σ2 (s)dB2 (s) A32 0 σ3 (s)dB3 (s) − 0 + , + A22 t t A22 t (24)  a11 m¯ + a31 n¯ = a21 , where m, ¯ n¯ satisfy the equations > which implies m¯ = AA12 22 a13 m¯ + a33 n¯ = a23 , > 0. Noting that A2 −μ2 = A32 R > 0, we have that if H (C0 )∗ ≤ μ3 , 0, n¯ = AA32 22 then H (C0 )∗
0, there is a T such that when t > T, H (C0 (t)) > H (C0 )∗ − ε and 4 A2

max

 m¯ t

ln

 ε x1 (t) n¯ x3 (t) . , ln , |m¯ f 1 (t) − f 2 (t) + n¯ f 3 (t)| < x1 (0) t x3 (0) 4 A22

It follows from (24) that when t > T , ln

 t  A2 H (C0 )∗ + ε A A12 t x2 (t, ω) A2 −  t− x2 (s, ω)ds + σ1 (s)dB1 (s) ≤ x2 (0) A22 A22 0 A22 0  t  A32 t − σ2 (s)dB2 (s) + σ3 (s)dB3 (s). A22 0 0

According to Lemma 2.4, we obtain that x2 (ω)∗ ≤

A2 H (C0 )∗ + ε A2 −  . A

(25)

Discussing as above, from (17) and Lemmas 2.3–2.4, we can have that if H (C0 )∗ ≤ μ3 , then x3 (ω)∗ ≤

A3 H (C0 )∗ + ε A3 −  . A

(26)

Substituting (25) and (26) into (20), we can easily get that σ12 A2 H (C0 )∗ + ε) a12 (A2 −  − r11 H (C0 )∗ − 2 A  a13 (A3 − A3 H (C0 )∗ + ε) − A A1 H (C0 )∗ ) (a13 + a12 ) a11 (A1 −  = − ε. (27) A A

a11 x1 (ω)∗ ≥ r1e −

123

Survival and Stationary Distribution Analysis of a. . .

Noting that

A1 −μ1  A1

=

         σ2 σ2 σ2 σ2 r1e − 21 r21 −r11 r2e − 22 A21 − r31 r1e − 21 −r11 r3e − 23 A31 r11  A1

> 0 and (11), we have that if H (C0 )∗ ≤ μ3 , then H (C0 )∗ < (ω)∗

A1 .  A1

Since ε is arbi-

trary, it follows from (27) that x1 > 0, which contradicts with x1 (ω)∗ = 0. In ∗ other words, x1  > 0 a.s. Next, we consider species x2 . If H (C0 )∗ < μ2 , then from above analysis, we know that x1 ∗ > 0. By the definition of H (C0 )∗ we have that for any ε > 0, there is a T such that when t > T, H (C0 (t)) > H (C0 )∗ − 2rε11 , f 1 (t) ≤ 2ε . It follows from (12) that    σ12 1 x1 (t) 1 t ln ≤ r1e − σ1 (s)dB1 (s). − r11 H (C0 )∗ + ε − a11 x1 (t) + t x1 (0) 2 t 0 (28) According to Lemma 2.4, we obtain that ∗

x1  ≤

r1e −

σ12 2

− r11 H (C0 )∗ + ε . a11

(29)

Then, using Lemma 2.3, it follows from (13) and (29) that 

 σ2 1 x2 (t) ln ≥ r2e − 2 − r21 H (C0 )∗ − 3j=1 a2 j x j ∗ 0 ≥ lim sup t x2 (0) 2 t→∞ σ2

r1e − 21 − r11 H (C0 )∗ + ε σ2 ≥ r2e − 2 − r21 H (C0 )∗ − a21 2 a11 − a22 x2 ∗ − a23 x3 ∗ A2 )33 H (C0 )∗ (A2 )33 − (  a21 = − ε − a22 x2 ∗ − a23 x3 ∗ . a11 a11 That is, a22 x2 ∗ + a23 x3 ∗ ≥

A2 )33 H (C0 )∗ (A2 )33 − (  a21 − ε. a11 a11

(30)

Since ε is arbitrary, it follows from (30) that if H (C0 )∗ < μ2 , then a22 x2 ∗ + a23 x3 ∗ > 0. We claim that x2 ∗ > 0 a.s. In fact, for any ω ∈ {x2 (ω)∗ = 0}, we must have that x3 (ω)∗ > 0. From (11) and H (C0 )∗ < μ2 , we need consider, respectively, the following two cases. (a) If μ3 < H (C0 )∗ < μ2 , then it follows from (18) that limt→∞ x3 (t, ω) = 0, / (μ3 , μ2 ). which contradicts with x3 (ω)∗ > 0. Therefore, H (C0 )∗ ∈ (b) If H (C0 )∗ ≤ μ3 , we first estimate x1 . Multiplying both sides of (12), (13) and (16) by 1, −m, −n, respectively, and then adding these three equalities, we obtain

123

Y. Zhao et al.

m x2 (t, ω) n x3 (t, ω) A1 −  1 x1 (t, ω) A1 H (C0 (t)) ln = ln + ln + t x1 (0) t x2 (0) t x3 (0) A11 t σ1 (s)dB1 (s) A − x1 (t, ω) + f 1 (t) − m f 2 (t) − n f 3 (t) + 0 A11 t t t A31 0 σ3 (s)dB3 (s) A21 0 σ2 (s)dB2 (s) − , (31) + A11 t A11 t 

a22 m + a32 n = a12 , which implies m = AA21 > 11 a23 m + a33 n = a13 , > 0. Discussing as that for (25), using Lemmas 2.3 and 2.4, it follows 0, n = − AA31 11 from (31) that if H (C0 )∗ ≤ μ3 , then where m, n satisfy the equations

x1 (ω)∗ ≤

A1 H (C0 )∗ + ε A1 −  . A

(32)

Taking upper limit on both sides of (13) and substituting (26) and (32) into it, we can easily get that 

 σ2 1 x2 (t) 0 ≥ lim sup ln ≥ r2e − 2 − r21 H (C0 )∗ − 3j=1 a2 j x j (ω)∗ t x2 (0) 2 t→∞ σ22 A1 H (C0 )∗ + ε) a21 (A1 −  ≥ r2e − − r21 H (C0 )∗ − 2 A  (A − A H (C ) + ε) a 23 3 3 0 ∗ − a22 x2 (ω)∗ − A A2 H (C0 )∗ ) (a21 + a23 ) a22 (A2 −  − ε − a22 x2 (ω)∗ . = A A Since ε is arbitrary, it follows from H (C0 )∗ < x2 (ω)∗ ≥

A2  A2

that

A2 H (C0 (t))∗ A2 −  > 0, A

which contradicts with x2 (ω)∗ = 0. Therefore, H (C0 )∗ ∈ / (0, μ3 ]. In other words, x2 ∗ > 0. Finally, we consider species x3 . If H (C0 )∗ < μ3 , then from above analysis, we know that x1 ∗ > 0 and x2 ∗ > 0. Taking upper limit on both sides of (16) and substituting (25) and (32) into it, using Lemma 2.3, we can get that σ2 a33 x3 ∗ ≥ r3e − 3 − r31 H (C0 )∗ − a31 x1 ∗ − a32 x2 ∗ 2 σ32 A1 H (C0 )∗ + ε) a31 (A1 −  − r31 H (C0 )∗ − ≥ r3e − 2 A A2 H (C0 )∗ + ε) a32 (A2 −  − A

123

Survival and Stationary Distribution Analysis of a. . .

=

A3 H (C0 )∗ ) (a31 + a32 ) a33 (A3 −  − ε. A A

Since ε is arbitrary, it follows that if H (C0 )∗ < μ3 , then x3 ∗ ≥

A3 H (C0 )∗ A3 −  > 0. A



Remark 3.1 From Theorems 3.1–3.3, we know that if H (C0 )∗ > μi , then xi will go to extinction; if H (C0 )∗ = μi , then xi will be nonpersistent in the mean; if H (C0 )∗ < μi , then xi will be weakly persistent in the mean. That is to say, H (C0 )∗ − μi is the threshold between weak persistence in the mean and extinction for the ith species. More precisely, (a) If R, R31 > 0, then μ1 > μ2 > μ3 . If H (C0 )∗ > μ1 , then all species x1 , x2 , x3 will go to extinction; if μ2 < H (C0 )∗ < μ1 , then x2 , x3 will go to extinction and x1 will be weakly persistent in the mean; if μ3 < H (C0 )∗ < μ2 , then x3 will go to extinction and x1 , x2 will be weakly persistent in the mean; and if H (C0 )∗ < μ3 , then x1 , x2 , x3 will be weakly persistent in the mean. From r −

σ12

r −

σ22

r −

σ32

a biological viewpoint, since 1er11 2 > 2er21 2 ≥ 3er31 2 , which implies that x1 owns the largest intrinsic growth rate, the strongest antinoise ability and the smallest dose–response parameter to the organismal toxicant concentration, then it is most possible for x1 to survival, followed by species x2 and then x3 in order. (b) If R > 0, R31 = 0, then μ1 = μ2 > μ3 . This implies that species x1 and x2 own the same antitoxic and antinoise ability; thus, they have the same survival probability and more than species x3 . In the above, we have established the threshold between weak persistence in the mean and extinction for species xi , i = 1, 2, 3. Now, we shall strengthen the conditions to give some other persistence results. Theorem 3.4 If limt→+∞ H (C0 (t)) exists, and suppose that for i = 1, 2, 3 i lim H (C0 (t)) > 0. Ai − A t→+∞

then species xi , i = 1, 2, 3 will all be strongly persistent in the mean almost surely, and moreover

lim xi (t) =

t→+∞

i lim H (C0 (t)) Ai − A t→+∞

A

,

i = 1, 2, 3.

Proof If limt→∞ H (C0 (t)) exists, then for any  > 0, there exists a T such that when t > T, H (C0 (t)) > limt→+∞ H (C0 (t)) − A11 . Moreover, from Lemma 2.3, 4 A1

we know that for the  given above, there exists a T1 (> T ) such that when t >

123

Y. Zhao et al. 2 (t) 3 (t) T1 , mt lnx2x(0) < 4 , nt lnx3x(0) < (31) that when t > T1 ,

 4

and f 1 (t) − m f 2 (t) − n f 3 (t)
0 is arbitrary, according to Lemma 2.4, we can have that ∗

x1  ≤

A1 lim H (C0 (t)) A1 −  t→+∞

A

.

(33)

.

(34)

.

(35)

Similarly, from (24) and (17), we can have that ∗

x2  ≤

A2 lim H (C0 (t)) A2 −  t→+∞

A

and ∗

x3  ≤

A3 −  A3 lim H (C0 (t)) t→+∞

A

Then, it follows from (33)–(35) and the existence of limt→+∞ H (C0 (t)) that for any  > 0, there exists a T¯ such that when t > T¯ , a1i xi (t) < a1i

i lim H (C0 (t)) Ai − A t→+∞

A

 + , 4

i = 1, 2, 3.

and  H (C0 (t)) < lim H (C0 (t)) + . t→+∞ 4 It follows from (12) that when t > T ,   σ12 1 x1 (t) ln > r1e − − r11 lim H (C0 (t)) t→+∞ t x1 (0) 2    A1 − A1 lim H (C0 (t)) t→+∞ −  − a11 A    A2 − A2 lim H (C0 (t)) t→+∞ − a12 A

123

Survival and Stationary Distribution Analysis of a. . .



 A3 −  A3 lim H (C0 (t)) t→+∞

− a13 A t σ1 (s)dB1 (s) + 0 t t σ1 (s)dB1 (s) , = − + 0 t 1 t→∞ t

which implies that lim inf

ln x1 (t) ≥ −. Then, we can derive that lim

1 ln x1 (t) t→∞ t 0 and limt→+∞ 1t ln x3 (t) =

0. Similarly, we can derive that limt→+∞ 1t ln x2 (t) = Taking limits on both sides of (31), (24) and (17), we get i lim H (C0 (t)) Ai − A t→+∞

lim xi (t) =

A

t→+∞

> 0,

i = 1, 2, 3.

= 0.



Remark 3.2 Theorem 3.4 tells us that if we know the ultimate concentration of toxicant in the mean, then we could obtain the ultimate scale of each species, which provided a methods to estimate the fluctuation range of each species under toxicant and noise disturbance. On the other hand, permanence indicates that if all species are initially present, even in a low abundances, their abundances reach and remain henceforth over a sizeable threshold (Levin et al. 1984). For stochastic system (8), we can further study its stochastic permanence. Theorem 3.5 If min

i=1,2,3

rie − ri1 lim supt→+∞ H (C0 (t))

!

>

3 max 2 i=1,2,3

! σi2 , then

model (8) will be stochastically permanent. Proof The proof is inspired by the method of Liu et al. (2011). We shall divide the whole proof into two parts. Firstly, we claim that for any arbitrary  > 0, there exists a constant γ > 0 such 3 that P∗ {|x(t)| ≥ γ } ≥ 1 − . Define V1 (x) = ( i=1 xi )−1 , for x ∈ R3+ , then Itô’s formula implies that

dV1 (x) =

−V12 (x)

3 

dxi +

i=1

V13

 3 

2 dxi

i=1

⎧ ⎡ ⎛ ⎞⎤ 3 3 ⎨   = −V12 (x)⎣ xi ⎝rie + (ri0 (0) − rie ) e−αi t − ri1 H (C0 (t)) − aij x j⎠⎦ ⎩ i=1 j=1 * 3 * 3 + +⎫ ⎬   + V13 (x) σi2 (t)xi2 σi (t)xi d Bi (t) . dt − V12 (x) ⎭ i=1

i=1

123

Y. Zhao et al.

Since min

i=1,2,3

! rie − ri1 lim supt→+∞ H (C0 (t)) >

positive constant θ such that

3 max 2 i=1,2,3

! σi2 , we can choose a

   θ  3 ! rie − ri1 lim sup H (C0 (t)) − max σi2 > max σi2 . i=1,2,3 2 i=1,2,3 2 i=1,2,3 t→+∞ min

Define V2 (x) = (1 + V1 (x))θ . Making use of Itô’s formula again leads to θ (θ − 1) (1 + V1 (x))θ −2 d(V1 (x))2 2 * 3   2 xi rie + (ri0 (0) − rie ) e−αi t − (1 + V1 (x)) V1 (x)

dV2 (x) = θ (1 + V1 (x))θ −1 dV1 (x) + / = θ (1 + V1 (x))θ −2

−ri1 H (C0 (t)) −

+

θ +1 4 V1 (x) 2



3 

aij x j ⎠⎦ + (1 +

j=1 3 

− θ (1 + V1 (x))

*

V1 (x))V13 (x)

σi2 (t)xi2 *

V12 (x)

3 

= θ (1 + V1 (x))

+ σi2 (t)xi2

dt + σi (t)xi dBi (t)

i=1 θ −2

3  i=1

0

i=1 θ −1

i=1

⎞⎤

* θ −1

G(x)dt − θ (1 + V1 (x))

V12 (x)

3 

+ σi (t)xi dBi (t) ,

i=1

(36) where  3 G(x) = −(1 + V1 (x))V12 (x) xi (rie + (ri0 (0) − rie )e−αi t − ri1 H (C0 (t)) i=1    3 3 3 2 2 − aij x j ) + (1 + V1 (x))V1 (x) σ (t)xi j=1 i=1 i   3 θ +1 4 σ 2 (t)xi2 V1 (x) + i=1 i 2  3 * 3    2 −1 xi ) xi rie + (ri0 (0) − rie ) e−αi t ≤ −V1 (x)(1 + i=1

i=1

+

− ri1 lim sup H (C0 (t)) − ε t→+∞

⎛

+

V12 (x) ⎝1 +



3 

−1 ⎞ ⎠ max xi

i=1

+ (1 + V1 (x))V1 (x) max

i=1,2,3

123

i, j=1,2,3

aij

 ! 3 i=1

2 xi

! θ +1 ! σi2 + max σ 2 V 2 (x) 2 i=1,2,3 i 1

Survival and Stationary Distribution Analysis of a. . .



! θ +3 min {rie − ri1 lim sup H (C0 (t)) − } − max σi2 i=1,2,3 2 i=1,2,3 t→+∞    ! ! −αi t 2 + (ri0 (0) − rie )e +V1 (x) max σi + max aij + max

≤ −V12 (x)

i=1,2,3

i, j=1,2,3

i, j=1,2,3

! aij .

Substituting the above inequalities into (36) results in   dV2 (x) ≤ θ (1 + V1 (x))θ−2 − V12 (x) min {rie − ri1 lim sup H (C0 (t)) − } i=1,2,3 t→+∞   ! ! θ +3 − max σi2 + (ri0 (0) − rie )e−αi t + V1 (x) max σi2 i=1,2,3 2 i=1,2,3  1 ! ! + max aij + max aij dt i, j=1,2,3

i, j=1,2,3

*

− θ (1 + V1 (x))θ−1 V12 (x)

3 

+ σi (t)xi dBi (t) .

i=1

Choose a positive constant κ sufficiently small such that  1 ! κ θ +3 max σi2 > > 0. min rie − ri1 lim sup H (C0 (t)) −  − i=1,2,3 2 i=1,2,3 θ t→+∞ Define V3 = exp{κt}V2 (x(t)), then applying Itô’s formula again results in dV3 (x) = κ exp{κt}V2 (x)dt + exp{κt}dV2 (x)  ≤ exp{κt}(1 + V1 (x))θ−2 κ(1 + V1 (x))2  1 ! θ +3 max σi2 − min rie − ri1 lim sup H (C0 (t)) −  − i=1,2,3 i=1,2,3 2 t→+∞   ! ! + (ri0 (0) − rie )e−αi t + V1 (x)θ max σi2 + max aij 

V12 (x)θ

i=1,2,3

+ θ max

i, j=1,2,3

i, j=1,2,3

* 3 + 1  ! θ−1 2 aij dt−θ exp{κt}(1+V1 (x)) V1 (x) σi (t)xi dBi (t) *

= exp{κt}H (x)dt − θ exp{κt}(1 + V1 (x))

θ−1

V12 (x)

i=1 3 

+

σi (t)xi dBi (t)

i=1

(37)

123

Y. Zhao et al.

for sufficiently large t, where    1 H (x) = θ (1 + V1 (x))θ−2 − V12 (x) min rie − ri1 lim sup H (C0 (t)) −  i=1,2,3

t→+∞

    ! κ ! θ +3 max σi2 − + V1 (x) max σi2 + ri0 (0) − rie e−αi t − i=1,2,3 2 i=1,2,3 θ  1 ! 2κ ! κ + max aij + + max aij + . i, j=1,2,3 i, j=1,2,3 θ θ It is easy to know that H (x) is upper bounded in R3+ , H1  supx∈R3 H (x) < +∞. + Consequently, * dV3 (x(t)) ≤ H1 exp{κt}dt − θ exp{κt}(1 + V1 (x))θ−1 V12 (x)

3 

+ σi (t)xi dBi (t)

i=1

for sufficiently large t. Integrating both sides of the above inequality and taking expectations result in E[V3 (x(t))] = E[exp{κt}(1 + V1 (x))θ ] ≤ (1 + V1 (x(0)))θ +

H1 exp{κt}, κ

which implies that lim supt→+∞ E[V1θ (x(t))] ≤ lim supt→+∞ E[(1 + V1 (x))θ ] ≤ H1 κ . θ

3 x (t))θ ≤ (3 max 2 θ θ Notice that ( i=1 i i=1,2,3 x i (t)) = 3 (maxi=1,2,3 x i (t)) 2 ≤ 3θ |x(t)|θ . Therefore,

 lim sup E t→+∞

For any  > 0, let γ =

1 |x(t)|θ



≤ 3θ

H1  H. κ

1

θ

1

Hθ

, by Chebyshev’s inequality, we get that 

P{|x(t) < γ |} = P

1 1 > θ |x(t)|θ γ

1

≤ γθE

That is P ∗ {|x(t) < γ |} ≤ γ θ H = . Consequently, P∗ {|x(t) ≥ γ |} ≥ 1 − .

123



 1 , |x(t)|θ

Survival and Stationary Distribution Analysis of a. . .

Next, we prove that for any ε > 0, there exists χ > 0 such that P∗ {|x(t) ≤ χ |} ≥ 3 x p (t), 1 − . The proof is similar to Liu et al. (2011), one can define V (x) = i=1 i then it follows from Itô’s formula that dV (x(t)) =

3 

pxi [rie + (ri0 (0) − rie )e−αi t − ri1 H (C0 (t)) p

i=1

−

3  j=1

≤ p

3 

 p−1 2 p σi (t)]dt + pσi (t)xi d Bi (t) 2 3

aij x j +

i=1

xi [rie + (ri0 (0) − rie )e−αi t − aii xi + p

i=1

+ p

3 

p 2 σ ]dt 2 i

p

σi (t)xi d Bi (t)

(38)

i=1

Let η0 besufficiently large such that every component of xi (0) is contained within  the interval η10 , η0 . For each integer η ≥ η0 , defining the stopping time   1     1 1 1 , η or x2 (t) ∈ , η or x3 (t) ∈ ,η . τη = inf t ≥ 0 : x1 (t) ∈ / / / η η η Obviously, τη → ∞ almost surely as η → ∞. Applying Itô’s formula again to exp{t}V (x) and taking expectations on both sides yields E[exp{t ∧ τη }x p (t ∧ τη )] − x p (0)    t∧τη  3 p 2 1 p s −αi s ds ≤E e pxi (s) rie + (ri0 (0) − rie )e − aii xi (s) + σi + 2 p 0 i=1 +  t∧τη * 3  ≤E (ri0 (0) − rie )e−(αi −1)s ds es (L 1 + L 2 + L 3 ) + 0

i=1

3  ri0 (0) − rie 1 − e−(αi −1)t , ≤ (L 1 + L 2 + L 3 ) (exp{t} − 1) + αi − 1 i=1

where L i = L i ( p) = η → ∞ leads to

p2 σ 2

[1+rie p+ 2 i ] p+1 p ( p+1) p+1 aii

(i = 1, 2, 3) is positive constants. Letting

exp{t}E[x p (t)] ≤ x p (0) + (L 1 + L 2 + L 3 )(exp{t} − 1) 3  ri0 (0) − rie 1 − e−(αi −1)t . + αi − 1 i=1

123

Y. Zhao et al.

That is to say lim sup E[x p (t)] ≤ L 1 + L 2 + L 3 .

(39)

t→∞

1

Then, for any ε > 0, set χ =

(L 1 +L 2 +L 3 ) p

we can get that

1

εp

, by means of Chebyshev’s inequality,

! ! E[|x(t)| p ] P |x(t)| > χ = P |x(t)| p > χ p ≤ . χp In other words, we can obtain that lim sup P{|x(t)| > χ } ≤ lim sup t→∞

t→∞

E[|x(t)| p ] ≤ ε, χp 

which is the desired assertion.

Remark 3.3 Theorem 3.5 tells us all the species will be stochastically permanent provided that the intensity of white noise, dose–response parameters and the concentration of toxicant are sufficiently small.

4 Existence of Stationary Distribution In this section, we will prove the existence of stationary distribution of model (8). We make the following basic assumption. Assumption 4.1 The limit of u(t) when t tends to infinity exists, i.e., limt→∞ u(t) = ue. Obviously, we can have the following lemma, whose proof can be found in Doanh et al. (2010). Lemma 4.1 If Assumption 4.1 holds, then for system (6), we have lim C E (t) =

t→∞

ue ≡ C Ee , h

lim C0 (t) =

t→∞

ku e ≡ C0e . h(g + m)

Noting that model (8) is a non-autonomous system, we cannot prove directly the existence of its stationary distribution. For this reason, we first apply aggregation method to transform model (6) into the following fast timescale system /

123

dC0 (t) dt dC E (t) dt

= kC E (t) − (g + m)C0 (t), = −hC E (t) + u(t)

(40)

Survival and Stationary Distribution Analysis of a. . .

and the following slower timescale system ⎧ * + 3 ⎪

e  ⎪ ⎪ d x¯1 (t) = x¯1 (t) r1e − r11 H C0 − a1 j x¯ j (t) dt + σ1 x¯1 (t)dB1 (t), ⎪ ⎪ ⎪ j=1 ⎪ ⎪ * + ⎪ ⎨ 3

e  d x¯2 (t) = x¯2 (t) r2e − r21 H C0 − a2 j x¯ j (t) dt + σ2 x¯2 (t)dB2 (t), (41) ⎪ j=1 ⎪ ⎪ * + ⎪ ⎪ 3 ⎪

e  ⎪ ⎪ a3 j x¯ j (t) dt + σ3 x¯3 (t)dB3 (t). ⎪ ⎩ d x¯3 (t) = x¯3 (t) r3e − r31 H C0 − j=1

Since system (41) is an autonomous system, we can use the well-known Hasminskii’s methods to prove the existence of its stationary distribution. Our proof can be outlined as follows. First, we prove that for any arbitrary solution of slower system (41), it is attractive to any solution of model (8). Then, by means of Hasminskii’s methods, we prove the existence of stationary distribution of slower system (41). Since all the solution of model (8) will tend to the stationary distribution of system (41), accordingly they have the same stationary distribution. Let us first make some preparations. Definition 4.1 Let X (t) be an arbitrary solution of system (8) and X¯ (t) be an arbitrary solution of system (41) with the same initial value (9). If lim |xi (t) − x¯i (t)| = 0,

t→∞

i = 1, 2, 3

for almost all ω ∈ , then we say X¯ (t) is globally attractive. Lemma 4.2 Let X (t) = (x1 (t), x2 (t), x3 (t)) be a solution of system (8) with initial value (9), then for all p > 1 and i = 1, 2, 3, p

lim sup E[xi (t)] ≤ L i ( p). t→∞

Proof The proof is similar to that of (39) in Theorem 3.5, we omit it here. p Lemma 4.2 tells us that there is a T > 0, such that E(xi (t)) ≤ 23 L i ( p) for p all t ≥ T . At the same time, it follows from the continuity of E(xi (t)) that there p exists L˜ i ( p) > 0, i = 1, 2, 3, such that E(xi (t)) ≤ L˜ i ( p) for t ≤ T . Denote p 3 L¯ i ( p) = max{ 2 L i ( p), L˜ i ( p)}, then for all t ≥ 0, E(xi (t)) ≤ L¯ i ( p).  Lemma 4.3 (See Karatzas and Shreve 1991) Suppose that a stochastic process X (t) on t ≥ 0 satisfies the condition E|x(t) − x(s)| ≤ c|t − s|1+ , 0 ≤ s, t < ∞, for some positive constants ,  and c, then there exists a continuous modification  X (t) of X (t), which has the property that for every ζ ∈ (0,  ), and a positive random variable h(ω) such that

123

Y. Zhao et al.

P ω:

! | X (t, ω) −  X (s, ω)| 2 = 1. ≤ ζ −ζ |t − s| 1−2 0 0 such that 3  i, j=1

aij ξi ξ j =

3 

αi2 xi2 ξi2 ≥ M|ξ |2 ,

i=1

123

Y. Zhao et al.

for all (x1 , x2 , x3 ) ∈ U¯ , ξ ∈ R3 , which shows that condition (i) of Assumption 6.1 is also satisfied. By Theorem 6.1, we can draw the conclusion that stochastic system (41) has a stationary distribution μ(·) and it is ergodic.  Combining Theorems 4.1 and 4.2, we can have the following result on the existence of stationary distribution of model (8). Theorem 4.3 Assume that the conditions in Theorem 4.2 hold, then there is a stationary distribution μ(·) for model (8) with initial value (9) and it has ergodic property.

5 Numerical Simulations To illustrate the theoretical results obtained above, we numerically simulate the solution of model (8) using the following Milstein scheme (Higham 2001). Consider the discretization of model (8) for t = 0, t, 2t, . . . , nt: ⎧ k+1



  k k k x1 = x1k + x1k [r1e + r10 (0) − r1e e−α1 (kt) ⎪ ⎪

− r11 H C  0 (kt) − a11 x1 − a12 x2 − a13 x3 ]t ⎪ 2 ⎪ 2 1−e−2α1 (kt)

√ 

 

⎪ σ 2 ⎪ ⎪ σ1 1 − e−2α1 (kt) x1k tξk + 1 x1k ξk2 t − t , ⎪ ⎪  −α (kt) 2



 ⎪ k+1 ⎪ k k k k k 2 ⎨ x2 = x2 + x2 [r2e + r20 (0) − r2e e

− r21 H C  0 (kt) − a21 x1 − a22 x2 − a23 x3 ]t 2

 k√  σ22 1−e−2α2 (kt) k 2 2 ⎪ σ2 1 − e−2α2 (kt) x2 tζk + x2 ζk t − t , ⎪ ⎪ ⎪  −α (kt) 2



 ⎪ k+1 k + x k [r + r (0) − r ⎪ 3 e x = x − r H C (kt) − a31 x1k − a32 x2k − a33 x3k ]t ⎪ 3e 30 3e 31 0 3 3 3 ⎪

 ⎪ 2

2 −2α (kt) ⎪ √  



3 σ 1−e ⎪ 2 ⎩ σ3 1 − e−2α3 (kt) x3k tηk + 3 x3k ηk2 t − t , 2

where ξk , ζk , ηk are the N (0, 1)-distributed independent Gaussian random variables which can be generated numerically by pseudo-random number generators. Let in model (8) or (41) r1e = 0.70, r2e = 0.65, r3e = 0.6, r10 (0) = r20 (0) = r30 (0) = 0.5, α1 = α2 = α3 = 0.3, r11 = r21 = r31 = 1, a11 = 1.2, a12 = 0.45, a13 = 0.5, a21 = 0.5, a22 = 1.1, a23 = 0.5, a31 = 0.5, a32 = 0.5, a33 = 1.1, σ1 = 0.2, σ2 = 0.25, σ3 = 0.3. By straightforward computations, we have that μ1 = 0.680, μ2 = 0.575, μ3 = 0.4740. Moreover, it is easy to check that Assumptions 3.1–3.4 are satisfied. The following example concerns the survival results of each species. To approximately estimate the probability density function (PDF), the frequency histogram corresponding to a sample path is also included in each numerical simulation figures. Example 5.1 The only difference between the following cases is the form of C0 (t). (i) In Fig. 1, we choose H (C0 (t)) = 0.75 + 0.1 sin t. It is easy to compute that H (C0 )∗ = 0.75 > μi , i = 1, 2, 3. In view of Theorem 3.1, we know that all the species will go to extinction, which are consistent with the simulation results as shown in Fig.1. (ii) In Fig. 2, we choose H (C0 (t)) = 0.6+0.1 sin t. We can compute that H (C0 )∗ = 0.6. Hence, we have that H (C0 )∗ < μ1 and H (C0 )∗ > max{μ2 , μ3 }. In view of Theorems 3.1 and 3.3, we know that the species x1 will be weakly persistent in the mean, and species x2 and x3 will go to extinction. Figure2 confirms these results.

123

0.05 0

0

500

1000

t

2

x (t)

0.1 0.05 0

0

500

1000

relative frequency density

x1(t)

0.1

relative frequency density

Survival and Stationary Distribution Analysis of a. . . ×104 600 400 200 0 −0.01 4

600

0

(a)

0

500

t

1000

relative frequency density

3

x (t)

0.05

×10

0.01

0.02

0.03

0.04

x1(t) at time 1000

400 200 0 −0.01

t 0.1

0

0 4

400

×10

0.01

0.02

0.03

0.04

x2(t) at time 1000

200

0 −0.02

(b)

0

0.02

0.04

0.06

x3(t) at time 1000

Fig. 1 a Time evolution of model (8) with initial value (0.1, 0.1, 0.1), using the parameter values in Example 5.1. All the species x1 , x2 , x3 goes to extinction, b frequency histogram of model (8) at time 1000 (Color Figure Online)

(iii) In Fig. 3, we choose H (C0 (t)) = 0.5 + 0.1 sin t. It is easy to compute that H (C0 )∗ = 0.5 < μ1 , μ2 and H (C0 )∗ > μ3 . According to Theorems 3.1 and 3.3, species x1 and x2 will be weakly persistent in the mean, and species x3 will go to extinction. Our simulation supports this result as shown in Fig. 3. (iv) In Fig. 4, we choose H (C0 (t)) = 0.3 + 0.1 sin t. It is easy to compute that the parameters satisfy H (C0 )∗ = 0.3 < μi , i = 1, 2, 3. In view of Theorem 3.3, species x1 , x2 and x3 will all be weakly persistent in the mean, please see Fig. 4. From Figs. 1, 2, 3 and 4, we can see that in a polluted environment, whether species xi will be persistent or go to extinction is determined uniquely by its survival threshold H (C0 )∗ − μi . Noting that H (C0 (t)) used in Example 5.1 have no limits when t tends to infinity since they are only the linear functions of sin t. Therefore, we can not determine whether model (8) has a stationary distribution in each case [Figs. 2, 3 and 4b in Example 5.1 are, respectively, the frequency histograms of a sample path of model (8)]. To further illustrate the existence of stationary distribution for model (8), we consider the following example. C0 (t) and C0 (t) = 0.3 + Example 5.2 Choose H (C0 (t)) = 2+C 0 (t) . compute that limt→∞ H (C0 (t)) = 0.1304 < μi and

0.001 sin t . t

It is easy to

123

0.07

0.06

0

200

400

600

800

1000

0.2

2

x (t)

t

0.1 0

0

200

400

600

800

1000

relative frequency density

x (t) 1

0.08

relative frequency density

Y. Zhao et al. 4

200

100

0 0.05

400

600

800

1000

t

(a)

relative frequency density

x (t) 3

0.05

200

0.07

0.08

0.09

0.1

x1(t) at time 1000

×10 300 200 100 0 −0.02

0

0.02

0.04

0.06

0.08

0.1

0.12

0.1

0.12

x (t) at time 1000

0.1

0

0.06 4

t

0

×10

2

4

400

×10

300 200 100 0 −0.02

0

0.02

0.04

0.06

0.08

x (t) at time 1000

(b)

3

0.12 0.1

0

200

400

600

800

1000

t x (t) 2

0.08

0.07

0.06

0

200

400

600

800

1000

relative frequency density

0.14

1

x (t)

0.16

relative frequency density

Fig. 2 a Time evolution of model (8) with initial value (0.2, 0.2, 0.1), using the parameter values in Example 5.1. The species x1 will be weakly persistent, and species x2 and x3 will go to extinction, b frequency histogram of model (8) at time 1000 (Color Figure Online) ×104 100

50

0 0.1

300

600

t

800

1000

relative frequency density

x3(t)

(a)

400

0.13

0.14

0.15

0.075

0.08

100 0 0.055

0.06

0.065

0.07

x (t) at time 1000

0.05

200

×10

2

4

0.1

0

0.12

x1(t) at time 1000

200

t

0

0.11 4

400

×10

300 200 100

(b)

0 −0.02

0

0.02

0.04

0.06

0.08

0.1

0.12

x3(t) at time 1000

Fig. 3 a Time evolution of model (8) with (x1 (0), x2 (0), x3 (0) = (0.2, 0.2, 0.1), using the parameter values in Example 5.1. The species x1 and x2 will be weakly persistent in the mean, and species x3 will go to extinction, b frequency histogram of model (8) at time 1000 (Color Figure Online)

123

1

x (t)

0.26 0.24 0.22 0.2

0

500

1000

x2(t)

t 0.19 0.18 0.17 0.16 0.15 0

500

1000

relative frequency density

0.28

relative frequency density

Survival and Stationary Distribution Analysis of a. . .

40 20

0 0.18

60

0.08

0

(a)

500

1000

t

relative frequency density

x3(t)

0.09

0.2

0.22

0.24

0.26

x1(t) at time 1000

×104

40 20 0

t 0.1

×104

300

0.14 ×104

0.16

0.18

0.2

x2(t) at time 1000

200 100 0 0.075

(b)

0.08

0.085

0.09

0.095

0.1

x (t) at time 1000 3

Fig. 4 a Time evolution of model (8) with initial value (0.2, 0.15, 0.1), using the parameter values in Example 5.1. All the species x1 , x2 , x3 will be weakly persistent in the mean, b frequency histogram of model (8) at time 1000 (Color Figure Online)

A1 −  A1 lim H (C0 (t)) = 0.27, t→+∞

A2 −  A2 lim H (C0 (t)) = 0.209, t→+∞

A3 −  A3 lim H (C0 (t)) = 0.1443. t→+∞

Conditions in Theorem 3.4 are satisfied, and therefore, all species x1 , x2 , x3 will be strongly persistent in the mean (see Fig. 5a). Furthermore, we can compute that lim x1 (t) = 0.3114,

t→+∞

lim x2 (t) = 0.2411,

t→+∞

lim x3 (t) = 0.1664.

t→+∞

On the other hand, we can compute from (45) that x1e = 0.3080, x2e = 0.2548, x3e = 0.1715. Moreover, if we choose  = 0.5, then we can further verify that

  ri1 H C0e 2  e 0.0451 = δ < min Ci xi + i=1,2,3 2Ci = min{0.3403, 0.1868, 0.0721} = 0.0721,

123

1

x (t)

0.35

0.3

0.25

0

200

400

600

800

1000

2

x (t)

t

0.25

0.2

0

200

400

600

800

1000

t x3(t)

0.2 0.18 0.16

0

200

(a)

400

600

t

800

1000

relative frequency density relative frequency density relative frequency density

Y. Zhao et al.

30

× 104 Frequncy histogram The PDF of x (t) 1

20 10 0

40

0.25

0.3

0.35

0.4

x1(t) at time 1000

× 104

Frequncy histogram The PDF of x2(t)

30 20 10 0 0.2 80

0.25

0.3

0.35

x2(t) at time 1000

× 104

Frequncy histogram The PDF of x3(t)

60 40 20

(b)

0

0.15

0.16

0.17

0.18

0.19

0.2

0.21

0.22

x (t) at time 1000 3

Fig. 5 a The species x1 , x2 , x3 will be strongly persistent in the mean with the parameter values in Example 5.2. b Blue line frequency histogram of model (8) at time 1000; Red line the probability density function of its corresponding stationary distribution simulated by 3000 sample trajectories of system (41) (Color figure online)

which means that condition (44) in Theorem 4.3 or 4.2 is satisfied. Accordingly, model (8) has a stationary distribution (see Fig. 5b). From the blue lines in Fig. 5b, we can see that the frequency histogram of a sample path of model (8) agrees well with the probability density function simulated by 3000 sample paths of system (41). These results indicate that model (8) has the same stationary distribution with system (41) as stated in Theorem 4.3.

6 Discussion From the ecological viewpoint, species persistence in polluted, stochastically fluctuating environments is an very interesting topic. In this paper, assuming that the species’ growth rate is subjected to both the white noise disturbance and the effect of toxicant, we develop and study a stochastic competitive model of three species with a generalized dose–response function in a polluted environment. It is shown that the extinction, non-persistence and the weak persistence of the ith species xi (t) of model (8) are determined by the threshold H (C0 )∗ − μi . If H (C0 )∗ > μi , then the ith species will go to extinction; if H (C0 )∗ = μi , then the ith species will be nonpersistent; if H (C0 )∗ < μi , then the ith species will be weakly persistent (see Theorems 3.1–3.3). Moreover, sufficient conditions for the ith species xi to be strongly persistent and model (8) to be stochastically permanent are also obtained (see Theorems 3.4, 3.5). Using of Hasminskii’s methods and constructing Lyapunov function, sufficient conditions ensuring the existence of stationary distribution for each

123

Survival and Stationary Distribution Analysis of a. . .

1.5

3

1

The threshold of x

The threshold of x

3

1.5

0.5 0 −0.5

1

0.5

0

−0.5 1

−1 1 1

σ

0.5

3

(a)

1

0.5

σ3

σ

2

0

0

(b)

0.5 0.5 0

σ1

0

Fig. 6 The sensitivity analysis of the threshold H (C0 )∗ − μ3 of x3 with respect to a σ1 , σ2 a σ2 , σ3 . Here H (C0 )∗ = 0.4, other parameters take the same values as that in Examples 5.1 and 5.2 (Color Figure Online)

populations are also given (see Theorem 4.3). Numerical simulations are carried out to support our theoretical results. Our results show that the parameters such as intensity of white noise, dose–response rates of toxicant, the mean stress measure of toxicant in organism, interaction coefficients may have significant impacts on the survival for each species. Taking species x3 as an example, we perform some sensitivity analysis of its survival probability to these factors. • Firstly, we find an interesting phenomenon: With the increase of σ1 and σ2 , the threshold H (C0 )∗ −μ3 of x3 decreases (see Fig. 6a), while with the increase of σ3 , the threshold increases (see Fig. 6b). That is, increasing the noise intensities σ1 and σ2 is advantage to the survival of species x3 , and the increase of σ3 is disadvantage to species x3 . In fact, we can see from the expression of the threshold H (C0 )∗ −μi that the white noise of the ith species is advantage to the survival of the jth species (i = j) and disadvantage to itself, and the large amplitude white noise will lead the species to extinction. In addition, the ith species’ critical amplitude of white noise can be get from its expression of threshold by simple computations. From biological viewpoint, each species competes with others for the resources, habitats and others. Since the ith species’s intrinsic growth rate is subjected to the white noise disturbance, then the jth species will get more resources; in other words, if the white noise is not so strong, it is helpful to the coexistence of all species.

123

Y. Zhao et al.

0.01

4

0

3

The threshold of x

The threshold of x3

3.5 3 2.5 2 1.5 1 0.5 0.4

−0.01

−0.02

−0.03

−0.04 0 0.4

0

0.2

r31

(a)

0.2 0

0

r21

r

11

(b)

0.2 0.2

r

11

0.4

Fig. 7 The sensitivity analysis of the threshold H (C0 )∗ − μ3 of x3 with respect to a r11 , r21 , b r21 , r31 . Here H (C0 )∗ = 0.16, other parameters take the same values as that in Examples 5.1 and 5.2 (Color Figure Online)

• Secondly, the dose–response rates of toxicant mean the capacity of toxicantresistance for each species. If the ith species has a stronger antitoxic ability than the jth species’s, then the ith species will get more chances to survival. More precisely, the dose–response rates r31 are harmful to the survival of itself and other dose–response rates r11 and r21 are helpful to the survival of species x3 (see Fig. 7). • Thirdly, the mean stress measure of toxicant in organism H (C0 ) may have great impacts on the persistence ability of each species for model (8). From numerical results in Examples 5.1–5.2, we can see that if the exogenous rate of input of toxicant into the environment is decreasing, the populations exposed in polluted environment will get more chances to be persistent. In summary, increasing environmental pollution is unfavorable for the survival ability of each populations. • Fourthly, we can see from Fig. 8 that the survival threshold of x3 will increase if a31 and/or a32 increase. That is, increasing the values of a31 and/or a32 will be harmful to the survival of species x3 . Therefore, the interaction coefficient plays an vital role in determining persistence or extinction for each species of model (8). On the other hand, noting that model (8) is a non-autonomous system, it is difficult to directly prove the existence of its stationary distribution. In this paper, we developed a new method for the proof. Using aggregation methods, we first transform system (6) into a fast timescale system (40) and a slower timescale system (41). Since system (41) is a autonomous system, we prove the existence of its stationary distribution using Hasminskii’s methods. We further prove that the solution of system (41) is globally

123

Survival and Stationary Distribution Analysis of a. . .

0.15

The threshold of x3

0.1 0.05 0 −0.05 −0.1 −0.15 −0.2 0.8 0.6

a32

0.4 0.2 0

0

0.1

0.2

0.3

0.5

0.4

0.6

a

31

Fig. 8 The sensitivity analysis of the threshold H (C0 )∗ − μ3 of x3 with respect to a31 , a32 . Here H (C0 )∗ = 0.4, other parameters take the same values as that in Examples 5.1 and 5.2 (Color Figure Online)

attractive with respect to arbitrary solution of model (8). Accordingly, model (8) has the same stationary distribution as system (41). The existence of stationary distribution indicates that the population is dynamic and the solution of model (8) will not trend to a steady state of fixed size. Furthermore, the stationary distribution reflects the spreading of the population distribution around steady state of its corresponding slower system. Unfortunately, it is not possible to predict the exact size of each species, and only their probability density can be estimated by its stationary distribution. Our results have interesting biological implications. The persistence conditions provide insights to protect an endangered species or to prevent the spread of invasive species such as blue-green algae. For example, we can appropriately cut down the pollutant outflow u(t) by closuring of egregious polluters or changing the dose–response rates of toxicant to the species by the use of some artificially methods, such as dilution. The model parameters can be estimated using similar statistical methods as in Jiang et al. (2007). These estimates will, for example, help to accurately estimate economic losses in fishery. More realistic model can be developed by further incorporating the randomness of u(t) produced by the randomly varying environment. To model the random changes in the rate of pollutant input into the environment, we can consider the following mean-reverting Ornstein-Uhlenbeck process: du(t) = α(u e − u(t))dt + βdB(t),

(46)

where B(t) is a standard Brownian motion. All parameters appeared in (46) are positive constants. α is the reversion rate, u e is the mean reversion level or long-run equilibrium of the rate of pollutant input into the environment u(t), and β is the intensity of volatility. We leave this for future consideration.

123

Y. Zhao et al. Acknowledgments The authors would like to thank the two referees for their very important and helpful comments and suggestions which have greatly improved the presentation of this paper. Also, the authors would like to thank Professor Daqing Jiang for his valuable suggestions in preparing the manuscript. Research is supported by the National Natural Science Foundation of China (11271260), the Shanghai Leading Academic Discipline Project (XTKX2012), the Hujiang Foundation of China (B14005), the Innovation Program of Shanghai Municipal Education Committee (13ZZ116), the Natural Science Foundation of Ningxia (NZ13212) and the University Scientific Research Project in Ningxia (NGY2013108).

Appendix: Hasminskii’s Method For the readers’ convenience, we briefly introduce the Hasminskii’s methods used to prove the existence of stationary distribution of system (41) in this paper (see Hasminskii 1980). Let X (t) be a homogeneous Markov process defined in the El (which is a l-dimensional Euclidean space) and be described by the following stochastic differential equation: dX (t) = b(X )dt +

k 

fr (X )dBr (t).

(47)

r =1

The diffusion matrix is defined as follows: D(x) = (aij (x)),

aij =

k 

j

fri (x) fr (x).

(48)

r =1

Assumption 6.1 There exists a bounded domain U ∈ El with regular boundary , having the following properties: (i) In the domain U and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix D(x) is bounded away from zero. (ii) If x ∈ El \U , the mean time τ at which a path emerging from x reaches the set U is finite, and supx∈S E x τ < ∞ for every compact subset S ∈ El . Theorem 6.1 If Assumption 6.1 holds, then the Markov process X (t) has a stationary distribution μ(·). Let g(·) be a function integrable with respect to the measure μ. Then Px

lim

T →∞

1 T



T 0



! g(x)μ(dx) = 1,

g(X (t))dt = El

for all x ∈ El . Remark 6.1 The proof of Theorem 6.1 is given in Hasminskii (1980). Exactly, the existence of stationary distribution with density is referred to Theorem 4.1 at p. 119 and Lemma 9.4 at p. 138 in Hasminskii (1980). To validate Assumption 6.1 (i), it is sufficient to prove that F is uniformly elliptical in U , where Fu = b(x) · u x + [tr (D(x)u x x )]/2, that is, there is a positive number M such that

123

Survival and Stationary Distribution Analysis of a. . . k 

aij (x)ξi ξ j ≥ M|ξ |2 , x ∈ U, ξ ∈ R k .

i, j=1

(see Gard 1988, Chapter 3, p. 103 and Rayleigh’s principle in Strang 1988, Chapter 6, p. 349). To verify Assumption 6.1 (ii), it is enough to show that there exists some neighborhood U and a nonnegative C 2 -function V such that for any El \U, L V is negative (see p. 1163 of Zhu and Yin 2007 for details).

References Bandyopadhyay M, Chattopadhyay J (2005) Ratio-dependent predator–prey model: effect of environmental fluctuation and stability. Nonlinearity 18:913–936 Butler GC (1979) Principles of ecotoxicology. Wiley, New York Chen L, Chen J (1993) Nonlinear biological dynamical system. Science Press, Beijing Dong Y, Wang L (1997) The threshold of survival for system of three-competitive in a polluted environment. J Syst Sci Math Sci 17:221–225 Duan L, Lu Q, Yang Z, Chen L (2004) Effects of diffusion on a stage-structured population in a polluted environment. Appl Math Comput 154:347–359 Danga NH, Dub NH, Yin G (2014) Existence of stationary distributions for Kolmogorov systems of competitive type under telegraph noise. J Differ Equa 257:2078–2101 Doanh NN, Rafael BP, Miguel AZ, Pierre A (2010) Competition and species coexistence in a metapopulation model: can fast asymmetric migration reverse the outcome of competition in a homogeneous environment? J Theor Biol 266:256–263 Dubeya B, Hussainb J (2006) Modelling the survival of species dependent on a resource in a polluted environment. Nonlinear Anal Real World Appl 7:187–210 Duffie D (1996) Dynamic asset pricing theory, 2nd edn. Princeton University Press, Princeton Filov VA, Golubev AA, Liublina EI, Tolokontsev NA (1979) Quantitative toxicology. Wiley, Chichester Gard TC (1992) Stochastic models for toxicant-stressed population. Bull Math Biol 54:827–837 Gard TC (1988) Introduction to stochastic differential equation. Marcel Dekker Inc, New York Hallam TG, Svobada LJ, Gard TC (1979) Persistence and extinction in three species Lotka–Volterra competitive systems. Math Biosci 46:117–124 Hallam TG, Clark CE, Jordan GS (1983) Effects of toxicants on populations: a qualitative approach. First order kinetics. J Math Biol 18:25–27 Hallam TG, De Luna JT (1984) Extinction and persistence in models of population–toxicant interactions. Ecol Model 22:13–20 Hallam TG, Ma Z (1986) Persistence in population models with demographic fluctuations. J Math Biol 24:327–339 He J, Wang K (2009) The survival analysis for a population in a polluted environment. Nonlinear Anal Real World Appl 10:1555–1571 Higham DJ (2001) An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Rev 43:525–546 Hasminskii RZ (1980) Stochastic stability of differential equations. In: Mechanicians analysis, monographs and textbooks on mechanics of solids and fluids. Alphen aan den Rijn, Sijthoff and Noordhoff, Netherlands Ji C, Jiang D, Shi N (2011) A note on predator-prey model with modified Leslie–Gower and Holling-type II schemes with stochastic perturbation. J Math Anal Appl 377:435–440 Jiang D, Ji C, Li X, O’Regan D (2012) Analysis of autonomous Lotka–Volterra competition systems with random perturbation. J Math Anal Appl 390:582–595 Jiang D, Zhang B, Wang D, Shi N (2007) Existence, uniqueness, and global attractive of positive solutions and MLE of the parameters to the logistic equation with random perturbation. Sci China Math 50:977– 986 Karatzas I, Shreve SE (1991) Brownian motion and stochastic calculus. Springer, Berlin Levin SA, Kimball KD, McDowell WH, Kimball SF (1984) New perspectives in ecotoxicology. Environ Manag 8:375–442

123

Y. Zhao et al. Liu H, Ma Z (1991) The threshold of survival for system of two species in a polluted environment. J Math Biol 30:49–61 Liu M, Wang K, Wu Q (2011) Survival analysis of stochastic competitive models in a polluted environment and stochastic competitive exclusion principle. Bull Math Biol 73:1969–2012 Liu M, Wang K (2012) Stationary distribution, ergodicity and extinction of a stochastic generalized logistic system. Appl Math Lett 25:1980–1985 Liu M, Wang K (2011) Survival analysis of a stochastic cooperation system in a polluted environment. J Biol Syst 19:183–204 Li DS (2013) The stationary distribution and ergodicity of a stochastic generalized logistic system. Stat Prob Lett 83:580–583 May RM (1973) Stability and complexity in model ecosystems. Princeton University Press, Princeton Ma Z, Hallam TG (1987) Effects of parameter fluctuations on community survival. Math Biosci 86:35–49 Mandal PS, Banerjee M (2012) Stochastic persistence and stationary distribution in a Holling–Tanner type prey–predator model. Phys A Stat Mech Appl 391:1216–1233 Mao X (2007) Stochastic differential equations and application. Horwood Publishing, Chichester Mao X (2011) Stationary distribution of stochastic population systems. Syst Control Lett 60:398–405 Melbourne BA, Hastings A (2008) Extinction risk depends strongly on factors contributing to stochasticity. Nature 3:100–103 Mezquita F, Hernandez R, Rueda J (1999) Ecology and distribution of ostracods in a polluted Mediterranean river. Palaeogeogr Palaeoclimatol Palaeoecol 148:87–103 Pan J, Jin Z, Ma Z (2000) Thresholds of survival for an n-dimensional volterra mutualistic system in a polluted environment. J Math Anal Appl 252:519–531 Rudnicki R, Pichor K (2007) Influence of stochastic perturbation on prey–predator systems. Math Biosci 206:108–119 Rudnicki R (2003) Long-time behaviour of a stochastic prey–predator model. Stoch Process Appl 108:93– 107 Strang G (1988) Linear algebra and its applications. Thomson Learning Inc, London Wang K (2010) Random mathematical biology model. Science Press, Beijing Wang W, Ma Z (1994) Permanence of populations in a polluted environment. Math Biosci 122:235–248 Wu F, Mao X, Chen K (2008) A highly sensitive mean-reverting process in finance and the Euler–Maruyama approximations. J Math Anal Appl 348:540–554 Yang S, Liu P (2010) Strategy of water pollution prevention in Taihu Lake and its effects analysis. J Great Lakes Res 1:150–158 Zhao Z, Chen L, Song X (2009) Extinction and permanence of chemostat model with pulsed input in a polluted environment. Commun Nonlinear Sci Numer Simul 14:1737–1745 Zhu C, Yin G (2007) Asymptotic properties of hybrid diffusion system. SIAM J Control Optim 46:1155– 1179

123