Analysis of a Delayed Stochastic One-Predator Two-Prey Population Model in a Polluted Environment

Mengwei Li, Yuanfu Shao, Yafei Yang

College of Physics, Guilin University of Technology, Guilin, China.

DOI: 10.4236/jamp.2019.710154                

Abstract

This paper is concerned with the dynamics of a delayed stochastic one-predator two-prey population model in a polluted environment. We show that there exists a unique positive solution that is permanent in time average under certain conditions. Moreover, the global attractively of system is studied. Finally, some numerical simulations are given to illustrate the main results.

Keywords

Random Perturbations, Time Delays, Pollution, Permanent in Time Average

1. Introduction

With the rapid development of economy, environmental pollution has gradually become the major social problem today. With a growing number of toxicant and contaminants entering into the ecosystem, the quality of our living environment has declined. Then many species have been extinct, and some of them are on the edge of extinction. Therefore, controlling environment pollution has become a major topic in many countries, which draws researchers to investigate the influence of environment pollution.

In the 1980s, Hallam et al. [1] [2] [3] firstly proposed the deterministic models to study the impact of environment toxicant on the survival of biological population. Their studies have provided useful bases about protecting species for us. However, population system is often affected by environmental noise, and there are many scholars who have studied the dynamics of stochastic models with toxicant [4] [5] [6] [7] .

On the other hand, more realistic models of population interactions should take the effects of time delay into account [8] [9] . Further, in the natural world, it is a common phenomenon that a predator feeds on some competing preys [10] [11] [12] . However, there is little research on the delayed stochastic one-predator two-prey model in a polluted environment. Thus we consider a stochastic delayed one-predator two-prey model with toxicant input in this article.

The rest is organized as follows. In Section 2, we show some notations and introduce a stochastic delayed one-predator two-prey model in polluted environment. In Section 3, we show that the system (SM) has a unique global positive solution. In Section 4, we give the main theorems and their proof. In Section 5, the attractively global system is investigated. In Section 6, we present numerical simulations to illustrate our mathematical findings.

2. The Model and Notations

In this section, we will give some notations on stochastic one-predator-two-prey system. The stochastic predator-prey system in a polluted environment takes the following form:

