Dynamical analysis of a modified Leslie–Gower Holling-type II predator-prey stochastic model in polluted environments with interspecific competition and impulsive toxicant input

Abstract

In this paper, we use a mean-reverting Ornstein-Uhlenbeck process to simulate the stochastic perturbations in the environment, and then a modified Leslie–Gower Holling-type II predator-prey stochastic model in a polluted environment with interspecific competition and pulse toxicant input is proposed. Through constructing V-function and applying $\mathrm{It}{\hat{\mathrm{o}}}^{\prime }\mathrm{s}$ formula, the sharp sufficient conditions including strongly persistent in the mean, persistent in the mean and extinction are established. In addition, the theoretical results are verified by numerical simulation.

Keywords:

• Stochastic models
• strongly persistent in the mean
• persistence in the mean
• extinction
• intensity of white noise
• toxicant effect

MATHEMATICS Subject Classification (2000):

• 60H10
• 60H30
• 92D25

1. Introduction

In recent years, with the continuous development of industry, agriculture and fisheries as well as the discharge of waste water, waste gas and wastes, the ecological environment has been seriously damaged. Today's society emphasizes the importance of protecting the environment because these toxins in the polluted environment pose a serious threat to the survival of species such as plankton, animal populations and even human beings themselves. All this has led more scholars to investigate the effects of toxins on populations in polluted environments. For example, Nelson [34] studied the early oil pollution of European seas, Jensen and Marshall [16] described laboratory populations of Daphnia pulex exposed to chronic radiation stress, and Shukla et al. [37] proposed a dynamical model to assess the effects of industrialization on forest biomass degradation in a particular area.

At the same time, many researchers have studied the effects of toxic substances and pollutants on ecological populations from a realistic perspective. In particular, Hallam and his colleagues [9–11] have opened the door to the study of environmental toxins by publishing three papers in a row that explored the effects of toxins on deterministic models of ecosystems. From then on, many deterministic single or multiple population models with toxic effects have been proposed and analysed (see e.g. [5,7,12,28,32,36]), however the underlying assumption is that exogenous inputs of toxins are continuous. Particularly, in real life, toxins are often released into the environment in pulses. A typical example is the regular spraying of pesticides in agriculture and forestry [21], and another example is the massive discharge of heavy metals in industry [20]. Therefore, on the basis of a deterministic population model, some authors have investigated the effects of toxins with pulsed inputs on the population in polluted environments(see e.g. [22,26,31,42]).

Those important and useful studies of deterministic models provide a great way to understand the effects of pollution, but in the real world, populations do not exist in isolation. Therefore, the study of the persistence and extinction of two or more interacting populations under the action of toxic substances has more important biological significance. In recent years, several scholars have proposed more realistic predation models which should take the functional response into account [4,38]. Holling type II functional response is used to measure the average feeding rate of a predator when it spends time looking for and handing the prey, which has served as the basis of many predator-prey theories. As a result, Aziz-Alaoui and Okiye [1] proposed a famous predator-prey model incorporating a modified version of the Leslie–Gower functional response and a Holling-type II one, the model can be shown as follows: $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x(t)}{\mathrm{d}t}}=x(t)(r_1-ax(t)-{\displaystyle \frac{cy(t)}{h+x(t)}}),\\ {\displaystyle \frac{\mathrm{d}y(t)}{\mathrm{d}t}}=y(t)(r_2-{\displaystyle \frac{fy(t)}{h+x(t)}}),\end{array}$where $x(t)$ and $y(t)$ stand for the size of the prey population and the predator population respectively; $r_1$ and $r_2$ denote the intrinsic growth rate of the prey and the predator respectively. a is the internal competition coefficient of the prey; c represents the per capita reduction rate of the prey; h indicates the degree of environmental protection to predator $y(t)$, that is to say, it represents the capacity of food provision of the environment except prey $x(t)$; f is the per capita reduction rate of the predator. A growing number of scholars based on the above model have studied the possibility of a reciprocal relationship between the decline of predator populations and the per capita availability of prey(see e.g. [17,18,29,35,39,41]).

Note that the above model is based on a deterministic model independent of random environment fluctuations. In fact, the growth of populations is inevitably affected by environmental noise. For example, May [33] pointed out that the growth rate in population system should be stochastic and the solution of system fluctuated randomly around some mean value. On this basis, many studies focussed on random population models in polluted environment and some good results were obtained(see e.g. [13,24,25,30,40,43,44]). In particular, Liu et al. [27] proposed the following stochastic predator-prey model in a polluted environment with pulse toxicant input (1) $\{\begin{array}{l}{\displaystyle \mathrm{d}x(t)=x(t)[r_{10}-r_{11}{C}_{10}(t)-ax(t)-{\displaystyle \frac{cy(t)}{h+x(t)}}]\mathrm{d}t+{\sigma }_{1}x(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}y(t)=y(t)[r_{20}-r_{21}{C}_{20}(t)-{\displaystyle \frac{fy(t)}{h+x(t)}}]\mathrm{d}t+{\sigma }_{2}y(t)\mathrm{d}{B}_2(t),}\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+}.}\end{array}$(1) where $C_{10}(t)$ and $C_{20}(t)$ denote the concentration of the toxicant in the organism of the prey species, the predator species at time t, respectively; $\mathrm{\Delta }\zeta (t)=\zeta (t^{+})-\zeta (t)$; $k_i$ stands for the organisms net uptake rate of toxicant from the environment; $l_i$ denotes the loss rate of the toxicant in the organism due to the processes such as egestion; $C_e(t)$ denotes the concentrations of the toxicant in the environment at time t; g is the loss rate of toxicant from environment; q is the toxicant input amount at every time; $\gamma (>0)$ stands for the period of the impulsive input of toxicant. ${\sigma }_i^2$, i = 1, 2 is the intensity of white noise; $B_1(t)$ and $B_2(t)$ are mutually independent standard Brownian motions defined on a complete probability space $(\mathrm{\Omega },\mathcal{F},\mathcal{P})$ with a filtration $\{{\mathcal{F}}_t{\}}_{t\in R_{+}}$.