$\{\begin{array}{l}\text{d}{x}_1\left(t\right)=x_1\left(t\right)\left[r_1-a_1{C}_1\left(t\right)-c_{11}{x}_1\left(t\right)-c_{12}{x}_2\left(t\right)\right]\text{d}t,\\ \text{d}{x}_2\left(t\right)=x_2\left(t\right)\left[-r_2-a_2{C}_2\left(t\right)+c_{21}{x}_1\left(t\right)-c_{22}{x}_2\left(t\right)\right]\text{d}t,\\ \text{d}{C}_1\left(t\right)=\left[k_1{C}_e\left(t\right)-\left(g_1+m_1\right)C_1\left(t\right)\right]\text{d}t,\\ \text{d}{C}_2\left(t\right)=\left[k_1{C}_e\left(t\right)-\left(g_2+m_2\right)C_2\left(t\right)\right]\text{d}t,\\ \text{d}{C}_e\left(t\right)=\left[-h{C}_e\left(t\right)+u\left(t\right)\right]\text{d}t.\end{array}$ (2.1)

The above model does not incorporate the effect of time delay, but for a long time, it has been recognized that delays can have a complex effect on the dynamics of a system [9] [13] . In the same time, the natural growth of many populations is inevitably affected by many random disturbances. Considering the effects of random disturbances, we assume the growth rate of prey and the death rate of predator are perturbed with

$r_i\to r_i+{\beta }_i{B}_i\left(t\right),i=1,2,3.$

where $B_i\left(t\right)\left(i=1,2,3\right)$ is mutually independent one-dimensional standard Brownian motions with $B_i\left(0\right)=0$ and ${\beta }_i>0\left(i=1,2,3\right)$ being the intensities of white noises. Stochastic version corresponding to deterministic system with time delays can be rewritten as:

$\{\begin{array}{l}\text{d}{x}_1\left(t\right)=x_1\left(t\right)\left[r_1-a_1{C}_1\left(t\right)-c_{11}{x}_1\left(t\right)-c_{12}{x}_2\left(t-\tau \right)\right]\text{d}t+{\beta }_1{x}_1\left(t\right)\text{d}{B}_1\left(t\right),\\ \text{d}{x}_2\left(t\right)=x_2\left(t\right)\left[-r_2-a_2{C}_2\left(t\right)+c_{21}{x}_1\left(t-\tau \right)-c_{22}{x}_2\left(t\right)\right]\text{d}t+{\beta }_2{x}_2\left(t\right)\text{d}{B}_2\left(t\right),\\ \text{d}{C}_1\left(t\right)=\left[k_1{C}_e\left(t\right)-\left(g_1+m_1\right)C_1\left(t\right)\right]\text{d}t,\\ \text{d}{C}_2\left(t\right)=\left[k_1{C}_e\left(t\right)-\left(g_2+m_2\right)C_2\left(t\right)\right]\text{d}t,\\ \text{d}{C}_e\left(t\right)=\left[-h{C}_e\left(t\right)+u\left(t\right)\right]\text{d}t.\end{array}$ (2.2)

In past few decades, delay population systems with one predator and two competing preys have received great attention and have been investigated widely. However, as far as the authors concerned, no one has yet explored the predator-prey system with time delays and toxicant inputs in the same time. Therefore, on the basis of article [14] , we establish the following delayed stochastic one-predator two-prey model in a polluted environment:

$\{\begin{array}{l}\text{d}{x}_1\left(t\right)=x_1\left(t\right)\left[r_1-a_1{C}_0\left(t\right)-c_{11}{x}_1\left(t\right)-c_{12}{x}_2\left(t-{\tau }_{12}\right)-c_{13}{x}_3\left(t-{\tau }_{13}\right)\right]\text{d}t+{\beta }_1{x}_1\left(t\right)\text{d}{B}_1\left(t\right),\\ \text{d}{x}_2\left(t\right)=x_2\left(t\right)\left[r_2-a_2{C}_0\left(t\right)-c_{21}{x}_1\left(t-{\tau }_{21}\right)-c_{22}{x}_2\left(t\right)-c_{23}{x}_3\left(t-{\tau }_{23}\right)\right]\text{d}t+{\beta }_2{x}_2\left(t\right)\text{d}{B}_2\left(t\right),\\ \text{d}{x}_3\left(t\right)=x_3\left(t\right)\left[-r_3-a_3{C}_0\left(t\right)+c_{31}{x}_1\left(t-{\tau }_{31}\right)+c_{32}{x}_2\left(t-{\tau }_{32}\right)-c_{33}{x}_3\left(t\right)\right]\text{d}t+{\beta }_3{x}_3\left(t\right)\text{d}{B}_3\left(t\right),\\ \text{d}{C}_0\left(t\right)=\left[k_1{C}_e\left(t\right)+{\rho }_1\theta \beta /k_1-\left(g_1+m_1\right)C_0\left(t\right)\right]\text{d}t,\\ \text{d}{C}_e\left(t\right)=\left[-h{C}_e\left(t\right)+u\left(t\right)\right]\text{d}t.\end{array}$ (2.3)

with initial data

$x_i\left(\theta \right)={\xi }_i\left(\theta \right),\theta \in \left[-\tau ,0\right],\tau =\max \left\{{\tau }_{ij}\right\},i,j=1,2,3.$

where $x_i\left(t\right)$ is the size of the the prey i, $i=1,2$, and $x_3\left(t\right)$ is the size of the predator; $r_i$ is the growth rate of the i the species, $i=1,2$, $r_3$ is the death rate of the predator; $c_{ii}$ is the intra-specific competition rate, $i=1,2,3$. $c_{12}$ and $c_{21}$ stand for the inter-specific competition rates between species 1 and 2, $c_{13}$ and $c_{23}$ stand for the capture rates, $c_{31}$ and $c_{32}$ are the efficiency of food conversion. $C_0\left(t\right)$ and $C_e\left(t\right)$ denote the concentrations of the toxicant in the organism of species and the environment at time t, respectively. $a_i$ stand for dose-response of the prey and predator to the organismal toxicant, and $-g_i$ and $-m_i$ denote the excretion and depuration rates of the toxicant, $i=1,2,3$, respectively. $k_1$ and ${\rho }_1$ represent the absorption of toxicant per unit of mass by the environment and by food, respectively. $\theta$ is the concentration of toxicant in the environment; $\beta$ is the uptake rate of food per unit mass. Parameter h reflects the ability of the environment to clean up toxicant. $u\left(t\right)$ denotes the exogenous rate of toxicant input into the environment and it is supposed to be bounded and $0\le U_1\le u\left(t\right)\le U_2<\infty$. All coefficients mentioned above are positive constants. ${\tau }_{ij}>0$ represents the time delay. $\xi \left(\theta \right)={\left({\xi }_1\left(\theta \right),{\xi }_2\left(\theta \right),{\xi }_3\left(\theta \right)\right)}^{{\rm T}}\in U$, where $U$ represents the space of all the continued functions from $\left[-\tau ,0\right]$ to $R_{+}^3=\left\{x=\left(x_1,x_2,x_3\right)\in R^3|x_i>0,i=1,2,3\right\}$.

Although the model is a five-dimensional system, because the explicit solutions of the latter two equations are easy to get, it is actually only necessary to study the first three stochastic differential equations of the model, which is called model $(SM)$ in this paper.

For the sake of simplification, we define some notations:

$b_i=r_i-\frac{{\beta }_i^2}{2},i=1,2,b_3=r_3+\frac{{\beta }_3^2}{2},$

$d_i=b_i-a_i{C}_0,i=1,2,3,\Delta =c_{22}{c}_{31}{b}_1+c_{11}{c}_{32}{b}_2-c_{11}{c}_{22}{b}_3,$

$\Gamma =\left|\begin{array}{ccc}{c}_{11}& r_1& {\beta }_1^2/2\\ c_{21}& r_2& {\beta }_2^2/2\\ -c_{31}& -r_3& {\beta }_3^2/2\end{array}\right|,C=\left|\begin{array}{ccc}{c}_{11}& c_{12}& c_{13}\\ c_{21}& c_{22}& c_{23}\\ -c_{31}& -c_{32}& c_{33}\end{array}\right|,$

$C_1=\left|\begin{array}{ccc}{r}_1& c_{12}& c_{13}\\ r_2& c_{22}& c_{23}\\ -r_3& -c_{32}& c_{33}\end{array}\right|,C_2=\left|\begin{array}{ccc}{c}_{11}& r_1& c_{13}\\ c_{21}& r_2& c_{23}\\ -c_{31}& -r_3& c_{33}\end{array}\right|,C_3=\left|\begin{array}{ccc}{c}_{11}& c_{12}& r_1\\ c_{21}& c_{22}& r_2\\ -c_{31}& -c_{32}& -r_3\end{array}\right|,$

${\tilde{C}}_1=\left|\begin{array}{ccc}{\beta }_1^2/2& c_{12}& c_{13}\\ {\beta }_2^2/2& c_{22}& c_{23}\\ {\beta }_3^2/2& -c_{32}& c_{33}\end{array}\right|,{\tilde{C}}_2=\left|\begin{array}{ccc}{c}_{11}& {\beta }_1^2/2& c_{13}\\ c_{21}& {\beta }_2^2/2& c_{23}\\ -c_{31}& {\beta }_3^2/2& c_{33}\end{array}\right|,{\tilde{C}}_3=\left|\begin{array}{ccc}{c}_{11}& c_{12}& {\beta }_1^2/2\\ c_{21}& c_{22}& {\beta }_2^2/2\\ -c_{31}& -c_{32}& {\beta }_3^2/2\end{array}\right|.$

For a function g, we denote the following notations:

$\overline{g\left(t\right)}=\frac{1}{t}{\displaystyle {\int }_0^{t}g\left(s\right)\text{d}s},{\overline{g\left(t\right)}}^{*}=\underset{t\to \infty }{\lim \sup }\frac{1}{t}{\displaystyle {\int }_0^{t}g\left(s\right)\text{d}s},{\overline{g\left(t\right)}}_{\ast }=\underset{t\to \infty }{\lim \inf }\frac{1}{t}{\displaystyle {\int }_0^{t}g\left(s\right)\text{d}s}.$

3. Existence and Uniqueness of the Global Positive Solution

In order to make the model be sense, we need to show the solution is non-negative and global.

Lemma 3.1 ( [15] ) For model (2.3), if $0<k_1+{\rho }_1\theta \beta /k_1<g_1+m_1$, $U_2\le h$, then $0\le C_0\left(t\right)\le 1$, $0\le C_e\left(t\right)\le 1$.

This paper assumes that condition $0<k_1+{\rho }_1\theta \beta /k_1<g_1+m_1$, $U_2\le h$ is always true in model (2.3), then the solution process of model (SM) should be non-negative.

Lemma 3.2 For any initial value $\xi \left(\theta \right)={\left({\xi }_1\left(\theta \right),{\xi }_2\left(\theta \right),{\xi }_3\left(\theta \right)\right)}^{{\rm T}}\in R_{+}^3$, there is a unique global positive solution $\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)\right)\in R_{+}^3,a.s.$. Moreover, there is a positive constants k such that

$\underset{t\to \infty }{\lim \sup }\text{E}\left(x_i\left(t\right)\right)\le k,i=1,2,3.$ (3.1)

Proof. Consider the following system:

$\{\begin{array}{l}\text{d}{N}_1\left(t\right)=\left[r_1-a_1{C}_0\left(t\right)-c_{11}{\text{e}}^{N_1\left(t\right)}-c_{12}{\text{e}}^{N_2\left(t-{\tau }_{12}\right)}-c_{13}{\text{e}}^{N_3\left(t-{\tau }_{13}\right)}\right]\text{d}t+{\beta }_1\text{d}{B}_1\left(t\right),\\ \text{d}{N}_2\left(t\right)=\left[r_2-a_2{C}_0\left(t\right)-c_{21}{\text{e}}^{N_1\left(t-{\tau }_{21}\right)}-c_{22}{\text{e}}^{N_2\left(t\right)}-c_{23}{\text{e}}^{N_1\left(t-{\tau }_{23}\right)}\right]\text{d}t+{\beta }_2\text{d}{B}_2\left(t\right),\\ \text{d}{N}_3\left(t\right)=\left[-r_3-a_3{C}_0\left(t\right)+c_{31}{\text{e}}^{N_1\left(t-{\tau }_{31}\right)}+c_{32}{\text{e}}^{N_2\left(t-{\tau }_{32}\right)}-c_{33}{\text{e}}^{N_3\left(t\right)}\right]\text{d}t+{\beta }_3\text{d}{B}_3\left(t\right).\end{array}$ (3.2)

with initial value $N_i\left(\theta \right)=\log {\xi }_i\left(\theta \right),i=1,2,3$. Since the coefficients of (3.2) obey the local Lipstchiz condition, then (3.2) has a unique local positive solution $N\left(t\right)$ on $\left[0,{\tau }_e\right)$, where ${\tau }_e$ stands for the explosion time. Hence it follows from Itô’s formula that (SM) has the following unique positive global solution

$x\left(t\right)={\left(x_1\left(t\right)={\text{e}}^{N_1\left(t\right)},x_2\left(t\right)={\text{e}}^{N_2\left(t\right)},x_3\left(t\right)={\text{e}}^{N_3\left(t\right)}\right)}^{{\rm T}}.$

Now we show that $x\left(t\right)$ is global, i.e., ${\tau }_e=\infty$. Consider the following system:

$\{\begin{array}{l}\text{d}{y}_1\left(t\right)=y_1\left(t\right)\left[r_1-a_1{C}_0\left(t\right)-c_{11}{y}_1\left(t\right)\right]\text{d}t+{\beta }_1{y}_1\left(t\right)\text{d}{B}_1\left(t\right),\\ \text{d}{y}_2\left(t\right)=y_2\left(t\right)\left[r_2-a_2{C}_0\left(t\right)-c_{22}{y}_2\left(t\right)\right]\text{d}t+{\beta }_2{y}_2\left(t\right)\text{d}{B}_2\left(t\right),\\ \text{d}{y}_3\left(t\right)=y_3\left(t\right)\left[-r_3-a_3{C}_0\left(t\right)+c_{31}{y}_1\left(t-{\tau }_{31}\right)+c_{32}{y}_3\left(t-{\tau }_{32}\right)-c_{33}{y}_3\left(t\right)\right]\text{d}t+{\beta }_3{y}_3\left(t\right)\text{d}{B}_3\left(t\right).\end{array}$ (3.3)

with initial value $y_i\left(\theta \right)={\xi }_i\left(\theta \right),i=1,2,3$. By the stochastic comparison theorem [16] , one can see that for $t\in \left[0,{\tau }_e\right)$,

$x_i\left(t\right)\le y_i\left(t\right),a.s.,i=1,2,3.$

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

$\{\begin{array}{l}{y}_1\left(t\right)=\frac{\exp \left\{d_{1}t+{\beta }_1{B}_1\left(t\right)\right\}}{y_1^{-1}\left(0\right)+c_{11}{\displaystyle {\int }_0^t\exp \left\{d_{1}s+{\beta }_1{B}_1\left(s\right)\right\}\text{d}s}},\\ y_2\left(t\right)=\frac{\exp \left\{d_{2}t+{\beta }_2{B}_2\left(t\right)\right\}}{y_2^{-1}\left(0\right)+c_{22}{\displaystyle {\int }_0^t\exp \left\{d_{2}s+{\beta }_2{B}_2\left(s\right)\right\}\text{d}s}},\\ y_3\left(t\right)=\frac{\exp \left\{-d_{3}t+{\displaystyle {\int }_0^t\left(c_{31}{y}_1\left(s-{\tau }_{31}\right)+c_{32}{y}_2\left(s-{\tau }_{32}\right)\right)\text{d}s}+{\beta }_3{B}_3\left(t\right)\right\}}{y_2^{-1}\left(0\right)+c_{33}{\displaystyle {\int }_0^t\exp \left\{-d_{3}s+{\displaystyle {\int }_0^s\left(c_{31}{y}_1\left(u-{\tau }_{31}\right)+c_{32}{y}_2\left(u-{\tau }_{32}\right)\right)\text{d}u}+{\beta }_3{B}_3\left(s\right)\right\}\text{d}s}}.\end{array}$

Note that $y_1\left(t\right),y_2\left(t\right)$ and $y_3\left(t\right)$ are existent on $t\ge 0$, hence ${\tau }_e=+\infty$.

Before we state the main theorem of this paper, we need to introduce several hypotheses.

Hypothesis 1. $C>0,C_i>0,i=1,2,3$. which imply that all the populations coexist if model (SM) frees from stochastic noises.

Hypothesis 2. $c_{11}>c_{12}+c_{13},c_{22}>c_{21}+c_{23},c_{33}>c_{31}+c_{32}$.

4. Permanence in Time Average

In this section, we study the permanent in time average of systems (2.3) and (SM). We firstly do some preparation.

Definition 4.1. System (2.3) is said to be permanent in time average if there are positive constants $s_i$ and $v_i$ ( $i=1,2,3$ ) such that

$v_i\le \underset{t\to \infty }{\lim \inf }{\bar{x}}_i\left(t\right)\le \underset{t\to \infty }{\lim \sup }{\bar{x}}_i\left(t\right)\le s_i,i=1,2,3.$

holds for any solution $\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)\right)$ of system (2.3) with initial condition $\xi \left(t\right)=\left\{\left({\xi }_1\left(t\right),{\xi }_2\left(t\right),{\xi }_3\left(t\right)\right):-\tau \le t\le 0\right\}\in C\left(\left[-\tau ,0\right]:R_{+}^2\right).$

Lemma 4.1 ( [18] ). Suppose $x\left(t\right)$ be a continuous function from $\Omega \times \left[0,+\infty \right)$ to $R_{+}$.

1) If there exist constants $\lambda ,{\lambda }_0>0$ and $T>0$ such that

$\log x\left(t\right)\le \lambda t-{\lambda }_0{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\beta }_i}{B}_i\left(t\right),$

for $t\ge T$, where $B_i\left(t\right)$ are independent standard Brownian motions and ${\beta }_i$ are constants, $1\le i\le n$, then we have:

$\{\begin{array}{l}{\bar{x}}^{*}\le \lambda /{\lambda }_0,a.s.,\text{if}\lambda \ge 0,\\ \underset{t\to \infty }{\lim }x\left(t\right)=0,a.s.,\text{if}\lambda <0.\end{array}$

2) If there exist positive constants $\lambda ,{\lambda }_0$ and T such that

$\log x\left(t\right)\ge \lambda t-{\lambda }_0{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\beta }_i}{B}_i\left(t\right),$

for all $t\ge T$, where $B_i\left(t\right)$ are independent standard Brownian motions and ${\beta }_i$ are constants, $1\le i\le n$, then if $\lambda \ge 0$, we have ${\bar{x}}_{*}\ge \lambda /{\lambda }_0,a.s.$.

Lemma 4.2 For arbitrary $\tau \ge 0,$ we have

$\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_{t-\tau }^t{y}_i\left(s\right)\text{d}s}=0,a.s.,i=1,2.$ (4.1)

Proof. Consider the first two equations in model (3.3). By [19] , we have

$\{\begin{array}{l}\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t{y}_i\left(s\right)\text{d}s}=\frac{a_i}{c_{ii}}a.s.,\text{if}{b}_i\ge \frac{{\beta }_i^2}{2},i=1,2,\\ \underset{t\to \infty }{\lim }{y}_i\left(t\right)=0a.s.,\text{if}{b}_i<\frac{{\beta }_i^2}{2},i=1,2.\end{array}$

Consequently, if $b_i<{\beta }_i^2/2$, then

$\underset{t\to +\infty }{\lim }\frac{1}{t}{\displaystyle {\int }_{t-\tau }^t{y}_i}\left(s\right)\text{d}s=\underset{t\to +\infty }{\lim }\frac{1}{t}\left({\displaystyle {\int }_0^t{y}_i}\left(s\right)\text{d}s-{\displaystyle {\int }_0^{t-\tau }{y}_i}\left(s\right)\text{d}s\right)=0.$

If $b_i\ge {\beta }_i^2/2$, then

$\underset{t\to +\infty }{\lim }\frac{1}{t}{\displaystyle {\int }_{t-\tau }^t{y}_i}\left(s\right)\text{d}s=\underset{t\to +\infty }{\lim }\frac{1}{t}\left({\displaystyle {\int }_0^t{y}_i}\left(s\right)\text{d}s-{\displaystyle {\int }_0^{t-\tau }{y}_i}\left(s\right)\text{d}s\right)=\frac{a_i}{c_{ii}}-\frac{a_i}{c_{ii}}=0.$

We now consider the following stochastic equation with delays

$\{\begin{array}{l}\text{d}{y}_1\left(t\right)=y_1\left(t\right)\left[\left(r_1-c_{11}{y}_1\left(t\right)\right)\text{d}t+{\beta }_1\text{d}{B}_1\left(t\right)\right],\\ \text{d}{y}_2\left(t\right)=y_2\left(t\right)\left[\left(r_2-c_{22}{y}_2\left(t\right)\right)\text{d}t+{\beta }_2\text{d}{B}_2\left(t\right)\right],\\ \text{d}{y}_3\left(t\right)=y_3\left(t\right)\left[\left(-r_3+c_{31}{y}_1\left(t-{\tau }_{31}\right)+c_{32}{y}_2\left(t-{\tau }_{32}\right)-c_{33}{y}_3\left(t\right)\right)\text{d}t+{\beta }_3\text{d}{B}_3\left(t\right)\right],\\ y_1\left(t\right)={\xi }_1\left(t\right)\in C\left(\left[-\tau ,0\right];R_{+}\right),\\ y_2\left(t\right)={\xi }_2\left(t\right)\in C\left(\left[-\tau ,0\right];R_{+}\right),\\ y_3\left(t\right)={\xi }_3\left(t\right)\in C\left(\left[-\tau ,0\right];R_{+}\right).\end{array}$ (4.2)

Lemma 4.3 Assume that $b_1\ge 0,b_2\ge 0$ and $c_{31}{b}_1/c_{11}+c_{32}{b}_2/c_{22}-b_3\ge 0$, then the solution of (4.2) has the following properties:

$\underset{t\to \infty }{\lim }{\bar{y}}_1\left(t\right)=\frac{b_1}{c_{11}},\underset{t\to \infty }{\lim }{\bar{y}}_2\left(t\right)=\frac{b_2}{c_{22}},\underset{t\to \infty }{\lim }{\bar{y}}_3\left(t\right)=\frac{\Delta }{c_{11}{c}_{22}{c}_{33}},$

$\underset{t\to \infty }{\lim }\frac{\log y_i\left(t\right)}{t}=0,a.s.\left(i=1,2,3\right).$

Proof. The solution of (4.2) has the property that ( [4] )

$\underset{t\to \infty }{\lim }\frac{\log y_1\left(t\right)}{t}=0,\underset{t\to \infty }{\lim }{\bar{y}}_1\left(t\right)=\frac{b_1}{c_{11}},$ (4.3)

By Itô’s formula, we get

$\log y_1\left(t\right)-\log y_1\left(0\right)=b_{1}t-c_{11}{\displaystyle {\int }_0^t{y}_1\left(s\right)\text{d}s}+{\beta }_1{B}_1\left(t\right),$ (4.4)

$\log y_2\left(t\right)-\log y_2\left(0\right)=b_{2}t-c_{22}{\displaystyle {\int }_0^t{y}_2\left(s\right)\text{d}s}+{\beta }_2{B}_2\left(t\right),$ (4.5)

$\begin{array}{l}\log y_3\left(t\right)-\log y_3\left(0\right)\\ =-b_{3}t+c_{31}{\displaystyle {\int }_0^t{y}_1\left(s-{\tau }_{31}\right)\text{d}s}+c_{32}{\displaystyle {\int }_0^t{y}_2\left(s-{\tau }_{32}\right)\text{d}s}-c_{33}{\displaystyle {\int }_0^t{y}_3\left(s\right)\text{d}s}+{\beta }_2{B}_2\left(t\right)\\ =-b_{3}t+c_{31}{\displaystyle {\int }_0^t{y}_1\left(s\right)\text{d}s}+c_{32}{\displaystyle {\int }_0^t{y}_2\left(s\right)\text{d}s}-c_{33}{\displaystyle {\int }_0^t{y}_3\left(s\right)\text{d}s}-c_{31}[{\displaystyle {\int }_{t-{\tau }_{31}}^t{y}_1\left(s\right)\text{d}s}\\ -{\displaystyle {\int }_{-{\tau }_{31}}^0{y}_1\left(s\right)\text{d}s}]-c_{32}\left[{\displaystyle {\int }_{t-{\tau }_{32}}^t{y}_2\left(s\right)\text{d}s}-{\displaystyle {\int }_{-{\tau }_{32}}^t{y}_2\left(s\right)\text{d}s}\right]+{\beta }_3{B}_3\left(t\right).\end{array}$ (4.6)

Dividing both sides of (4.4), (4.5) and (4.6) by t, one can see that

$\frac{1}{t}\log \frac{y_1\left(t\right)}{y_1\left(0\right)}=b_1-c_{11}{\bar{y}}_1\left(t\right)+\frac{{\beta }_1{B}_1\left(t\right)}{t},$ (4.7)

$\frac{1}{t}\log \frac{y_2\left(t\right)}{y_2\left(0\right)}=b_2-c_{22}{\bar{y}}_2\left(t\right)+\frac{{\beta }_2{B}_2\left(t\right)}{t},$ (4.8)

$\begin{array}{l}\frac{1}{t}\log \frac{y_3\left(t\right)}{y_3\left(0\right)}\\ =-b_3+c_{31}{\bar{y}}_1\left(t\right)+c_{32}{\bar{y}}_2\left(t\right)-c_{33}{\bar{y}}_3\left(t\right)-\frac{c_{31}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{31}}^t{y}_1\left(s\right)\text{d}s}-{\displaystyle {\int }_{-{\tau }_{31}}^0{y}_1\left(s\right)\text{d}s}\right]\\ -\frac{c_{32}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{32}}^t{y}_2\left(s\right)\text{d}s}-{\displaystyle {\int }_{-{\tau }_{32}}^0{y}_2\left(s\right)\text{d}s}\right]+\frac{{\beta }_3{B}_3\left(t\right)}{t}.\end{array}$ (4.9)

And it is well known that

$\underset{t\to \infty }{\lim }\frac{B_i\left(t\right)}{t}=0,a.s.,i=1,2,3.$ (4.10)