Obviously, there is a limitation in model (1(1) $\{\begin{array}{l}{\displaystyle \mathrm{d}x(t)=x(t)[r_{10}-r_{11}{C}_{10}(t)-ax(t)-{\displaystyle \frac{cy(t)}{h+x(t)}}]\mathrm{d}t+{\sigma }_{1}x(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}y(t)=y(t)[r_{20}-r_{21}{C}_{20}(t)-{\displaystyle \frac{fy(t)}{h+x(t)}}]\mathrm{d}t+{\sigma }_{2}y(t)\mathrm{d}{B}_2(t),}\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+}.}\end{array}$(1) ), that is, in model (1(1) $\{\begin{array}{l}{\displaystyle \mathrm{d}x(t)=x(t)[r_{10}-r_{11}{C}_{10}(t)-ax(t)-{\displaystyle \frac{cy(t)}{h+x(t)}}]\mathrm{d}t+{\sigma }_{1}x(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}y(t)=y(t)[r_{20}-r_{21}{C}_{20}(t)-{\displaystyle \frac{fy(t)}{h+x(t)}}]\mathrm{d}t+{\sigma }_{2}y(t)\mathrm{d}{B}_2(t),}\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+}.}\end{array}$(1) ) there exists an assumption that the growth rate $r_{i0}(i=1,2)$ is linear with respect to the Gaussian white noise ${\hat{r}}_{i0}(t):=r_{i0}+{\sigma }_i\frac{\mathrm{d}{B}_i(t)}{\mathrm{d}t},i=1,2.$Integrating on the interval $[0,T]$ results in ${\overline{r}}_{i0}:=\frac{1}{T}{\int }_0^T{\hat{r}}_{i0}(t)\mathrm{d}t=r_{i0}+{\sigma }_i\frac{B_i(T)}{T}\sim N(r_{i0},{\sigma }_i^2/T).$Hence, the variance of the average per capita growth rate ${\overline{r}}_{i0}$ over an interval of length T tends to ∞ as $T\to 0$. It may not be enough to describe the real situation. In recent years, some authors [6,45] have claimed that using the mean-reverting Ornstein-Uhlenbeck process is a more appropriate approach to incorporate the environment perturbations, and due to this approach, one has $\mathrm{d}{\hat{r}}_{i0}(t)={\alpha }_i(r_{i0}-{\hat{r}}_{i0}(t))\mathrm{d}t+{\xi }_i\mathrm{d}{B}_i(t),i=1,2,$i.e. $\begin{array}{rl}{\hat{r}}_{i0}(t)& =r_{i0}+({\tilde{r}}_{i0}-r_{i0})e^{-{\alpha }_{i}t}+{\xi }_i{\int }_0^{te-{\alpha }_i(t-s)}\mathrm{d}{B}_i(s)\\ & =r_{i0}+({\tilde{r}}_{i0}-r_{i0})e^{-{\alpha }_{i}t}+{\sigma }_i(t)\frac{\mathrm{d}{B}_i(t)}{\mathrm{d}t},i=1,2,\end{array}$where ${\tilde{r}}_{i0}={\hat{r}}_{i0}(0)$, ${\sigma }_i(t)=\frac{{\xi }_i}{\sqrt{2{\alpha }_i}}\sqrt{1-e^{-2{\alpha }_{i}t}}$, ${\alpha }_i>0$ represents the speed of reversion, ${\xi }_i^2$ is the intensity of stochastic perturbations. Based on the ideas above, we can derive a new version as follows: (2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \mathrm{d}x(t)=x(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-ax(t)-\frac{cy(t)}{h+x(t)}]\mathrm{d}t}\\ +{\sigma }_1(t)x(t)\mathrm{d}{B}_1(t),\\ {\displaystyle \mathrm{d}y(t)=y(t)[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-\frac{fy(t)}{h+x(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)x(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+}.}\end{array}\end{array}$(2) In the real ecological environment, the predator-prey relationship is a complex biological relationship. Therefore, Gao et al. [8] studied the two-predator one-prey model with pulsed toxin input based on model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \mathrm{d}x(t)=x(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-ax(t)-\frac{cy(t)}{h+x(t)}]\mathrm{d}t}\\ +{\sigma }_1(t)x(t)\mathrm{d}{B}_1(t),\\ {\displaystyle \mathrm{d}y(t)=y(t)[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-\frac{fy(t)}{h+x(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)x(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+}.}\end{array}\end{array}$(2) ). In the model in paper [8], there is no competition between the two predators. In this paper, we study an impulsive stochastic model with one prey and two predators where the two predators have a competitive relationship. The model is as follows: (3) $\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,3,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\end{array}$(3) where $y_1(t)$ denotes the size of the prey population; $y_2(t)$ and $y_3(t)$ denote the size of the predators population; $r_{i0},i=1,2,3$ denotes the growth rate of population $y_i$; $c_2$ and $c_3$ mean the per capita reduction rate of the prey due to the capture by $y_2$ and $y_3$ respectively; $f_2$ and $f_3$ have similar meanings to $c_2$ and $c_3$ respectively; $h_2$ and $h_3$ represent the degree of environmental protection to predators $y_2(t)$ and $y_3(t)$, respectively.

Remark 1.1

In model (3(3) $\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,3,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\end{array}$(3) ), the parameters $r_{i0},{\tilde{r}}_{i0},k_i,l_i(i=1,2,3)$, $f_i,h_i,c_i(i=2,3)$ and a satisfy $0<r_{i0},{\tilde{r}}_{i0},k_i,l_i,f_i,h_i,c_i,a\le 1$. The parameter ${\sigma }_i^2$ is the intensity of the white noise on the growth rate of species i, thus ${\sigma }_i^2>0,i=1,2,3$. The parameter $\gamma ,q>0$. $r_{i1}$ represents the decreasing rate of the intrinsic growth rate associated with the uptake of the toxicant. Thus $r_{i1}>0$. $a_{23},a_{32}$ denote the inter-specific competition coefficients of two predators. Thus $a_{23}>0,a_{32}>0$. In addition, $C_{i0}$ and $C_e$ stand for the concentrations of the toxicant, therefore $0\le C_{i0}(t)\le 1$ and $0\le C_e(t)\le 1$ for $t\ge 0$. As a result, the following conditions need to be met: $k_i\le l_i$, $q\le 1-e^{-g\gamma }$, i = 1, 2, 3.

As far as we know, a very little amount of work has been done about the stochastic three species predation model with toxicants effect, and little is known of the survival impact of competing predators in a polluted environment. In essence, the proof methods and techniques for having two predator species are similar to those for having more predator species. Motivated by these, we deeply analyse the properties of model (3(3) $\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,3,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\end{array}$(3) ).

The arrangement of this paper is as follows. In Section 2, the threshold of persistence and extinction for each species is obtained. In Section 3, we carry out some numerical simulations to verify the theoretical results. Finally, we give some conclusions in Section 4.

2. Main results

For the sake of convenience and simplicity, we define the following notations: $\begin{array}{rl}& R_{+}^3=\{z\in R^3|z_i>0,i=1,2,3\},b_i(t)=r_{i0}-\frac{{\xi }_i^2}{4{\alpha }_i}+\frac{{\xi }_i^2}{4{\alpha }_i}{e}^{-2{\alpha }_{i}t}-r_{i1}{C}_{i0}(t),\\ & K_i=\frac{q{k}_i}{g{l}_i},{\overline{b}}_i=\underset{t\to +\mathrm{\infty }}{lim}{t}^{-1}{\int }_0^t{b}_i(s)\mathrm{d}s=r_{i0}-\frac{{\xi }_i^2}{4{\alpha }_i}-\frac{r_{i1}{K}_i}{\gamma },i=1,2,3.\\ & ⟨y(t)⟩=t^{-1}{\int }_0^{t}y(s)\mathrm{d}s,y^{\ast }=\underset{t\to +\mathrm{\infty }}{lim sup}y(t),y_{\ast }=\underset{t\to +\mathrm{\infty }}{lim inf}y(t)\end{array}$

Lemma 2.1

[26,42]

Consider the following subsystem of (3(3) $\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,3,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\end{array}$(3) ): (4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{C}_{10}(t)}{\mathrm{d}t}=k_1{C}_e(t)-l_1{C}_{10}(t),}\\ {\displaystyle \frac{\mathrm{d}{C}_{20}(t)}{\mathrm{d}t}=k_2{C}_e(t)-l_2{C}_{20}(t),}\\ {\displaystyle \frac{\mathrm{d}{C}_{30}(t)}{\mathrm{d}t}=k_3{C}_e(t)-l_3{C}_{30}(t),}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_{i0}(t)=0,\mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\\ {\displaystyle 0\le C_{i0}(0)\le 1,0\le C_e(0)\le 1,i=1,2,3.}\end{array}\end{array}$(4)

The model has a unique positive γ-periodic solution$({\tilde{C}}_{10}(t),{\tilde{C}}_{20}(t),{\tilde{C}}_{30}(t),{\tilde{C}}_e(t){)}^T$ which satisfies $\underset{t\to +\mathrm{\infty }}{lim}{t}^{-1}{\int }_0^t{\tilde{C}}_{i0}(s)\mathrm{d}s=\frac{k_{i}q}{g{l}_i\gamma }=\frac{K_i}{\gamma },i=1,2,3.$Notice that in model (4(4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{C}_{10}(t)}{\mathrm{d}t}=k_1{C}_e(t)-l_1{C}_{10}(t),}\\ {\displaystyle \frac{\mathrm{d}{C}_{20}(t)}{\mathrm{d}t}=k_2{C}_e(t)-l_2{C}_{20}(t),}\\ {\displaystyle \frac{\mathrm{d}{C}_{30}(t)}{\mathrm{d}t}=k_3{C}_e(t)-l_3{C}_{30}(t),}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_{i0}(t)=0,\mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\\ {\displaystyle 0\le C_{i0}(0)\le 1,0\le C_e(0)\le 1,i=1,2,3.}\end{array}\end{array}$(4) ), $C_{i0}(t)$ and $C_e(t)$ can be explicitly solved. Then we only need to consider the first three equations in model (3(3) $\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,3,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\end{array}$(3) ), which is expressed as follows (5) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t).\end{array}\end{array}$(5)

Lemma 2.2

For arbitrary$(y_1(0),y_2(0),y_3(0))\in R_{+}^3$, model (5(5) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t).\end{array}\end{array}$(5) ) possesses a unique positive solution $(y_1(t),y_2(t),y_3(t))\in R_{+}^3$ for all $t\ge 0$ a.s.(almost surely).