Utilizing Lemma 4.1, $b_1=r_1-{\beta }_1^2/2\ge 0$ and $b_2=r_2-{\beta }_2^2/2\ge 0$, it is easy to derive that

$\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t{y}_1\left(s\right)\text{d}s}=\frac{b_1}{c_{11}},a.s\mathrm{..}$ (4.11)

And

$\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t{y}_2\left(s\right)\text{d}s}=\frac{b_2}{c_{22}},a.s\mathrm{..}$ (4.12)

Plugging (4.11), (4.12) into (4.7), (4.8), respectively, then combining (4.10) leads to

$\underset{t\to \infty }{\lim }\frac{\log y_1\left(t\right)}{t}=0,a.s\mathrm{..}$ (4.13)

and

$\underset{t\to \infty }{\lim }\frac{\log y_2\left(t\right)}{t}=0,a.s\mathrm{..}$ (4.14)

Besides, computing $\left(4.7\right)\times c_{22}{c}_{31}+\left(4.8\right)\times c_{11}{c}_{32}+\left(4.9\right)\times c_{11}{c}_{22}$, one can derive that

$\begin{array}{l}\frac{c_{22}{c}_{31}}{t}\log \frac{y_1\left(t\right)}{y_1\left(0\right)}+\frac{c_{11}{c}_{32}}{t}\log \frac{y_2\left(t\right)}{y_2\left(0\right)}+\frac{c_{11}{c}_{22}}{t}\log \frac{y_3\left(t\right)}{y_3\left(0\right)}\\ +\frac{c_{11}{c}_{22}{c}_{31}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{31}}^t{y}_1\left(s\right)\text{d}s}-{\displaystyle {\int }_{-{\tau }_{31}}^0{y}_1\left(s\right)\text{d}s}\right]\\ +\frac{c_{11}{c}_{22}{c}_{32}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{32}}^t{y}_2\left(s\right)\text{d}s}-{\displaystyle {\int }_{-{\tau }_{32}}^0{y}_2\left(s\right)\text{d}s}\right]\\ =c_{22}{c}_{31}\left[b_1-c_{11}{\bar{y}}_1\left(t\right)+\frac{{\beta }_1{B}_1\left(t\right)}{t}\right]+c_{11}{c}_{32}\left[b_2-c_{22}{\bar{y}}_2\left(t\right)+\frac{{\beta }_2{B}_2\left(t\right)}{t}\right]\\ +c_{11}{c}_{22}\left[-b_3+c_{31}{\bar{y}}_1\left(t\right)+c_{32}{\bar{y}}_2\left(t\right)-c_{33}{\bar{y}}_3\left(t\right)+\frac{{\beta }_3{B}_3\left(t\right)}{t}\right]\end{array}$

$\begin{array}{l}=c_{22}{c}_{31}{b}_1+c_{11}{c}_{32}{b}_2-c_{11}{c}_{22}{b}_3-c_{11}{c}_{22}{c}_{31}{\bar{y}}_1\left(t\right)-c_{11}{c}_{22}{c}_{32}{\bar{y}}_2\left(t\right)\\ +c_{11}{c}_{22}{c}_{31}{\bar{y}}_1\left(t\right)+c_{11}{c}_{22}{c}_{32}{\bar{y}}_2\left(t\right)-c_{11}{c}_{22}{c}_{33}{\bar{y}}_3\left(t\right)+c_{22}{c}_{31}\frac{{\beta }_1{B}_1\left(t\right)}{t}\\ +c_{11}{c}_{31}\frac{{\beta }_2{B}_2\left(t\right)}{t}+c_{11}{c}_{22}\frac{{\beta }_3{B}_3\left(t\right)}{t}\\ =\Delta -c_{11}{c}_{22}{c}_{33}{\bar{y}}_3\left(t\right)+c_{22}{c}_{31}\frac{{\beta }_1{B}_1\left(t\right)}{t}+c_{11}{c}_{31}\frac{{\beta }_2{B}_2\left(t\right)}{t}+c_{11}{c}_{22}\frac{{\beta }_3{B}_3\left(t\right)}{t}.\end{array}$ (4.15)

In view of (4.1), (4.13), (4.14), (4.15) and Lemma 4.1, we can gain that

$\{\begin{array}{l}\underset{t\to \infty }{\lim }{y}_3\left(t\right)=0,a.s.,\text{if}{c}_{31}{b}_1/c_{11}+c_{32}{b}_2/c_{22}-b_3<0,\\ \underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t{y}_3\left(s\right)\text{d}s}=\frac{\Delta }{c_{11}{c}_{22}{c}_{33}},a.s.,\text{if}{c}_{31}{b}_1/c_{11}+c_{32}{b}_2/c_{22}-b_3\ge 0.\end{array}$ (4.16)

By (3.11)-(3.16), we have

$\underset{t\to \infty }{\lim }\frac{\log y_3\left(t\right)}{t}=0,a.s\mathrm{..}$

This completes the proof.

Lemma 4.4 For arbitrarily $\tau \ge 0$,

$\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_{t-\tau }^t{y}_3\left(s\right)\text{d}s}=0,a.s\mathrm{..}$ (4.17)

Proof. From Lemma 4.3, we can see that either $\underset{t\to \infty }{\lim }{y}_3\left(t\right)=0$ or $\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t{y}_3\left(s\right)\text{d}s}=a$ (a constant).

If $\underset{t\to \infty }{\lim }{y}_3\left(t\right)=0$, then

$\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_{t-\tau }^t{y}_3\left(s\right)\text{d}s}=\underset{t\to \infty }{\lim }\frac{1}{t}\left({\displaystyle {\int }_0^t{y}_3\left(s\right)\text{d}s}-{\displaystyle {\int }_0^{t-\tau }{y}_3\left(s\right)\text{d}s}\right)=0,a.s\mathrm{..}$

If $\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t{y}_3\left(s\right)\text{d}s}=a$, then it is easy to see that

$\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_{t-\tau }^t{y}_3\left(s\right)\text{d}s}=\underset{t\to \infty }{\lim }\frac{1}{t}\left({\displaystyle {\int }_0^t{y}_3\left(s\right)\text{d}s}-{\displaystyle {\int }_0^{t-\tau }{y}_3\left(s\right)\text{d}s}\right)=0,a.s\mathrm{..}$

This completes the proof.

Lemma 4.5 For any initial data $\xi \left(t\right)=\left\{\left({\xi }_1\left(t\right),{\xi }_2\left(t\right),{\xi }_3\left(t\right)\right):-\tau \le t\le 0\right\}\in C\left(\left[-\tau ,0\right];R_{+}^2\right)$, the solution

$\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)\right)$ of system (SM) has the property that

$x_i\left(t\right)\le y_i\left(t\right),a.s.$ for $i=1,2,3$.

where $\left(y_1\left(t\right),y_2\left(t\right),y_3\left(t\right)\right)$ is the solution of system (4.2).

Proof. Let $z_1\left(t\right)=\frac{1}{x_1\left(t\right)}$. Then, by Itô’s formula, we have

$\begin{array}{c}\text{d}{z}_1\left(t\right)=\text{d}\left(\frac{1}{x_1\left(t\right)}\right)=\{\left(-\frac{1}{x_1\left(t\right)}\right)[r_1-c_{11}{x}_1\left(t\right)-c_{12}{x}_2\left(t-{\tau }_{12}\right)\\ -c_{13}{x}_3\left(t-{\tau }_{13}\right)]+\frac{{\beta }_1^2}{x\left(t\right)}\}\text{d}t-\frac{{\beta }_1}{x_1\left(t\right)}\text{d}{B}_1\left(t\right)\\ =-\left[\frac{r_1}{x_1\left(t\right)}-c_{11}-\frac{c_{12}{x}_2\left(t-{\tau }_{12}\right)}{x_1\left(t\right)}-\frac{c_{13}{x}_3\left(t-{\tau }_{13}\right)}{x_1\left(t\right)}\right]\text{d}t+\frac{{\beta }_1^2}{x\left(t\right)}\text{d}t-\frac{{\beta }_1}{x_1\left(t\right)}\text{d}{B}_1\left(t\right)\\ =\left[\left({\beta }_1^2-r_1\right)z_1\left(t\right)+c_{11}\right]\text{d}t-{\beta }_1{z}_1\left(t\right)\text{d}{B}_1\left(t\right)+\left(\frac{c_{12}{x}_2\left(t-{\tau }_{12}\right)}{x_1\left(t\right)}+\frac{c_{13}{x}_3\left(t-{\tau }_{13}\right)}{x_1\left(t\right)}\right)\text{d}t\\ =\left[\left({\beta }_1^2-r_1\right)\text{d}t-{\beta }_1\text{d}{B}_1\left(t\right)\right]z_1\left(t\right)+\left(c_{11}+\frac{c_{12}{x}_2\left(t-{\tau }_{12}\right)}{x_1\left(t\right)}+\frac{c_{13}{x}_3\left(t-{\tau }_{13}\right)}{x_1\left(t\right)}\right)\text{d}t.\end{array}$

This is

$\text{d}{z}_1\left(t\right)=\left[\left({\beta }_1^2-r_1\right)\text{d}t-{\beta }_1\text{d}{B}_1\left(t\right)\right]z_1\left(t\right)+\left(c_{11}+\frac{c_{12}{x}_2\left(t-{\tau }_{12}\right)}{x_1\left(t\right)}+\frac{c_{13}{x}_3\left(t-{\tau }_{13}\right)}{x_1\left(t\right)}\right)\text{d}t.$

Then

$\begin{array}{c}{z}_1\left(t\right)={\text{e}}^{{\displaystyle {\int }_0^t\left(\frac{{\beta }_1^2}{2}-r_1\right)\text{d}s-{\beta }_1\text{d}{B}_1\left(s\right)}}[\frac{1}{x_1\left(0\right)}+{\displaystyle {\int }_0^t\left(c_{11}+\frac{c_{12}{x}_2\left(t-{\tau }_{12}\right)}{x_1\left(t\right)}+\frac{c_{13}{x}_3\left(t-{\tau }_{13}\right)}{x_1\left(t\right)}\right)}{\text{e}}^{{\displaystyle {\int }_0^s\left[r_1-\frac{{\beta }_1^2}{2}\right]\text{d}\tau +{\beta }_1\text{d}{B}_1\left(\tau \right)}}\text{d}s]\\ ={\text{e}}^{\left(\frac{{\beta }_1^2}{2}-r_1\right)t-{\beta }_1\text{d}{B}_1\left(t\right)}\left[\frac{1}{x_1\left(0\right)}+{\displaystyle {\int }_0^t\left(c_{11}+\frac{c_{12}{x}_2\left(t-{\tau }_{12}\right)}{x_1\left(t\right)}+\frac{c_{13}{x}_3\left(t-{\tau }_{13}\right)}{x_1\left(t\right)}\right){\text{e}}^{\left(r_1-\frac{{\beta }_1^2}{2}\right)s+{\beta }_1\text{d}{B}_1\left(s\right)}\text{d}s}\right]\\ \ge {\text{e}}^{\left(\frac{{\beta }_1^2}{2}-r_1\right)t-{\beta }_1\text{d}{B}_1\left(t\right)}\left[\frac{1}{x_1\left(0\right)}+{\displaystyle {\int }_0^t{c}_{11}{\text{e}}^{{}^{\left(r_1-\frac{{\beta }_1^2}{2}\right)s+{\beta }_1\text{d}{B}_1\left(s\right)}}\text{d}s}\right]=y_1^{-1}\left(t\right).\end{array}$

By Lemma 3.2, we obtain that $y_1\left(t\right)$ is the solution of the following equation

$\text{d}{y}_1\left(t\right)=y_1\left(t\right)\left[\left(r_1-c_{11}{y}_1\left(t\right)\right)\text{d}t+{\beta }_1\text{d}{B}_1\left(t\right)\right].$

Hence, we have

$x_1\left(t\right)\le y_1\left(t\right),a.s\mathrm{..}$

In the same way, we can get

$x_2\left(t\right)\le y_2\left(t\right),a.s\mathrm{..}$

On the other hand, let $z_3\left(t\right)=\frac{1}{x_3\left(t\right)}$. Then, by Itô’s formula, we derive that

$\begin{array}{c}\text{d}{z}_3\left(t\right)=\text{d}\left(\frac{1}{x_3\left(t\right)}\right)=-[\left(-\frac{r_3}{x_3\left(t\right)}+\frac{c_{31}{x}_1\left(t-{\tau }_{31}\right)}{x_3\left(t\right)}+\frac{c_{32}{x}_2\left(t-{\tau }_{32}\right)}{x_3\left(t\right)}-c_{33}\right)\text{d}t\\ -\frac{{\beta }_3}{x_3\left(t\right)}\text{d}{B}_3\left(t\right)]+\frac{{\beta }_3^2}{x_3\left(t\right)}\text{d}t\\ =\left[\left({\beta }_3^2+r_3\right)z_3\left(t\right)+c_{33}-c_{31}{x}_1\left(t-{\tau }_{31}\right)z_3\left(t\right)-c_{32}{x}_2\left(t-{\tau }_{32}\right)z_3\left(t\right)\right]\text{d}t\\ +{\beta }_3{z}_3\left(t\right)\text{d}{B}_3\left(t\right)\\ =\left[\left({\beta }_3^2+r_3-c_{31}{x}_1\left(t-{\tau }_{31}\right)-c_{32}{x}_2\left(t-{\tau }_{32}\right)\right)\text{d}t+{\beta }_3\text{d}{B}_3\left(t\right)\right]z_3\left(t\right)+c_{33}\text{d}t.\end{array}$

We then have

$\begin{array}{c}{z}_3\left(t\right)=\frac{1}{x_3\left(0\right)}{\text{e}}^{\left(\frac{{\beta }_3^2}{2}+r_3\right)t+{\beta }_3{B}_3\left(t\right)-c_{31}{\displaystyle {\int }_0^t{x}_1\left(s-{\tau }_{31}\right)\text{d}s-c_{32}{\displaystyle {\int }_0^t{x}_2\left(s-{\tau }_{32}\right)\text{d}s}}}\\ +c_{33}{\displaystyle {\int }_0^t{\text{e}}^{\left(\frac{{\beta }_3^2}{2}+r_3\right)\left(t-s\right)+{\beta }_3\left(B_3\left(t\right)-B_3\left(s\right)\right)-c_{31}{\displaystyle {\int }_s^t{x}_1\left(u-{\tau }_{31}\right)\text{d}u-c_{32}{\displaystyle {\int }_s^t{x}_2\left(u-{\tau }_{32}\right)\text{d}u}}}}\text{d}s\\ \ge \frac{1}{x_3\left(0\right)}{\text{e}}^{\left(\frac{{\beta }_3^2}{2}+r_3\right)t+{\beta }_3{B}_3\left(t\right)-c_{31}{\displaystyle {\int }_0^t{y}_1\left(s-{\tau }_{31}\right)\text{d}s-c_{32}{\displaystyle {\int }_0^t{y}_2\left(s-{\tau }_{32}\right)\text{d}s}}}\\ +c_{33}{\displaystyle {\int }_0^t{\text{e}}^{\left(\frac{{\beta }_3^2}{2}+r_3\right)\left(t-s\right)+{\beta }_3\left(B_3\left(t\right)-B_3\left(s\right)\right)-c_{31}{\displaystyle {\int }_s^t{y}_1\left(u-{\tau }_{31}\right)\text{d}u-c_{32}{\displaystyle {\int }_s^t{y}_2\left(u-{\tau }_{32}\right)\text{d}u}}}}\text{d}s\\ =y_3^{-1}\left(t\right).\end{array}$

Therefore

$x_3\left(t\right)\le y_3\left(t\right),a.s\mathrm{..}$

By Lemma 3.2, we know that $y_3\left(t\right)$ is a solution of the following equation

$\text{d}{y}_3\left(t\right)=y_3\left(t\right)\left[\left(-r_3+c_{31}{y}_1\left(t-{\tau }_{31}\right)+c_{32}{y}_2\left(t-{\tau }_{32}\right)-c_{33}{y}_3\left(t\right)\right)\text{d}t-{\beta }_3\text{d}{B}_3\left(t\right)\right].$

Hence the proof is completed.

Lemma 4.6 The solution $x\left(t\right)$ of (SM) satisfies

$\underset{t\to \infty }{\lim \sup }\log x_i\left(t\right)/t\le 0,a.s.,i=1,2,3.$ (4.18)

Proof. We know that either $\underset{t\to \infty }{\lim }{y}_i\left(t\right)=0,a.s.$ or $\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t{y}_i\left(s\right)\text{d}s}=a$ constant, $i=1,2,3$. If $\underset{t\to \infty }{\lim }{y}_i\left(t\right)=0$, then

$\underset{t\to \infty }{\lim \sup }\log x_i\left(t\right)/t\le \underset{t\to \infty }{\lim \sup }\log y_i\left(t\right)/t\le 0,a.s\mathrm{..}$

If $\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t{y}_i\left(s\right)\text{d}s=a}$ constant, similar to the proof of [20] , we have $\underset{t\to \infty }{\lim }\log y_i\left(t\right)/t=0,a.s.$.

Hence $\underset{t\to \infty }{\lim \sup }\log x_i\left(t\right)/t\le \underset{t\to \infty }{\lim \sup }\log y_i\left(t\right)/t=0,a.s.$. The proof is completed.

Based on the above discussion, now we show the main results in this sectionm

Theorem 4.1 Assume that $\Gamma >0$ and $C_3/{\tilde{C}}_3>1$, then for any initial data $\xi \left(t\right)=\left\{\left({\xi }_1\left(t\right),{\xi }_2\left(t\right),{\xi }_3\left(t\right)\right):-\tau \le t\le 0\right\}\in C\left(\left[-\tau ,0\right];R_{+}^3\right)$, the solution $\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)\right)$ of (SM) has the following properties:

$\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_o^t{x}_i\left(s\right)\text{d}s}=\frac{C_i-{\tilde{C}}_i}{C},a.s.,i=1,2,3.$

Proof. Let $\left(y_1\left(t\right),y_2\left(t\right),y_3\left(t\right)\right)$ be the solution of system (4.2). By Lemma 4.5, we obtain

$\frac{{\displaystyle {\int }_t^{t+\tau }{x}_1\left(s-\tau \right)\text{d}s}}{t}\le \frac{{\displaystyle {\int }_t^{t+\tau }{y}_1\left(s-\tau \right)\text{d}s}}{t}.$

Since

$\underset{t\to \infty }{\lim }\frac{{\displaystyle {\int }_t^{t+\tau }{y}_1\left(s-\tau \right)\text{d}s}}{t}=0.$

Then

$\underset{t\to \infty }{\lim }\frac{{\displaystyle {\int }_t^{t+\tau }{x}_1\left(s-\tau \right)\text{d}s}}{t}=0.$ (4.19)

Similarly,