Proof.

Pay attention to the following system: (6) $\begin{array}{r}\{\begin{array}{l}\mathrm{d}u(t)=[b_1(t)+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-a{e}^{u(t)}-\frac{c_2{e}^{v(t)}}{h_2+e^{u(t)}}-\frac{c_3{e}^{w(t)}}{h_3+e^{u(t)}}]\mathrm{d}t\\ +{\sigma }_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}v(t)=[b_2(t)+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-a_{23}{e}^{w(t)}-\frac{f_2{e}^{v(t)}}{h_2+e^{u(t)}}]\mathrm{d}t\\ +{\sigma }_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}w(t)=[b_3(t)+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-a_{32}{e}^{v(t)}-\frac{f_3{e}^{w(t)}}{h_3+e^{u(t)}}]\mathrm{d}t\\ +{\sigma }_3(t)\mathrm{d}{B}_3(t),\end{array}\end{array}$(6) and $u(0)=ln{y}_1(0)$, $v(0)=ln{y}_2(0)$, $w(0)=ln{y}_3(0)$. Due to the fact that the coefficients of model (6(6) $\begin{array}{r}\{\begin{array}{l}\mathrm{d}u(t)=[b_1(t)+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-a{e}^{u(t)}-\frac{c_2{e}^{v(t)}}{h_2+e^{u(t)}}-\frac{c_3{e}^{w(t)}}{h_3+e^{u(t)}}]\mathrm{d}t\\ +{\sigma }_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}v(t)=[b_2(t)+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-a_{23}{e}^{w(t)}-\frac{f_2{e}^{v(t)}}{h_2+e^{u(t)}}]\mathrm{d}t\\ +{\sigma }_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}w(t)=[b_3(t)+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-a_{32}{e}^{v(t)}-\frac{f_3{e}^{w(t)}}{h_3+e^{u(t)}}]\mathrm{d}t\\ +{\sigma }_3(t)\mathrm{d}{B}_3(t),\end{array}\end{array}$(6) ) satisfy the local Lipschitz condition, model (6(6) $\begin{array}{r}\{\begin{array}{l}\mathrm{d}u(t)=[b_1(t)+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-a{e}^{u(t)}-\frac{c_2{e}^{v(t)}}{h_2+e^{u(t)}}-\frac{c_3{e}^{w(t)}}{h_3+e^{u(t)}}]\mathrm{d}t\\ +{\sigma }_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}v(t)=[b_2(t)+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-a_{23}{e}^{w(t)}-\frac{f_2{e}^{v(t)}}{h_2+e^{u(t)}}]\mathrm{d}t\\ +{\sigma }_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}w(t)=[b_3(t)+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-a_{32}{e}^{v(t)}-\frac{f_3{e}^{w(t)}}{h_3+e^{u(t)}}]\mathrm{d}t\\ +{\sigma }_3(t)\mathrm{d}{B}_3(t),\end{array}\end{array}$(6) ) possesses a unique local solution $(u(t),v(t),w(t){)}^T$ on $[0,{\tau }_{\ast })$, where ${\tau }_{\ast }$ means the explosion time. From the $\mathrm{It}\hat{o}\mathrm{\prime }s$ formula, we derive that $(y_1(t),y_2(t),y_3(t))=(e^{u(t)},e^{v(t)},e^{w(t)})$ is the unique local positive solution of model (5(5) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t).\end{array}\end{array}$(5) ).

Now let us verify that ${\tau }_{\ast }=+\mathrm{\infty }$. Consider the following systems: (7) $\begin{array}{rl}\mathrm{d}\varphi (t)& =\varphi (t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a\varphi (t)]\mathrm{d}t\\ & +{\sigma }_1\varphi (t)\mathrm{d}{B}_1(t),\varphi (0)=y_1(0);\end{array}$(7) (8) $\begin{array}{rl}\mathrm{d}M(t)& =M(t)[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-\frac{f_3}{h_3+\varphi (t)}M(t)]\mathrm{d}t\\ & +{\sigma }_{3}M(t)\mathrm{d}{B}_3(t),M(0)=y_3(0).\end{array}$(8) (9) $\begin{array}{rl}\mathrm{d}N(t)& =N(t)[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-\frac{f_2}{h_2+\varphi (t)}N(t)]\mathrm{d}t\\ & +{\sigma }_{2}N(t)\mathrm{d}{B}_2(t),N(0)=y_2(0);\end{array}$(9) (10) $\begin{array}{rl}\mathrm{d}n(t)& =n(t)[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-\frac{f_2}{h_2}n(t)-a_{23}M(t)]\mathrm{d}t\\ & +{\sigma }_{2}n(t)\mathrm{d}{B}_2(t),n(0)=y_2(0);\end{array}$(10) (11) $\begin{array}{rl}\mathrm{d}m(t)& =m(t)[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-\frac{f_3}{h_3}m(t)-a_{32}N(t)]\mathrm{d}t\\ & +{\sigma }_{3}m(t)\mathrm{d}{B}_3(t),m(0)=y_3(0);\end{array}$(11) On the basis of the comparison theorem for stochastic differential equations [15],we get for $t\in [0,{\tau }_{\ast })$, (12) $y_1(t)\le \varphi (t),n(t)\le y_2(t)\le N(t),m(t)\le y_3(t)\le M(t),a.s.$(12) According to Theorem 2.2 in Jiang and Shi [19], we get (13) $\begin{array}{rl}\varphi (t)& =\frac{e^{{\int }_0^t{b}_1(s)\mathrm{d}s-\frac{{\tilde{r}}_{10}-r_{10}}{{\alpha }_1}(e^{-{\alpha }_{1}t}-1)+{\int }_0^t{\sigma }_1(s)\mathrm{d}{B}_1(s)}}{y_1^{-1}(0)+a{\int }_0^t{e}^{{\int }_0^s{b}_1(\tau )d\tau -\frac{{\tilde{r}}_{10}-r_{10}}{{\alpha }_1}(e^{-{\alpha }_{1}s}-1)+{\int }_0^s{\sigma }_1(\tau )\mathrm{d}{B}_1(\tau )}\mathrm{d}s},\end{array}$(13) (14) $\begin{array}{rl}M(t)& =\frac{e^{{\int }_0^t{b}_3(s)\mathrm{d}s-\frac{{\tilde{r}}_{30}-r_{30}}{{\alpha }_3}(e^{-{\alpha }_{3}t}-1)+{\int }_0^t{\sigma }_3(s)\mathrm{d}{B}_3(s)}}{y_3^{-1}(0)+{\int }_0^t\frac{f_3}{h_3+\varphi (s)}{e}^{{\int }_0^s{b}_3(\tau )d\tau -\frac{{\tilde{r}}_{30}-r_{30}}{{\alpha }_3}(e^{-{\alpha }_{3}s}-1)+{\int }_0^s{\sigma }_3(\tau )\mathrm{d}{B}_3(\tau )}\mathrm{d}s}.\end{array}$(14) (15) $\begin{array}{rl}N(t)& =\frac{e^{{\int }_0^t{b}_2(s)\mathrm{d}s-\frac{{\tilde{r}}_{20}-r_{20}}{{\alpha }_2}(e^{-{\alpha }_{2}t}-1)+{\int }_0^t{\sigma }_2(s)\mathrm{d}{B}_2(s)}}{y_2^{-1}(0)+{\int }_0^t\frac{f_2}{h_2+\varphi (s)}{e}^{{\int }_0^s{b}_2(\tau )d\tau -\frac{{\tilde{r}}_{20}-r_{20}}{{\alpha }_2}(e^{-{\alpha }_{2}s}-1)+{\int }_0^s{\sigma }_2(\tau )\mathrm{d}{B}_2(\tau )}\mathrm{d}s},\end{array}$(15) (16) $\begin{array}{rl}n(t)& =\frac{e^{{\int }_0^t[b_2(s)-a_{23}M(s)]\mathrm{d}s-\frac{{\tilde{r}}_{20}-r_{20}}{{\alpha }_2}(e^{-{\alpha }_{2}t}-1)+{\int }_0^t{\sigma }_2(s)\mathrm{d}{B}_2(s)}}{y_2^{-1}(0)+\frac{f_2}{h_2}{\int }_0^t{e}^{{\int }_0^s[b_2(\tau )-a_{23}M(\tau )]d\tau -\frac{{\tilde{r}}_{20}-r_{20}}{{\alpha }_2}(e^{-{\alpha }_{2}s}-1)+{\int }_0^s{\sigma }_2(\tau )\mathrm{d}{B}_2(\tau )}\mathrm{d}s},\end{array}$(16) (17) $\begin{array}{rl}m(t)& =\frac{e^{{\int }_0^t[b_3(s)-a_{32}N(s)]\mathrm{d}s-\frac{{\tilde{r}}_{30}-r_{30}}{{\alpha }_3}(e^{-{\alpha }_{3}t}-1)+{\int }_0^t{\sigma }_3(s)\mathrm{d}{B}_3(s)}}{y_3^{-1}(0)+\frac{f_3}{h_3}{\int }_0^t{e}^{{\int }_0^s[b_3(\tau )-a_{32}N(\tau )]d\tau -\frac{{\tilde{r}}_{30}-r_{30}}{{\alpha }_3}(e^{-{\alpha }_{3}s}-1)+{\int }_0^s{\sigma }_3(\tau )\mathrm{d}{B}_3(\tau )}\mathrm{d}s},\end{array}$(17) Due to the fact that $\varphi (t),n(t),N(t),m(t)$ and $M(t)$ are global, we can know that ${\tau }_{\ast }=+\mathrm{\infty }$.