$\underset{t\to \infty }{\lim }\frac{{\displaystyle {\int }_t^{t+\tau }{x}_2\left(s-\tau \right)\text{d}s}}{t}=0,\underset{t\to \infty }{\lim }\frac{{\displaystyle {\int }_t^{t+\tau }{x}_3\left(s-\tau \right)\text{d}s}}{t}=0.$ (4.20)

It can be straightly shown by Lemma 4.3 that

${\left(\frac{\log x_i\left(t\right)}{t}\right)}^{*}\le \underset{t\to \infty }{\lim }\frac{\log y_i\left(t\right)}{t}=0,i=1,2,3.$ (4.21)

Applying Itô’s formula to (SM) results in

$\begin{array}{c}\log \frac{x_1\left(t\right)}{x_1\left(0\right)}=b_{1}t-c_{11}{\displaystyle {\int }_0^t{x}_1\left(s\right)\text{d}s}-c_{12}{\displaystyle {\int }_0^t{x}_2\left(s-{\tau }_{12}\right)\text{d}s}\\ -c_{13}{\displaystyle {\int }_0^t{x}_3\left(s-{\tau }_{13}\right)\text{d}s}+{\beta }_1{B}_1\left(t\right),\end{array}$ (4.22)

$\begin{array}{c}\log \frac{x_2\left(t\right)}{x_2\left(0\right)}=b_{2}t-c_{21}{\displaystyle {\int }_0^t{x}_1\left(s-{\tau }_{21}\right)\text{d}s}-c_{22}{\displaystyle {\int }_0^t{x}_2\left(s\right)\text{d}s}\\ -c_{23}{\displaystyle {\int }_0^t{x}_3\left(s-{\tau }_{23}\right)\text{d}s}+{\beta }_2{B}_2\left(t\right),\end{array}$ (4.23)

$\begin{array}{c}\log \frac{x_3\left(t\right)}{x_3\left(0\right)}=-b_{3}t+c_{31}{\displaystyle {\int }_0^t{x}_1\left(s-{\tau }_{31}\right)\text{d}s}+c_{32}{\displaystyle {\int }_0^t{x}_2\left(s-{\tau }_{32}\right)\text{d}s}\\ -c_{33}{\displaystyle {\int }_0^t{x}_3\left(s\right)\text{d}s}+{\beta }_3{B}_3\left(t\right).\end{array}$ (4.24)

Dividing both sides of (4.22), (4.23) and (4.24) by time t, one can obtain that

$\begin{array}{c}\frac{1}{t}\log \frac{x_1\left(t\right)}{x_1\left(0\right)}=b_1-c_{11}{\bar{x}}_1\left(t\right)-c_{12}{\bar{x}}_2\left(t\right)-c_{13}{\bar{x}}_3\left(t\right)+\frac{{\beta }_1{B}_1\left(t\right)}{t}\\ +\frac{c_{12}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{12}}^t{x}_2\left(s\right)\text{d}s}-{\displaystyle {\int }_{-{\tau }_{12}}^0{x}_2\left(s\right)\text{d}s}\right]\\ +\frac{c_{13}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{13}}^t{x}_3\left(s\right)\text{d}s}-{\displaystyle {\int }_{-{\tau }_{13}}^0{x}_3\left(s\right)\text{d}s}\right],\end{array}$ (4.25)

$\begin{array}{c}\frac{1}{t}\log \frac{x_2\left(t\right)}{x_2\left(0\right)}=b_2-c_{21}{\bar{x}}_1\left(t\right)-c_{22}{\bar{x}}_2\left(t\right)-c_{23}{\bar{x}}_3\left(t\right)+\frac{{\beta }_2{B}_2\left(t\right)}{t}\\ +\frac{c_{21}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{21}}^t{x}_1\left(s\right)\text{d}s-{\displaystyle {\int }_{-{\tau }_{21}}^0{x}_1\left(s\right)\text{d}s}}\right]\\ +\frac{c_{23}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{23}}^t{x}_3\left(s\right)\text{d}s-{\displaystyle {\int }_{-{\tau }_{23}}^0{x}_3\left(s\right)\text{d}s}}\right],\end{array}$ (4.26)

$\begin{array}{c}\frac{1}{t}\log \frac{x_3\left(t\right)}{x_3\left(0\right)}=-b_3+c_{31}{\bar{x}}_1\left(t\right)+c_{32}{\bar{x}}_2\left(t\right)-c_{33}{\bar{x}}_3\left(t\right)+\frac{{\beta }_3{B}_3\left(t\right)}{t}\\ -\frac{c_{31}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{31}}^t{x}_1\left(s\right)\text{d}s-{\displaystyle {\int }_{-{\tau }_{31}}^0{x}_1\left(s\right)\text{d}s}}\right]\\ -\frac{c_{32}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{32}}^t{x}_2\left(s\right)\text{d}s-{\displaystyle {\int }_{-{\tau }_{32}}^0{x}_2\left(s\right)\text{d}s}}\right].\end{array}$ (4.27)

Let m and n satisfy the following system:

$\{\begin{array}{l}{c}_{22}m-c_{32}n=c_{12},\\ c_{23}m+c_{33}n=c_{13}.\end{array}$

Then,

$m=\frac{C_{21}}{C_{11}}>0,n=-\frac{C_{31}}{C_{11}}>0.$

where $C_{ij}$ represents the complement minor of $c_{ij}$ in the determinant $C,i,j=1,2,3$.

After a simple calculation of $\left(4.25\right)\times \left(-1\right)+\left(4.26\right)\times m+\left(4.27\right)\times n$, we get

$\begin{array}{c}\frac{1}{t}\log \frac{x_1\left(t\right)}{x_1\left(0\right)}=\frac{m}{t}\log \frac{x_2\left(t\right)}{x_2\left(0\right)}+\frac{n}{t}\log \frac{x_3\left(t\right)}{x_3\left(0\right)}+\frac{C_1-{\tilde{C}}_1}{C_{11}}-\frac{C}{C_{11}}{\bar{x}}_1\left(t\right)\\ +\frac{c_{12}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{12}}^t{x}_2\left(s\right)\text{d}s-{\displaystyle {\int }_{-{\tau }_{12}}^0{x}_2\left(s\right)\text{d}s}}\right]\\ +\frac{c_{13}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{13}}^t{x}_3\left(s\right)\text{d}s-{\displaystyle {\int }_{-{\tau }_{13}}^0{x}_3\left(s\right)\text{d}s}}\right]\\ -\frac{c_{21}m}{t}\left[{\displaystyle {\int }_{t-{\tau }_{21}}^t{x}_1\left(s\right)\text{d}s-{\displaystyle {\int }_{-{\tau }_{21}}^0{x}_2\left(s\right)\text{d}s}}\right]\\ -\frac{c_{23}m}{t}\left[{\displaystyle {\int }_{t-{\tau }_{23}}^t{x}_3\left(s\right)\text{d}s-{\displaystyle {\int }_{-{\tau }_{23}}^0{x}_3\left(s\right)\text{d}s}}\right]\\ +\frac{c_{31}n}{t}\left[{\displaystyle {\int }_{t-{\tau }_{31}}^t{x}_1\left(s\right)\text{d}s-{\displaystyle {\int }_{-{\tau }_{31}}^0{x}_1\left(s\right)\text{d}s}}\right].\end{array}$ (4.28)

Substituting(4.18), (4.17), (4.1) into (4.28) gives

$\frac{1}{t}\log \frac{x_1\left(t\right)}{x_1\left(0\right)}\le \frac{C_1-{\tilde{C}}_1}{C_{11}}+\epsilon -\frac{C}{C_{11}}{\bar{x}}_1\left(t\right)+\frac{{\beta }_1{B}_1\left(t\right)-m{\beta }_2{B}_2\left(t\right)-n{\beta }_3{B}_3\left(t\right)}{t}.$ (4.29)

For sufficiently large t. By $C_1/{\tilde{C}}_1>2r_1/{\beta }_1^2>1$, Lemma 4.1 and the arbitrariness of $\epsilon$, we have

${\bar{x}}_1{\left(t\right)}^{*}\le \frac{C_1-{\tilde{C}}_1}{C},a.s\mathrm{..}$ (4.30)

Analogously, let $\tilde{m}$ and $\tilde{n}$ satisfy the following system

$\{\begin{array}{l}{c}_{11}\tilde{m}-c_{31}\tilde{n}=c_{21},\\ c_{13}\tilde{m}+c_{33}\tilde{n}=c_{23}.\end{array}$

Then,

$\tilde{m}=\frac{C_{21}}{C_{22}}>0,\tilde{n}=-\frac{C_{32}}{C_{22}}>0.$

Multiplying both sides of (4.25), (4.26) and (4.27) by $\tilde{m}$, (−1) and $\tilde{n}$, respectively, then adding these three equations, we get

${\bar{x}}_2{\left(t\right)}^{*}\le \frac{C_2-{\tilde{C}}_2}{C},a.s\mathrm{..}$ (4.31)

Plugging (4.1), (4.17), (4.30), (4.31) into (4.27) results in

$\frac{1}{t}\log \frac{x_3\left(t\right)}{x_3\left(0\right)}\le \frac{c_{33}\left(C_3-{\tilde{C}}_3\right)}{C}+\epsilon -c_{33}{\bar{x}}_3\left(t\right)+\frac{{\beta }_3{B}_3\left(t\right)}{t}.$ (4.32)

Noting that $C_3/{\tilde{C}}_3>1$. Then an application of Lemma 4.1 to (4.32) yields that

${\bar{x}}_3{\left(t\right)}^{*}\le \frac{C_3-{\tilde{C}}_3}{C},a.s\mathrm{..}$ (4.33)

Substituting (4.1), (4.17), (4.31) and (4.33) into (4.25) elicits that

$\begin{array}{c}\frac{1}{t}\log \frac{x_1\left(t\right)}{x_1\left(0\right)}\ge b_1-\epsilon -c_{11}{\bar{x}}_1(t)-\frac{c_{12}\left(C_2-{\tilde{C}}_2\right)}{C}-\frac{c_{13}\left(C_3-{\tilde{C}}_3\right)}{C}+\frac{{\beta }_1{B}_1\left(t\right)}{t}\\ +\frac{c_{12}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{12}}^t{x}_2\left(s\right)\text{d}s-{\displaystyle {\int }_{-{\tau }_{12}}^0{x}_2\left(s\right)\text{d}s}}\right]\\ +\frac{c_{13}}{t}\left[{\displaystyle {\int }_{t-{\tau }_{13}}^t{x}_3\left(s\right)\text{d}s-{\displaystyle {\int }_{-{\tau }_{13}}^0{x}_3\left(s\right)\text{d}s}}\right]\\ \ge \frac{c_{11}\left(C_1-{\tilde{C}}_1\right)}{C}-2\epsilon -c_{11}{\bar{x}}_1\left(t\right)+\frac{{\beta }_1{B}_1\left(t\right)}{t}.\end{array}$