Lemma 2.3

[30]

Let $X(t)\in C(\mathrm{\Omega }\times [0,+\mathrm{\infty }),R_{+})$, where $C(\mathrm{\Omega }\times [0,+\mathrm{\infty }),R_{+})$ denotes the family of all continuous positive-valued functions defined on $\mathrm{\Omega }\times [0,+\mathrm{\infty })$.

(i)	If there exist three positive constants $t_0$,β and ${\beta }_0$ such that for all $t\ge t_0$, $\ln X(t)\le \beta t-{\beta }_0{\int }_0^{t}X(s)\mathrm{d}s+F(t),$where $F(t)/t\to 0$ as $t\to +\mathrm{\infty }$, then $⟨X(t){⟩}^{\ast }\le \beta /{\beta }_0,a.s.$
(ii)	If there exist three positive constants $t_0$,β and ${\beta }_0$ such that for all $t\ge t_0$, $\ln X(t)\ge \beta t-{\beta }_0{\int }_0^{t}X(S)\mathrm{d}s+F(t),$where $F(t)/t\to 0$ as $t\to +\mathrm{\infty }$, then $⟨X(t){⟩}_{\ast }\ge \beta /{\beta }_0,a.s.$

Definition 2.1

(1) $y(t)$ is said to be extinct if $\underset{t\to +\mathrm{\infty }}{lim}y(t)=0$, a.s. (2) $y(t)$ is said to be strongly persistent in the mean if $⟨y(t){⟩}_{\ast }>0$, a.s. (3) $y(t)$ is said to be persistent in the mean if $\underset{t\to +\mathrm{\infty }}{lim}⟨y(t)⟩=C>0$, a.s.

Theorem 2.1

For model (5(5) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t).\end{array}\end{array}$(5) ), if ${\overline{b}}_i<0,$ then $y_i$ become extinct, $i.e.$ $\underset{t\to +\mathrm{\infty }}{lim}{y}_i(t)=0,$ a.s., $i=1,2,3.$

Proof.

Applying $\mathrm{It}{\hat{\mathrm{o}}}^{\prime }\mathrm{s}$ formula to model (5(5) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t).\end{array}\end{array}$(5) )gives $\begin{array}{rl}\mathrm{d}\ln y_1(t)& =(b_1(t)+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-a{y}_1(t)-\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)})\mathrm{d}t\\ & +{\sigma }_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}\ln y_2(t)& =(b_2(t)+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)})\mathrm{d}t+{\sigma }_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}\ln y_3(t)& =(b_3(t)+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)})\mathrm{d}t+{\sigma }_3(t)\mathrm{d}{B}_3(t).\end{array}$As a consequence, (18) $\begin{array}{rl}\ln y_1(t)-\ln y_1(0)& ={\int }_0^t{b}_1(s)\mathrm{d}s-\frac{{\tilde{r}}_{10}-r_{10}}{{\alpha }_1}(e^{-{\alpha }_{1}t}-1)-a{\int }_0^t{y}_1(s)\mathrm{d}s\\ & -c_2{\int }_0^t\frac{y_2(s)}{h_2+y_1(s)}\mathrm{d}s-c_3{\int }_0^t\frac{y_3(s)}{h_3+y_1(s)}\mathrm{d}s\\ & +{\int }_0^t{\sigma }_1(s)\mathrm{d}{B}_1(s),\end{array}$(18) (19) $\begin{array}{rl}\ln y_2(t)-\ln y_2(0)& ={\int }_0^t{b}_2(s)\mathrm{d}s-\frac{{\tilde{r}}_{20}-r_{20}}{{\alpha }_2}(e^{-{\alpha }_{2}t}-1)-a_{23}{\int }_0^t{y}_3(s)\mathrm{d}s\\ & -f_2{\int }_0^t\frac{y_2(s)}{h_2+y_1(s)}\mathrm{d}s+{\int }_0^t{\sigma }_2(s)\mathrm{d}{B}_2(s),\end{array}$(19) (20) $\begin{array}{rl}\ln y_3(t)-\ln y_3(0)& ={\int }_0^t{b}_3(s)\mathrm{d}s-\frac{{\tilde{r}}_{30}-r_{30}}{{\alpha }_3}(e^{-{\alpha }_{3}t}-1)-a_{32}{\int }_0^t{y}_2(s)\mathrm{d}s\\ & -f_3{\int }_0^t\frac{y_3(s)}{h_3+y_1(s)}\mathrm{d}s+{\int }_0^t{\sigma }_3(s)\mathrm{d}{B}_3(s).\end{array}$(20) We derive from (18(18) $\begin{array}{rl}\ln y_1(t)-\ln y_1(0)& ={\int }_0^t{b}_1(s)\mathrm{d}s-\frac{{\tilde{r}}_{10}-r_{10}}{{\alpha }_1}(e^{-{\alpha }_{1}t}-1)-a{\int }_0^t{y}_1(s)\mathrm{d}s\\ & -c_2{\int }_0^t\frac{y_2(s)}{h_2+y_1(s)}\mathrm{d}s-c_3{\int }_0^t\frac{y_3(s)}{h_3+y_1(s)}\mathrm{d}s\\ & +{\int }_0^t{\sigma }_1(s)\mathrm{d}{B}_1(s),\end{array}$(18) ) that, for sufficiently large t, (21) $t^{-1}\ln \frac{y_1(t)}{y_1(0)}\le {\overline{b}}_1+\epsilon +t^{-1}{\int }_0^t{\sigma }_1(s)\mathrm{d}{B}_1(s).$(21) Note that $\underset{t\to +\mathrm{\infty }}{lim}{t}^{-1}{\int }_0^t{\sigma }_1(s)\mathrm{d}{B}_1(s)=0$ and ${\overline{b}}_1+\epsilon <0$, thus $\underset{t\to +\mathrm{\infty }}{lim}{y}_1(t)=0,$ a.s. In the same way, ${\overline{b}}_2<0$ means that $\underset{t\to +\mathrm{\infty }}{lim}{y}_2(t)=0,$ a.s., ${\overline{b}}_3<0$ means that $\underset{t\to +\mathrm{\infty }}{lim}{y}_3(t)=0,$ a.s.

Theorem 2.2

For model (5(5) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t).\end{array}\end{array}$(5) ), if ${\overline{b}}_1<0,{\overline{b}}_i>0$ and ${\overline{b}}_j<0$, then $y_1,y_j$ become extinct and $y_i$ is persistent in the mean, $i.e.$ $\underset{t\to +\mathrm{\infty }}{lim}⟨y_i(t)⟩=\frac{h_i{\overline{b}}_i}{f_i},$ a.s. $i,j(i\ne j)=2,3$

Here we only prove the case of ${\overline{b}}_2>0$, and the proof of ${\overline{b}}_3>0$ can be obtained similarly.

Proof.

Since ${\overline{b}}_1<0,$ ${\overline{b}}_3<0$, we then deduce from Theorem 2.1 that $\underset{t\to +\mathrm{\infty }}{lim}{y}_1(t)=0,\underset{t\to +\mathrm{\infty }}{lim}{y}_3(t)=0,$ a.s. Then for sufficiently large t, (22) $\begin{array}{rl}& \ln y_2(t)\le ({\overline{b}}_2+\epsilon )t-\frac{f_2}{h_2+\epsilon }{\int }_0^t{y}_2(s)\mathrm{d}s+{\int }_0^t{\sigma }_2(s)\mathrm{d}{B}_2(s),\end{array}$(22) (23) $\begin{array}{rl}& \ln y_2(t)\ge ({\overline{b}}_2-\epsilon )t-\frac{f_2}{h_2-\epsilon }{\int }_0^t{y}_2(s)\mathrm{d}s+{\int }_0^t{\sigma }_2(s)\mathrm{d}{B}_2(s).\end{array}$(23) Making use of Lemma 2.3 to (22(22) $\begin{array}{rl}& \ln y_2(t)\le ({\overline{b}}_2+\epsilon )t-\frac{f_2}{h_2+\epsilon }{\int }_0^t{y}_2(s)\mathrm{d}s+{\int }_0^t{\sigma }_2(s)\mathrm{d}{B}_2(s),\end{array}$(22) ) and (23(23) $\begin{array}{rl}& \ln y_2(t)\ge ({\overline{b}}_2-\epsilon )t-\frac{f_2}{h_2-\epsilon }{\int }_0^t{y}_2(s)\mathrm{d}s+{\int }_0^t{\sigma }_2(s)\mathrm{d}{B}_2(s).\end{array}$(23) ) respectively, we obtain $\begin{array}{rl}⟨y_2(t){⟩}^{\ast }& \le \frac{(h_2+\epsilon )({\overline{b}}_2+\epsilon )}{f_2},\\ ⟨y_2(t){⟩}_{\ast }& \ge \frac{(h_2-\epsilon )({\overline{b}}_2-\epsilon )}{f_2}.\end{array}$Therefore, from the arbitrariness of ε, we can derive that $\underset{t\to +\mathrm{\infty }}{lim}⟨y_2(t)⟩=h_2{\overline{b}}_2/f_2,$ a.s.

Theorem 2.3

For model (5(5) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t).\end{array}\end{array}$(5) ), if ${\overline{b}}_1>0,{\overline{b}}_2<0$ and ${\overline{b}}_3<0$, then $y_2,y_3$ become extinct and $y_1$ is persistent in the mean, $i.e.$ $\underset{t\to +\mathrm{\infty }}{lim}⟨y_1(t)⟩=\frac{{\overline{b}}_1}{a},$ a.s.

Proof.

Since ${\overline{b}}_2<0,{\overline{b}}_3<0$, we then deduce from Theorem 2.1 that $\underset{t\to +\mathrm{\infty }}{lim}{y}_2(t)=0,\underset{t\to +\mathrm{\infty }}{lim}{y}_3(t)=0,$ a.s. Then for sufficiently large t, by (18(18) $\begin{array}{rl}\ln y_1(t)-\ln y_1(0)& ={\int }_0^t{b}_1(s)\mathrm{d}s-\frac{{\tilde{r}}_{10}-r_{10}}{{\alpha }_1}(e^{-{\alpha }_{1}t}-1)-a{\int }_0^t{y}_1(s)\mathrm{d}s\\ & -c_2{\int }_0^t\frac{y_2(s)}{h_2+y_1(s)}\mathrm{d}s-c_3{\int }_0^t\frac{y_3(s)}{h_3+y_1(s)}\mathrm{d}s\\ & +{\int }_0^t{\sigma }_1(s)\mathrm{d}{B}_1(s),\end{array}$(18) ), we know (24) $\begin{array}{rl}& \ln y_1(t)\le ({\overline{b}}_1+\epsilon )t-a{\int }_0^t{y}_1(s)\mathrm{d}s+{\int }_0^t{\sigma }_1(s)\mathrm{d}{B}_1(s),\end{array}$(24) (25) $\begin{array}{rl}& \ln y_1(t)\ge ({\overline{b}}_1-\epsilon )t-a{\int }_0^t{y}_1(s)\mathrm{d}s+{\int }_0^t{\sigma }_1(s)\mathrm{d}{B}_1(s).\end{array}$(25) Making use of Lemma 2.3 to (24(24) $\begin{array}{rl}& \ln y_1(t)\le ({\overline{b}}_1+\epsilon )t-a{\int }_0^t{y}_1(s)\mathrm{d}s+{\int }_0^t{\sigma }_1(s)\mathrm{d}{B}_1(s),\end{array}$(24) ) and (25(25) $\begin{array}{rl}& \ln y_1(t)\ge ({\overline{b}}_1-\epsilon )t-a{\int }_0^t{y}_1(s)\mathrm{d}s+{\int }_0^t{\sigma }_1(s)\mathrm{d}{B}_1(s).\end{array}$(25) ) respectively, we obtain $\begin{array}{rl}⟨y_1(t){⟩}^{\ast }& \le \frac{{\overline{b}}_1+\epsilon }{a},\\ ⟨y_1(t){⟩}_{\ast }& \ge \frac{{\overline{b}}_1-\epsilon }{a}.\end{array}$Therefore, from the arbitrariness of ε, we can derive that $\underset{t\to +\mathrm{\infty }}{lim}⟨y_1(t)⟩={\overline{b}}_1/a,$ a.s.

Theorem 2.4

For model (5(5) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t).\end{array}\end{array}$(5) ), if ${\overline{b}}_1<0,{\overline{b}}_2>max\{0,\frac{a_{23}{\overline{b}}_3{h}_3}{f_3}\}$ and ${\overline{b}}_3>max\{0,\frac{a_{32}{\overline{b}}_2{h}_2}{f_2}\}$, then $y_1$ become extinct and $y_2$, $y_3$ satisfies $\begin{array}{rl}& {\mathrm{\Delta }}_2\le ⟨y_2(t){⟩}_{\ast }\le ⟨y_2(t){⟩}^{\ast }\le \frac{{\overline{b}}_2{h}_2}{f_2},\\ & {\mathrm{\Delta }}_3\le ⟨y_3(t){⟩}_{\ast }\le ⟨y_3(t){⟩}^{\ast }\le \frac{{\overline{b}}_3{h}_3}{f_3},\end{array}$

where ${\mathrm{\Delta }}_2=\frac{{\overline{b}}_2{h}_2}{f_2}-\frac{a_{23}{\overline{b}}_3{h}_3{h}_2}{f_3{f}_2},{\mathrm{\Delta }}_3=\frac{{\overline{b}}_3{h}_3}{f_3}-\frac{a_{32}{\overline{b}}_2{h}_3{h}_2}{f_3{f}_2}.$That is to say, the populations $y_2,y_3$ are strongly persistent in the mean.

Proof.

Since ${\overline{b}}_1<0$, we then deduce from Theorem 2.1 that $\underset{t\to +\mathrm{\infty }}{lim}{y}_1(t)=0$ a.s.