for sufficiently large t. Hence, we can further get from Lemma 4.1 and the arbitrariness of $\epsilon$ that

${\bar{x}}_1{\left(t\right)}_{*}\ge \frac{C_1-{\tilde{C}}_1}{C},a.s\mathrm{..}$ (4.34)

In the same way, we can show that

${\bar{x}}_2{\left(t\right)}_{*}\ge \frac{C_2-{\tilde{C}}_2}{C},a.s\mathrm{..}$ (4.35)

Substituting (4.34) and (4.35) into (4.27), and then using (4.1) and (4.17), Lemma 4.1 and the arbitrariness of $\epsilon$, we get

${\bar{x}}_3{\left(t\right)}_{*}\ge \frac{C_3-{\tilde{C}}_3}{C},a.s\mathrm{..}$

This together with (4.30), (4.31), (4.33), (4.34) and (4.35), yields

$\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t{x}_i\left(s\right)\text{d}s}=\frac{C_i-{\tilde{C}}_i}{C},a.s.,i=1,2,3.$

This completes the proof.

As to model (2.3), by utilizing similar techniques to those employed in the proof of Theorem 4.1, one can obtain the similar result.

5. Global Attractivity

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

Proof. Let $g_i$ be the cofactor of the ith diagonal element of $L_C$, where

$L_C=\left(\begin{array}{ccc}{c}_{12}+c_{13}& -c_{12}& -c_{13}\\ -c_{21}& c_{21}+c_{23}& -c_{23}\\ -c_{31}& -c_{32}& c_{31}+c_{32}\end{array}\right).$

Making use of Kirchhoff’s Matrix Tree Theorem (See, e.g., [21] ), one has $g_i>0$.

Define

$V\left(t\right)={\displaystyle \underset{i=1}{\overset{3}{\sum }}{g}_i}\left[\left|\log x_i\left(t\right)-\log {\tilde{x}}_i\left(t\right)\right|+{\displaystyle \underset{j=1,j\ne i}{\overset{3}{\sum }}{c}_{ij}{\displaystyle {\int }_{t-{\tau }_{ij}}^t\left|x_j\left(s\right)-{\tilde{x}}_j\left(s\right)\right|\text{d}s}}\right].$

Reckoning the right differential $\text{d}V\left(t\right)$ deduces that

$\begin{array}{c}\text{d}V\left(t\right)={\displaystyle \underset{i=1}{\overset{3}{\sum }}{g}_i\mathrm{sgn}\left(x_i\left(t\right)-{\tilde{x}}_i\left(t\right)\right)\text{d}\left(x_i\left(t\right)-{\tilde{x}}_i\left(t\right)\right)}\\ +{\displaystyle \underset{i=1}{\overset{3}{\sum }}{\displaystyle \underset{j=1,j\ne i}{\overset{3}{\sum }}{g}_i{c}_{ij}\left(\left|x_j\left(t\right)-{\tilde{x}}_j\left(t\right)\right|-\left|x_j\left(t-{\tau }_{ij}\right)-{\tilde{x}}_j\left(t-{\tau }_{ij}\right)\right|\right)}}\\ \le {\displaystyle \underset{i=1}{\overset{3}{\sum }}\left[-g_{ij}{c}_{ii}\left|x_i\left(t\right)-{\tilde{x}}_i\left(t\right)\right|+{\displaystyle \underset{j=1,j\ne i}{\overset{3}{\sum }}{g}_{ij}{c}_{ij}\left|x_j\left(t\right)-{\tilde{x}}_j\left(t\right)\right|}\right]\text{d}t}.\end{array}$

According to Theorem 2.3 in [22] , we have

${\displaystyle \underset{i=1}{\overset{3}{\sum }}{\displaystyle \underset{j=1,j\ne i}{\overset{3}{\sum }}{g}_i{c}_{ij}\left|x_j\left(t\right)-{\tilde{x}}_j\left(t\right)\right|}}={\displaystyle \underset{i=1}{\overset{3}{\sum }}{\displaystyle \underset{j=1,j\ne i}{\overset{3}{\sum }}{g}_i{c}_{ij}\left|x_i\left(t\right)-{\tilde{x}}_i\left(t\right)\right|}}.$

As a consequence,

$E\left(V\left(t\right)\right)+{\displaystyle {\int }_0^t{\displaystyle \underset{i=1}{\overset{3}{\sum }}{g}_i[c_{ii}-{\displaystyle \underset{j=1,j\ne i}{\overset{3}{\sum }}\left|c_{ij}\right|}]}}E\left|x_i\left(s\right)-{\tilde{x}}_i\left(s\right)\right|\text{d}s\le V\left(0\right)<\infty .$

Therefore

$E\left|x_i\left(t\right)-{\tilde{x}}_i\left(t\right)\right|\in L^1\left[0,+\infty \right).$ (5.1)

According to (SM),

$\begin{array}{c}E\left(x_1\left(t\right)\right)=x_1\left(0\right)+{\displaystyle {\int }_0^t[E\left(x_1\left(s\right)\right)r_1-a_1{C}_0\left(t\right)E\left(x_1\left(s\right)\right)-c_{11}E{\left(x_1\left(s\right)\right)}^2}\\ -c_{12}E\left(x_1\left(s\right)x_2\left(s-{\tau }_{12}\right)\right)-c_{13}E\left(x_1\left(s\right)x_3\left(s-{\tau }_{13}\right)\right)]\text{d}s.\end{array}$

Clearly, $E\left(x_1\left(t\right)\right)$ is differentiable. By virtue of (3.1),

$\begin{array}{c}\frac{\text{d}E\left(x_1\left(t\right)\right)}{\text{d}t}=E\left(x_1\left(t\right)\right)r_1-a_1{C}_0\left(t\right)E\left(x_1\left(t\right)\right)-c{}_{11}E{\left(x_1\left(t\right)\right)}^2\\ -c_{12}E\left(x_1\left(t\right)x_2\left(t-{\tau }_{12}\right)\right)-c_{13}E\left(x_1\left(t\right)x_3\left(t-{\tau }_{13}\right)\right)\\ \le E\left(x_1\left(t\right)\right)r_1\le r_{1}K.\end{array}$

where $K>0$ is a constant. Consequently, $E\left(x_1\left(t\right)\right)$ is uniformly continuous. In the same way, $E\left(x_2\left(t\right)\right)$ and $E\left(x_3\left(t\right)\right)$ are also uniformly continuous functions. By virtue of (5.1) and Barbalat’s conclusion [23] , we derive the required assertion.

6. Numerical Simulations

Now we introduce some numerical figures to support Theorem 4.1 by using the Euler scheme [24] .

Considering the parameters as following:

$\{\begin{array}{l}\text{d}{x}_1\left(t\right)=x_1\left(t\right)\left[0.8-0.6x_1\left(t\right)-0.2x_2\left(t-1\right)-0.3x_3\left(t-2\right)\right]\text{d}t+0.3x_1\left(t\right)\text{d}{W}_1\left(t\right),\\ \text{d}{x}_2\left(t\right)=x_2\left(t\right)\left[0.8-0.3x_1\left(t-2\right)-0.6x_2\left(t\right)-0.2x_3\left(t-1\right)\right]\text{d}t+0.2x_1\left(t\right)\text{d}{W}_2\left(t\right),\\ \text{d}{x}_3\left(t\right)=x_3\left(t\right)\left[0.8-0.2x_1\left(t-3\right)-0.3x_2\left(t-2\right)-0.6x_3\left(t\right)\right]\text{d}t+0.1x_1\left(t\right)\text{d}{W}_3\left(t\right).\end{array}$ (6.1)

It is easy to get that:

$\begin{array}{l}C=0.1430,C_1=0.3340,C_2=0.0580,C_3=0.3340,\\ C_{23}=0.1400,C_{32}=0.0300,C_{31}=-0.1400,C_{33}=\mathrm{0.3000.}\end{array}$

Set ${\beta }_1^2/2=0.0450,{\beta }_2^2/2=0.0200$ and ${\beta }_3^2/2=0.0050$. Then by calculation, we have ${\tilde{C}}_1=0.0140>0$, ${\tilde{C}}_2=0.0074>0$, ${\tilde{C}}_3=0.0056>0$ and $\Gamma =4\times {10}^{-4}>0$. Furthermore, $C_3/{\tilde{C}}_3=59.6429>1$. Therefore all conditions in Theorem 4.1 have been checked. Then we have

$\begin{array}{l}\underset{t\to +\infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t{x}_1\left(s\right)\text{d}s}=\frac{C_1-{\tilde{C}}_1}{C}=2.3007,\\ \underset{t\to +\infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t{x}_2\left(s\right)\text{d}s}=\frac{C_2-{\tilde{C}}_2}{C}=0.3538,\\ \underset{t\to +\infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t{x}_3\left(s\right)\text{d}s}=\frac{C_3-{\tilde{C}}_3}{C}=\mathrm{2.2965.}\end{array}$

see Figure 1, which is obtained by applying the Milstein method [25] .

Figure 1 shows the simulations of the solutions of systems (SM), besides, ${\bar{x}}_1\left(t\right),{\bar{x}}_2\left(t\right)$ and ${\bar{x}}_3\left(t\right)$ are shown in Figure 1. We see that ${\bar{x}}_1\left(t\right),{\bar{x}}_2\left(t\right),{\bar{x}}_3\left(t\right)$ tend to constants, which is consistent with the results of Theorem 4.1. the solutions of systems (SM) fluctuate around a small zone. Thus, we think the system (SM) is permanent.

Since the parameters given above meet the hypothesis 2: $c_{11}>c_{12}+c_{13}$, $c_{22}>c_{21}+c_{23}$ and $c_{33}>c_{31}+c_{32}$. According to Theorem 5.1, we can get the system (SM) is global attractively, see Figure 2.

Figure 2 shows the simulations of the solutions of systems (SM), From Figure 2(a) and Figure 2(b), we can see that the solution of the system is globally attractively, whether with or without random perturbations.

7. Conclusions and Discussions

The dynamic relationship between predator and their preys has been and will continue to be one of the major themes in ecology due to its importance and

Figure 1. System (6.1) with ${\beta }_1^2/2=0.045,{\beta }_2^2/2=0.02$ and ${\beta }_3^2/2=0.005$. The red line and purplish red represent two prey population $x_1\left(t\right),x_2\left(t\right)$, respectively, the pink line represents the predator population $x_3\left(t\right)$. The yellow line represents ${\bar{x}}_1\left(t\right)$, the green line represent ${\bar{x}}_2\left(t\right)$ and the blue line represent ${\bar{x}}_3\left(t\right)$. Figure 1(a) and Figure 1(b) represent system with and without random perturbations, respectively. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

Figure 2. System (6.1) with $a_{11}=0.6>a_{12}+a_{13}=0.2+0.3=0.5$, $a_{22}=0.6>a_{21}+a_{23}=0.3+0.2=0.5$, $a_{33}=0.6>a_{31}+a_{32}=0.2+0.3=0.5$. Time series diagram of $x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)$. The green line represents one prey population $x_1(t),$ The purplish red line represents another prey population $x_2\left(t\right)$ and the blue line represents one predator population $x_3\left(t\right)$. Figure 2(a) and Figure 2(b) represent system with and without random perturbations, respectively.

universal existence. This paper is concerned with a delayed stochastic one-predator-two-prey population model in a polluted environment. Firstly, we show that there exists a unique positive solution in our system. Secondly, the permanence of this system is investigated. Conditions for the system to be permanent in time average are given. Our main result in part 4 reveals the impacts of stochastic perturbations on the persistence and extinction of every species. Finally, our results are confirmed by numerical simulation.

Some questions deserve further explorations. In the first place, it is significant to study the delay population system with other disturbance, such as Lévy jumps, or Markovian switching. Another problem is to consider population models with different functional responses, such as Holling II-IV type and Beddington-DeAngelis functional response. We leave these investigations for future work.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (11861027). This work was supported by the Natural Science Foundation of Guangxi (2016GXNSFAA380194).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Hallam, T.G., Clark, C.E. and Jordan, G.S. (1983) Effects of Toxicant on Population: A Qualitative Approach II. First Order Kinetics. Journal of Mathematical Biology, 18, 25-37. https://doi.org/10.1007/BF00275908
[2]	Hallam, T.G., Clark, C.E. and Lassiter, R.R. (1983) Effects of Toxicant on Population: A Qualitative Approach I. Equilibrium Environmental Exposure. Ecological Modelling, 18, 291-304. https://doi.org/10.1016/0304-3800(83)90019-4
[3]	Hallam, T.G. and Deluna, J.L. (1984) Effects of Toxicant on Populations: A Qualitative Approach III. Environmental and Food Chain Pathways. Journal of Theoretical Biology, 109, 411-429. https://doi.org/10.1016/S0022-5193(84)80090-9
[4]	Ji, C.Y., Jiang, D.Q. and Shi, N.Z. (2009) Analysis of a Predator-Prey Model with Modified Lesile-Gower and Holling-Type II. Schemes with Stochastic Perturbation. Journal of Mathematical Analysis and Applications, 359, 482-498. https://doi.org/10.1016/j.jmaa.2009.05.039
[5]	Ji, C.Y., Jiang, D.Q. and Li, X.Y. (2011) Qualitative Analysis of a Stochastic Ratio-Dependent Predator-Prey System. Journal of Computational and Applied Mathematics, 235, 1326-1341. https://doi.org/10.1016/j.cam.2010.08.021
[6]	Chi, M.N. and Zhao, W.C. (2018) Dynamical Analysis of Multi-Nutrient and Single Microorganism Chemostat Model in a Polluted Environment. Advances in Difference Equations, 2018, 120. https://doi.org/10.1186/s13662-018-1573-3
[7]	Li, Y.J. and Meng, X.Z. (2019) Dynamics of an Impulsive Stochastic Nonautonomous Chemostat Model with Two Different Growth Rates in a Polluted Environment. Discrete Dynamics in Nature and Society, 2019, Article ID: 5498569. https://doi.org/10.1155/2019/5498569
[8]	Liu, Y.L., Liu, Q. and Liu, Z.H. (2013) Dynamical Behaviors of a Stochastic Delay Logistic System with Impulsive Toxicant Input in a Polluted Environment. Journal of Theoretical Biology, 329, 1-5. https://doi.org/10.1016/j.jtbi.2013.03.005
[9]	Song, Y.L., Peng, Y.H. and Wei, J.J. (2008) Bifurcations for a Predator-Prey System with Two Delays. Journal of Mathematical Analysis and Applications, 337, 466-479. https://doi.org/10.1016/j.jmaa.2007.04.001
[10]	Liu, M. and Bai, C. (2016) Dynamics of a Stochastic One-Prey Two-Predator Model with Lévy Jumps. Applied Mathematics and Computation, 284, 308-321. https://doi.org/10.1016/j.amc.2016.02.033
[11]	Liu, M. and Fan, M. (2017) Stability in Distribution of a Three-Species Stochastic Cascade Predator-Prey System with Time Delays. IMA Journal of Applied Mathematics, 82, 396-423.
[12]	Liu, M., Bai, C. and Jin, Y. (2017) Population Dynamical Behavior of a Two-Predator One-Prey Stochastic Model with Time Delay. Discrete & Continuous Dynamical Systems, 37, 2513-2538. https://doi.org/10.3934/dcds.2017108
[13]	Kuang, Y. (1993) Delay Differential Equations with Applications in Population Dynamics. Academic Press, Boston.
[14]	Han, Q.X., Jiang, D.Q. and Ji, C.Y. (2014) Analysis of a Delayed Stochastic Predator-Prey Model in a Polluted Environment. Applied Mathematical Modelling, 33, 3067-3080. https://doi.org/10.1016/j.apm.2013.11.014
[15]	Wang, K. (2010) Stochastic Biomathematical Model. Science Press, Beijing, 174-175.
[16]	Ikeda, N. and Watanabe, S. (1977) A Comparison Theorem for Solutions of Stochastic Differential Equations and Its Applications. Osaka Journal of Mathematics, 14, 619-633.
[17]	Jiang, D. and Shi, N. (2005) A Note on Nonautonomous Logistic Equation with Random Perturbation. Journal of Mathematical Analysis and Applications, 303, 164-172. https://doi.org/10.1016/j.jmaa.2004.08.027
[18]	Liu, M., Wang, K. and Wu, Q. (2011) Survival Analysis of Stochastic Competitive Models in a Polluted Environment and Stochastic Competitive Exclusion Principle. Bulletin of Mathematical Biology, 73, 1969-2012. https://doi.org/10.1007/s11538-010-9569-5
[19]	Jiang, D., Ji, C., Li, X. and O’Regand, D. (2012) Analysis of Autonomous Lotka-Volterra Competition Systems with Random Perturbation. Journal of Mathematical Analysis and Applications, 390, 582-595. https://doi.org/10.1016/j.jmaa.2011.12.049
[20]	Bao, J., Mao, X., Yin, G. and Yuan, C. (2011) Competitive Lotka-Volterra Population Dynamics with Jumps. Nonlinear Analysis, 74, 6601-6616. https://doi.org/10.1016/j.na.2011.06.043
[21]	Moon, J. (1970) Counting Labelled Tress. Canadian Mathematical Congress, Montreal.
[22]	Li, M.Y. and Shuai, Z. (2010) Global-Stability Problem for Coupled Systems of Differential Equations on Networks. Journal of Differential Equations, 248, 1-20. https://doi.org/10.1016/j.jde.2009.09.003
[23]	Barbalat, I. (1959) Systèms d’équations différentielles d’oscillations non linéaires. Revue Roumaine de Mathématique Pures et Appliquées, 4, 267-270.
[24]	Protter, P. and Talay, D. (1997) The Euler Scheme for Lévy Driven Stochastic Differential Equations. Annals of Probability, 25, 393-423. https://doi.org/10.1214/aop/1024404293
[25]	Higham, D.J. (2001) An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations. SIAM Review, 43, 525-546. https://doi.org/10.1137/S0036144500378302