By (19(19) $\begin{array}{rl}\ln y_2(t)-\ln y_2(0)& ={\int }_0^t{b}_2(s)\mathrm{d}s-\frac{{\tilde{r}}_{20}-r_{20}}{{\alpha }_2}(e^{-{\alpha }_{2}t}-1)-a_{23}{\int }_0^t{y}_3(s)\mathrm{d}s\\ & -f_2{\int }_0^t\frac{y_2(s)}{h_2+y_1(s)}\mathrm{d}s+{\int }_0^t{\sigma }_2(s)\mathrm{d}{B}_2(s),\end{array}$(19) ) and (20(20) $\begin{array}{rl}\ln y_3(t)-\ln y_3(0)& ={\int }_0^t{b}_3(s)\mathrm{d}s-\frac{{\tilde{r}}_{30}-r_{30}}{{\alpha }_3}(e^{-{\alpha }_{3}t}-1)-a_{32}{\int }_0^t{y}_2(s)\mathrm{d}s\\ & -f_3{\int }_0^t\frac{y_3(s)}{h_3+y_1(s)}\mathrm{d}s+{\int }_0^t{\sigma }_3(s)\mathrm{d}{B}_3(s).\end{array}$(20) ), we know that for sufficiently large t, (26) $\begin{array}{rl}& \frac{\ln y_2(t)}{t}\le ({\overline{b}}_2+\epsilon )-\frac{f_2}{h_2+\epsilon }⟨y_2(s)⟩+t^{-1}{\int }_0^t{\sigma }_2(s)\mathrm{d}{B}_2(s),\end{array}$(26) (27) $\begin{array}{rl}& \frac{\ln y_2(t)}{t}\ge ({\overline{b}}_2-\epsilon )-\frac{f_2}{h_2-\epsilon }⟨y_2(s)⟩-a_{23}⟨y_3(s)⟩+t^{-1}{\int }_0^t{\sigma }_2(s)\mathrm{d}{B}_2(s),\end{array}$(27) (28) $\begin{array}{rl}& \frac{\ln y_3(t)}{t}\le ({\overline{b}}_3+\epsilon )-\frac{f_3}{h_3+\epsilon }⟨y_3(s)⟩+t^{-1}{\int }_0^t{\sigma }_3(s)\mathrm{d}{B}_3(s),\end{array}$(28) (29) $\begin{array}{rl}& \frac{\ln y_3(t)}{t}\ge ({\overline{b}}_3-\epsilon )-\frac{f_3}{h_3-\epsilon }⟨y_3(s)⟩-a_{32}⟨y_2(s)⟩+t^{-1}{\int }_0^t{\sigma }_3(s)\mathrm{d}{B}_3(s).\end{array}$(29) An application of (i) in Lemma 2.3 to the above inequalities and the arbitrariness of ε, we can deduce (30) $\underset{t\to +\mathrm{\infty }}{lim sup}{t}^{-1}{\int }_0^t{y}_2(s)\mathrm{d}s\le \frac{{\overline{b}}_2{h}_2}{f_2},$(30) (31) $\underset{t\to +\mathrm{\infty }}{lim sup}{t}^{-1}{\int }_0^t{y}_3(s)\mathrm{d}s\le \frac{{\overline{b}}_3{h}_3}{f_3}.$(31) Substituting (30(30) $\underset{t\to +\mathrm{\infty }}{lim sup}{t}^{-1}{\int }_0^t{y}_2(s)\mathrm{d}s\le \frac{{\overline{b}}_2{h}_2}{f_2},$(30) , (31(31) $\underset{t\to +\mathrm{\infty }}{lim sup}{t}^{-1}{\int }_0^t{y}_3(s)\mathrm{d}s\le \frac{{\overline{b}}_3{h}_3}{f_3}.$(31) ) into (29(29) $\begin{array}{rl}& \frac{\ln y_3(t)}{t}\ge ({\overline{b}}_3-\epsilon )-\frac{f_3}{h_3-\epsilon }⟨y_3(s)⟩-a_{32}⟨y_2(s)⟩+t^{-1}{\int }_0^t{\sigma }_3(s)\mathrm{d}{B}_3(s).\end{array}$(29) ), (27(27) $\begin{array}{rl}& \frac{\ln y_2(t)}{t}\ge ({\overline{b}}_2-\epsilon )-\frac{f_2}{h_2-\epsilon }⟨y_2(s)⟩-a_{23}⟨y_3(s)⟩+t^{-1}{\int }_0^t{\sigma }_2(s)\mathrm{d}{B}_2(s),\end{array}$(27) ) respectively results in $\begin{array}{rl}& \frac{\ln y_3(t)}{t}\ge {\overline{b}}_3-a_{32}\frac{{\overline{b}}_2{h}_2}{f_2}-\epsilon -\frac{f_3}{h_3-\epsilon }⟨y_3(s)⟩+t^{-1}{\int }_0^t{\sigma }_3\mathrm{d}{B}_3,\\ & \frac{\ln y_2(t)}{t}\ge {\overline{b}}_2-a_{23}\frac{{\overline{b}}_3{h}_3}{f_3}-\epsilon -\frac{f_2}{h_2-\epsilon }⟨y_2(s)⟩+t^{-1}{\int }_0^t{\sigma }_2\mathrm{d}{B}_2.\end{array}$From Lemma 2.3(ii) we also deduce that $\begin{array}{rl}& \underset{t\to +\mathrm{\infty }}{lim inf}{t}^{-1}{\int }_0^t{y}_2(s)\mathrm{d}s\ge \frac{{\overline{b}}_2{h}_2}{f_2}-\frac{a_{23}{\overline{b}}_3{h}_2{h}_3}{f_2{f}_3},\\ & \underset{t\to +\mathrm{\infty }}{lim inf}{t}^{-1}{\int }_0^t{y}_3(s)\mathrm{d}s\ge \frac{{\overline{b}}_3{h}_3}{f_3}-\frac{a_{32}{\overline{b}}_2{h}_2{h}_3}{f_2{f}_3},\end{array}$the proof is completed.

Remark 1.1 Theorems 2.1–2.4 reveals some important and interesting impacts of the speed of reversion, the intensity of stochastic perturbations and the toxicant impulsive period on the persistence and extinction of the species. In fact, we can see that the persistence or not of the resource species $y_1$, $y_2$ and $y_3$ depends on the size of ${\overline{b}}_i$,i = 1, 2, 3.

3. Numerical simulations

Now we use the Milstein method offered in [14] to verify the theoretical results numerically (here we only provide the functions of ${\xi }_i^2$ since the functions of ${\alpha }_i$ can be proffered analogously).

In Figures 1–5, choose $r_{10}=0.5$, $r_{20}=0.5$, $r_{30}=0.3$, ${\tilde{r}}_{10}=0.3$, ${\tilde{r}}_{20}=0.25$, ${\tilde{r}}_{30}=0.2$, $r_{11}=r_{21}=r_{31}=0.9$, a = 0.25, $c_2=0.36$, $c_3=0.4$, $f_2=0.4$, $f_3=0.47$, $a_{23}=0.2$, $a_{32}=0.15$, $h_2=h_3=1$, ${\alpha }_1=0.31$, ${\alpha }_2=0.35$, ${\alpha }_3=0.42$, $\gamma =2$, $k_1=0.2$, $k_2=0.26$ $k_3=0.21$, $l_1=0.8$, $l_2=0.7$, $l_3=0.7$, g = 0.9, q = 0.5. The only difference between conditions of Figures 1–5 is that the values of ${\xi }_1^2$, ${\xi }_2^2$ and ${\xi }_3^2$ are different, see Table 1.

Figure 1. Solutions of system (3(3) $\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,3,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\end{array}$(3) ) for $r_{10}=0.5,r_{20}=0.5,r_{30}=0.3,{\tilde{r}}_{10}=0.3,{\tilde{r}}_{20}=0.25,{\tilde{r}}_{30}=0.2,r_{11}=r_{21}=r_{31}=0.9,a=0.25,c_2=0.36,c_3=0.4,f_2=0.4,f_3=0.47,a_{23}=0.2,a_{32}=0.15,h_2=h_3=1,{\alpha }_1=0.31,{\alpha }_2=0.35,{\alpha }_3=0.42,\gamma =2,k_1=0.2,k_2=0.26,k_3=0.21,l_1=0.8,l_2=0.7,l_3=0.7,g=0.9,q=0.5,$ and ${\xi }_1^2=0.55,{\xi }_2^2=0.9,{\xi }_3^2=\mathrm{0.5.}$

Figure 5. Solutions of system (3(3) $\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,3,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\end{array}$(3) ) for $r_{10}=0.5,r_{20}=0.5,r_{30}=0.3,{\tilde{r}}_{10}=0.3,{\tilde{r}}_{20}=0.25,{\tilde{r}}_{30}=0.2,r_{11}=r_{21}=r_{31}=0.9,a=0.25,c_2=0.36,c_3=0.4,f_2=0.4,f_3=0.47,a_{23}=0.2,a_{32}=0.15,h_2=h_3=1,{\alpha }_1=0.31,{\alpha }_2=0.35,{\alpha }_3=0.42,\gamma =2,k_1=0.2,k_2=0.26,k_3=0.21,l_1=0.8,l_2=0.7,l_3=0.7,g=0.9,q=0.5,$ and ${\xi }_1^2=0.55,{\xi }_2^2=0.16,{\xi }_3^2=\mathrm{0.09.}$

Table 1. The values with ${\xi }_i^2$ and ${\overline{b}}_i$(i=1,2,3).

	${\xi }_1^2$	${\xi }_2^2$	${\xi }_3^2$	${\overline{b}}_1$	${\overline{b}}_2$	${\overline{b}}_3$
1	0.55	0.9	0.5	−0.006	−0.2357	−0.0726
2	0.55	0.16	0.5	−0.006	0.2929	−0.0726
3	0.55	0.9	0.09	−0.006	−0.2357	0.1714
4	0.15	0.9	0.5	0.3165	−0.2357	−0.0726
5	0.55	0.16	0.09	−0.006	0.2929	0.1714

From Table 1, we know the following cases.

• In Figure 1, ${\overline{b}}_1=-0.006<0$, ${\overline{b}}_2=-0.2357<0$, ${\overline{b}}_3=-0.0726<0$. According to Theorem 2.1, $y_1$, $y_2$ and $y_3$ die out. Figure 1 confirms these.

• In Figure 2, ${\overline{b}}_1=-0.006<0$, ${\overline{b}}_2=0.2929>0$, ${\overline{b}}_3=-0.0726<0$. By Theorem 2.2, $y_1$ and $y_3$ die out and $\underset{t\to +\mathrm{\infty }}{lim}{t}^{-1}{\int }_0^t{y}_2(s)\mathrm{d}s=\frac{h_2{\overline{b}}_2}{f_2}=\mathrm{0.7321.}$See Figure 2. $y_1$ becomes extinct and $y_2$ is persistent in the mean, that is to say, the prey $y_1$ is not the only food source for the predator $y_2$.

• In Figure 3, ${\overline{b}}_1=-0.006<0$, ${\overline{b}}_2=-0.2357<0$, ${\overline{b}}_3=0.1714>0$. By Theorem 2.2, $y_1$ and $y_2$ die out and $\underset{t\to +\mathrm{\infty }}{lim}{t}^{-1}{\int }_0^t{y}_3(s)\mathrm{d}s=\frac{h_3{\overline{b}}_3}{f_3}=\mathrm{0.3647.}$See Figure 3. $y_1$ becomes extinct and $y_3$ is persistent in the mean, that is to say, the prey $y_1$ is not the only food source for the predator $y_3$.

• In Figure 4, ${\overline{b}}_1=0.3165>0$, ${\overline{b}}_2=-0.2357<0$, ${\overline{b}}_3=-0.0726<0$. Thus according to Theorem 2.3, $y_2$ and $y_3$ die out and $\underset{t\to +\mathrm{\infty }}{lim}{t}^{-1}{\int }_0^t{y}_1(s)\mathrm{d}s=\frac{{\overline{b}}_1}{a}=\mathrm{1.2661.}$See Figure 4. This means that when predators $y_2$ and $y_3$ are extinct, the prey $y_1$ will survive without natural enemies.

• In Figure 5, ${\overline{b}}_1=-0.006<0$, ${\overline{b}}_2=0.2929>0$, ${\overline{b}}_3=0.1714>0$, $\frac{{\overline{b}}_2{h}_2}{f_2}=0.7321$, $\frac{{\overline{b}}_3{h}_3}{f_3}=0.3647$, ${\mathrm{\Delta }}_2=0.5498$, ${\mathrm{\Delta }}_3=0.1311$. Thus according to Theorem 2.4, $y_1$ dies out and $\begin{array}{r}0.5498={\mathrm{\Delta }}_2<\underset{t\to +\mathrm{\infty }}{lim}{t}^{-1}{\int }_0^t{y}_2(s)\mathrm{d}s<\frac{{\overline{b}}_2{h}_2}{f_2}=0.7321,\\ 0.1311={\mathrm{\Delta }}_3<\underset{t\to +\mathrm{\infty }}{lim}{t}^{-1}{\int }_0^t{y}_3(s)\mathrm{d}s<\frac{{\overline{b}}_3{h}_3}{f_3}=\mathrm{0.3647.}\end{array}$Figure 5 confirms these. $y_1$ dies out, $y_2$ and $y_3$ are strongly persistent in the mean. that is to say, the prey $y_1$ is not the only food source for predators $y_2$ and $y_3$.

  Compared with Figures 1–5, we find that with the increase of intensity of stochastic perturbations ${\xi }_1^2$, species $y_1$ tends to die out while the predator $y_2$ and $y_3$ remain persistence. That is to say, the stochastic perturbations on the growth rate of $y_1$ is harmful for the persistence of $y_1$ and does not imply the extinction of $y_2$ and $y_3$. Similarly, with the increase of ${\xi }_i^2$, $y_i$,i = 2, 3, tends to die out while $y_1$ tends to be persistent. From a biological point of view, this is reasonable. ${\xi }_i^2$ is intensity of stochastic perturbations, and growth rates are sensitive to these stochastic perturbations. Hence the stochastic perturbations on the growth rate of a species is harmful for the persistence of this species. In addition, although prey $y_1$ tends to die out, predators $y_2$ and $y_3$ can choose other foods, hence the parameter ${\xi }_1^2$ does not affect the survival of predators.

  In Figures 6 and 7, we choose $r_{10}=0.2$, $r_{20}=0.1$, $r_{30}=0.2$, ${\tilde{r}}_{10}=0.1$, ${\tilde{r}}_{20}=0.05$, ${\tilde{r}}_{30}=0.13$, $r_{11}=2$, $r_{21}=1$, $r_{31}=2.1$, a = 0.5, $c_2=0.2$, $c_3=0.2$, $f_2=0.85$, $f_3=0.77$, $a_{23}=0.2$, $a_{32}=0.15$, $h_2=h_3=1$, ${\alpha }_1=0.81$, ${\alpha }_2=0.85$, ${\alpha }_3=0.79$, $k_1=0.11$, $k_2=0.12$ $k_3=0.035$, $l_1=0.6$, $l_2=0.7$, $l_3=0.5$, g = 0.9, q = 0.16, ${\xi }_1^2=0.75$, ${\xi }_2^2=0.006$, ${\xi }_3^2=0.29$. The only difference between conditions of Figures 6 and 7 is that the values of γ are different, see Table 2.

  From Table 2, we know the following cases.

• In Figure 6, we choose $\gamma =3$. Then ${\overline{b}}_1=-0.0532<0$, ${\overline{b}}_2=0.0881>0$, ${\overline{b}}_3=0.0995>0$, $\frac{{\overline{b}}_2{h}_2}{f_2}=0.1036$, $\frac{{\overline{b}}_3{h}_3}{f_3}=0.1292$, ${\mathrm{\Delta }}_2=0.0732$, and ${\mathrm{\Delta }}_3=0.1091$. By Theorem 2.4, we have $y_1$ die out and $\begin{array}{r}0.0732={\mathrm{\Delta }}_2<\underset{t\to +\mathrm{\infty }}{lim}{t}^{-1}{\int }_0^t{y}_2(s)\mathrm{d}s<\frac{{\overline{b}}_2{h}_2}{f_2}=0.1036,\\ 0.1091={\mathrm{\Delta }}_3<\underset{t\to +\mathrm{\infty }}{lim}{t}^{-1}{\int }_0^t{y}_3(s)\mathrm{d}s<\frac{{\overline{b}}_3{h}_3}{f_3}=\mathrm{0.1292.}\end{array}$Figure 6 confirms these.

• In Figure 7, choose $\gamma =0.2$. Then ${\overline{b}}_1=-0.3574<0$, ${\overline{b}}_2=-0.0541<0$ and ${\overline{b}}_3=-0.0224<0$. Thus according to Theorem 2.1, $y_1$, $y_2$ and $y_3$ die out. See Figure 7.

Figure 2. Solutions of system (3(3) $\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,3,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\end{array}$(3) ) for $r_{10}=0.5,r_{20}=0.5,r_{30}=0.3,{\tilde{r}}_{10}=0.3,{\tilde{r}}_{20}=0.25,{\tilde{r}}_{30}=0.2,r_{11}=r_{21}=r_{31}=0.9,a=0.25,c_2=0.36,c_3=0.4,f_2=0.4,f_3=0.47,a_{23}=0.2,a_{32}=0.15,h_2=h_3=1,{\alpha }_1=0.31,{\alpha }_2=0.35,{\alpha }_3=0.42,\gamma =2,k_1=0.2,k_2=0.26,k_3=0.21,l_1=0.8,l_2=0.7,l_3=0.7,g=0.9,q=0.5,$ and ${\xi }_1^2=0.55,{\xi }_2^2=0.16,{\xi }_3^2=\mathrm{0.5.}$

Figure 3. Solutions of system (3(3) $\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,3,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\end{array}$(3) ) for $r_{10}=0.5,r_{20}=0.5,r_{30}=0.3,{\tilde{r}}_{10}=0.3,{\tilde{r}}_{20}=0.25,{\tilde{r}}_{30}=0.2,r_{11}=r_{21}=r_{31}=0.9,a=0.25,c_2=0.36,c_3=0.4,f_2=0.4,f_3=0.47,a_{23}=0.2,a_{32}=0.15,h_2=h_3=1,{\alpha }_1=0.31,{\alpha }_2=0.35,{\alpha }_3=0.42,\gamma =2,k_1=0.2,k_2=0.26,k_3=0.21,l_1=0.8,l_2=0.7,l_3=0.7,g=0.9,q=0.5,$ and ${\xi }_1^2=0.55,{\xi }_2^2=0.9,{\xi }_3^2=\mathrm{0.09.}$

Figure 4. Solutions of system (3(3) $\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,3,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\end{array}$(3) ) for $r_{10}=0.5,r_{20}=0.5,r_{30}=0.3,{\tilde{r}}_{10}=0.3,{\tilde{r}}_{20}=0.25,{\tilde{r}}_{30}=0.2,r_{11}=r_{21}=r_{31}=0.9,a=0.25,c_2=0.36,c_3=0.4,f_2=0.4,f_3=0.47,a_{23}=0.2,a_{32}=0.15,h_2=h_3=1,{\alpha }_1=0.31,{\alpha }_2=0.35,{\alpha }_3=0.42,\gamma =2,k_1=0.2,k_2=0.26,k_3=0.21,l_1=0.8,l_2=0.7,l_3=0.7,g=0.9,q=0.5,$ and ${\xi }_1^2=0.15,{\xi }_2^2=0.9,{\xi }_3^2=\mathrm{0.5.}$

Figure 6. Solutions of system (3(3) $\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,3,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\end{array}$(3) ) for $r_{10}=0.2,r_{20}=0.1,r_{30}=0.2,{\tilde{r}}_{10}=0.1,{\tilde{r}}_{20}=0.05,{\tilde{r}}_{30}=0.13,r_{11}=2,r_{21}=1,r_{31}=2.1,a=0.5,c_2=0.2,c_3=0.2,f_2=0.85,f_3=0.77,a_{23}=0.2,a_{32}=0.15,h_2=h_3=1,{\alpha }_1=0.81,{\alpha }_2=0.85,{\alpha }_3=0.79,\gamma =3,k_1=0.11,k_2=0.12,k_3=0.035,l_1=0.6,l_2=0.7,l_3=0.5,g=0.9,q=0.16,$ and ${\xi }_1^2=0.75,{\xi }_2^2=0.006,{\xi }_3^2=\mathrm{0.29.}$

Figure 7. Solutions of system (3(3) $\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,3,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\end{array}$(3) ) for $r_{10}=0.2,r_{20}=0.1,r_{30}=0.2,{\tilde{r}}_{10}=0.1,{\tilde{r}}_{20}=0.05,{\tilde{r}}_{30}=0.13,r_{11}=2,r_{21}=1,r_{31}=2.1,a=0.5,c_2=0.2,c_3=0.2,f_2=0.85,f_3=0.77,a_{23}=0.2,a_{32}=0.15,h_2=h_3=1,{\alpha }_1=0.81,{\alpha }_2=0.85,{\alpha }_3=0.79,\gamma =0.2,k_1=0.11,k_2=0.12,k_3=0.035,l_1=0.6,l_2=0.7,l_3=0.5,g=0.9,q=0.16,$ and ${\xi }_1^2=0.75,{\xi }_2^2=0.006,{\xi }_3^2=\mathrm{0.29.}$

Table 2. The values with γ and ${\overline{b}}_i$(i=1,2,3).

	γ	${\overline{b}}_1$	${\overline{b}}_2$	${\overline{b}}_3$
1	3	$-0.0532$	0.0881	0.0995
2	0.2	$-0.3574$	$-0.0541$	$-0.0224$

By comparing Figures 6 with 7, one can observe that the toxicant impulsive period γ in polluted environments affects species survival. Moreover,with the decrease of γ, species tends to die out.

4. Conclusions

In this paper, we take advantage of a mean-reverting Ornstein-Uhlenbeck process to portray the random perturbations in the environment and formulate a modified Leslie–Gower Holling-type II two-predator one-prey stochastic model in polluted environments with impulsive toxicant input. We consider interspecific competition between predators and we obtain sharp sufficient conditions for persistence in the mean, strongly persistent in the mean and extinction for each species of model (3(3) $\{\begin{array}{l}{\displaystyle \mathrm{d}{y}_1(t)=y_1(t)[r_{10}+({\tilde{r}}_{10}-r_{10})e^{-{\alpha }_{1}t}-r_{11}{C}_{10}(t)-a{y}_1(t)}\\ {\displaystyle -\frac{c_2{y}_2(t)}{h_2+y_1(t)}-\frac{c_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t+{\sigma }_1(t)y_1(t)\mathrm{d}{B}_1(t),}\\ {\displaystyle \mathrm{d}{y}_2(t)=y_2[r_{20}+({\tilde{r}}_{20}-r_{20})e^{-{\alpha }_{2}t}-r_{21}{C}_{20}(t)-a_{23}{y}_3(t)-\frac{f_2{y}_2(t)}{h_2+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_2(t)y_2(t)\mathrm{d}{B}_2(t),\\ {\displaystyle \mathrm{d}{y}_3(t)=y_3[r_{30}+({\tilde{r}}_{30}-r_{30})e^{-{\alpha }_{3}t}-r_{31}{C}_{30}(t)-a_{32}{y}_2(t)-\frac{f_3{y}_3(t)}{h_3+y_1(t)}]\mathrm{d}t}\\ +{\sigma }_3(t)y_3(t)\mathrm{d}{B}_3(t),\\ {\displaystyle \frac{\mathrm{d}{C}_{i0(t)}}{\mathrm{d}t}=k_i{C}_e(t)-l_i{C}_{i0}(t),i=1,2,3,}\\ {\displaystyle \frac{\mathrm{d}{C}_e(t)}{\mathrm{d}t}=-g{C}_e(t),t\ne n\gamma ,n\in Z^{+},}\\ {\displaystyle \mathrm{\Delta }{C}_e(t)=q,t=n\gamma ,n\in Z^{+},}\end{array}$(3) ).

The theorems that we obtain have some interesting biological interpretations. By these Theorems we can see that each species is either extinct, strongly persistent in the mean or persistent in the mean, relying on the sign of ${\overline{b}}_i(i=1,2,3)$.

We note that the intensity of the perturbation ${\xi }_i^2$ and the speed of reversion ${\alpha }_i$ are two key parameters in the Ornstein-Uhlenbeck process. Obviously, $\frac{\mathrm{d}{\overline{b}}_i}{\mathrm{d}{\alpha }_i}>0,\frac{\mathrm{d}{\overline{b}}_i}{\mathrm{d}{\xi }_i^2}<0.$Therefore, we conclude that along with the increase of ${\alpha }_i$ (respectively, ${\xi }_i^2$), ${\overline{b}}_i$ will increase (respectively, decrease). That is, increasing the reversion rate ${\alpha }_i$ is advantageous for the survival of species $y_i$, while increasing intensity ${\xi }_i^2$ is harmful for them.

Some interesting topics remain to be solved. For example, it would be interesting to dissect other random noises such as the telephone noise(see [23]), the Lévy noise(see [3]) or reaction diffusion(see [2]) etc. We leave these questions for future research.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The work is supported by National Science Foundation of China [grant number 12172376] and Scientific Research Project of Tianjin Municipal Education Commission [grant number 2019KJ131].