Transmission dynamics and optimal control of a Huanglongbing model with time delay

1.   Introduction

The relationship between predator and prey has long been one of the important topics concerned by scholars. In recent years, several scholars have proposed more realistic models which should take the functional response into account ^[1,2]. As a result, Aziz-Alaoui and Okiye ^[3] proposed the famous predator-prey model with modified Leslie-Gower and Holling-type Ⅱ schemes, which is described as follows:

where $x(t)$ and $y(t)$ stand for the sizes of the prey population and the predator population respectively; $a$ represents the intraspecific competition strength; $c$ means the per capita reduction rate of the prey due to the capture of the predator; $h$ stands for the safeguard of the environment; $f$ has the like signification of $c$. 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. ^{[4,5,6,7,8,9]}).

In recent years, the world economy has grown rapidly with the development of industry and agriculture. At the same time, environmental pollution is becoming more and more serious, and even poses a threat to the survival of biological populations and human beings. For example, serious soil erosion and exhaust emissions from cars on the road are destroying the biological population structure. In order to better control and understand the effects of toxic substances on species, we must assess the population survival risk of exposure to toxic substances.

Hallam and his colleagues ^[10,11,12] have opened the door to the study of environmental toxins by publishing three papers in a row that suggest the effects of toxins on deterministic models of ecosystems. From then on, many deterministic population models with toxic effects have been proposed and studied (see e.g. ^{[13,14,15,16,17,18]}). Particularly, consider that toxins are often released into the environment in pulses of regularity. For example, pesticides and heavy metals ^[19,20]. Therefore, based on the study of deterministic population models in polluted environments with impulsive toxin inputs, several authors explored the effects of toxins on population (see e.g. ^{[21,22,23,24]}).

In particular, suppose that the living organisms absorb environmental toxicant into their bodies. $C_{10}(t)$, $C_{20}(t)$ and $C_e(t)$ denote the concentration of the toxicant in the organism of the prey species, the predator species and the environment at time $t$, respectively. Suppose that the growth rate, $r_i$, is an affine function of $C_{i0}$:

Therefore, the following model of predator and prey with modified Leslie-Gower and Holling-type Ⅱ schemes in the presence of toxins is proposed.

In fact, the rate of species growth is often disturbed by random perturbations ^[25]. In general, random perturbations in the environment can be represented by white noise ^[26,27]. Therefore, we consider the perturbations of white noise to the population growth rate with $r_{i0}\rightarrow r_{i0}+\sigma_i\dot{B}_i(t)$, we obtain the following stochastic model:

where $\sigma_i^2$, $i = 1, 2$ is the intensity of white noise; $B_1(t)$ and $B_2(t)$ are mutually independent Brownian motions defined on a complete probability space $(\Omega, \mathcal{F}, \mathcal{P})$ with a filtration $\{\mathcal{F}_t\}_{t\in R_+}$. Liu et al. ^[28] probed into several dynamical characteristics of model $(1.2)$ and offered extinct and persistent conditions for the model $(1.2)$.

Model $(1.2)$ assumes that the growth rate in the random environments is linear with respect to the Gaussian white noise

Integrating on the interval $[0, T]$ results in

Hence, the variance of the average per capita growth rate $\bar{r}_{i0}$ over an interval of length $T$ tends to $\infty$ as $T\rightarrow0$. This is insufficient to describe the actual situation. Several authors ^[29,30] have claimed that using the mean-reverting Ornstein-Uhlenbeck process is a more appropriate approach to incorporate the environment perturbations. On account of this method ^[31], one has

i.e.

where $\tilde{r}_{i0} = \hat{r}_{i0}(0)$, $\sigma_i(t) = \frac{\xi_i}{\sqrt{2\alpha_i}}\sqrt{1-e^{-2\alpha_it}}$, $\alpha_i > 0$ represents the speed of reversion, $\xi_i^2$ is the intensity of stochastic perturbations. Based on the ideas above, a three-species predator-prey model can be expressed as follows:

where $y_1(t)$ is the population size of the prey at time $t$, $y_i(t), i = 2, 3$ is the population size of the predator at time $t$.

Now let us introduce the model of the concentration of toxicant. Suppose that $C_{i0}(t)$ satisfies the following model:

where $k_i$ stands for the uptake rate of toxicant from the environment; $l_i$ denotes the loss rate of the toxicant from the species. $C_e(t)$ denotes the concentrations of the toxicant in the environment at time $t$ and satisfies the following model:

where $\bigtriangleup\zeta(t) = \zeta(t^{+})-\zeta(t)$, $h$ is the loss rate of toxicant from environment, $q$ is the toxicant input amount at every time, $\gamma$ stands for the period of the impulsive input of toxicant.

Therefore, we have the following one-prey two-predator system in polluted environments with pulse toxicant input:

Remark 1 In model $(1.4)$, 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\leq1$. 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$. In addition, $C_{i0}$ and $C_e$ stand for the concentrations of the toxicant, therefore $0\leq C_{i0}(t)\leq1$ and $0\leq C_e(t)\leq1$ for $t\geq0$. As a result, the following conditions need to be met: $k_i\leq l_i$, $q\leq1-e^{-h\gamma}$, $i = 1, 2, 3$.

To the best of our knowledge, Liu et al. ^[28] only studied the persistence and extinction of two species system in polluted environment, little research has been done on the dynamics of corresponding three species system. Therefore, we deeply analyze the properties of model $(1.4)$ in the stochastic environment improved by the Ornstein-Uhlenbeck process.

The arrangement of this paper is as follows. In section 2, the persistence and extinction threshold for each species are proposed. 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:

Lemma 1 Consider the following subsystem of model $(1.4)$ ^[32]:

The model has a unique positive $\gamma$-periodic solution $(\widetilde{C}_{10}(t), \widetilde{C}_{20}(t), \widetilde{C}_{30}(t), \widetilde{C}_{e}(t))^T$ which satisfies

Lemma 2 For arbitrary $(y_1(0), y_2(0), y_3(0))\in R_{+}^{3}$, model $(1.3)$ possesses a unique solution $(y_1(t), y_2(t), y_3(t))\in R_{+}^{3}$ for all $t\geq0$ a.s.(almost surely).

Proof. Pay attention to the following system:

and $u(0) = lny_1(0)$, $v(0) = lny_2(0)$, $w(0) = lny_3(0)$. Due to the fact that the coefficients of model $(2.1)$ satisfy the local Lipschitz condition, model $(2.1)$ possesses a unique local solution $(u(t), v(t), w(t))^T$ on $[0, \tau_*)$, where $\tau_*$ means the explosion time. From the Itô'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 $(1.3)$.

Now let us verify that $\tau_* = +\infty$. Consider the following systems:

On the basis of the comparison theorem for stochastic differential equations ^[33], we get for $t\in[0, \tau_*)$,

According to Theorem 2.2 in Jiang and Shi ^[34], we get

Due to the fact that $\phi(t), n(t), N(t), m(t)$ and $M(t)$ are global, we can know that $\tau_* = +\infty$.

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

(i) If there exist three positive constants $t_0$, $\beta$ and $\beta_0$ such that for all $t\geq t_0$,

where $F(t)/t\rightarrow0$ as $t\rightarrow +\infty$, then

$\limsup_{t\rightarrow +\infty}t^{-1}\int_0^tX(s)ds\leq\beta/\beta_0$, a.s.

(ii) If there exist three positive constants $t_0$, $\beta$ and $\beta_0$ such that for all $t\geq t_0$,

where $F(t)/t\rightarrow0$ as $t\rightarrow +\infty$, then

$\liminf_{t\rightarrow +\infty}t^{-1}\int_0^tX(s)ds\geq\beta/\beta_0$ , a.s.

Lemma 4 Suppose that $\bar{b}_1 > 0$, then

(i) If $\bar{b}_2 > 0$, then $\lim_{t\rightarrow +\infty}\frac{\ln y_2(t)}{t} = 0$, a.s.

(ii) If $\bar{b}_3 > 0$, then $\lim_{t\rightarrow +\infty}\frac{\ln y_3(t)}{t} = 0$, a.s.

Proof (ⅰ). Choose sufficiently large T which fulfils that, for $t\geq T$,

For $t\geq T$, one can deduce from $(2.8)$ that

where

Note that $L_1(t)\geq1$, consequently,

where

Then $(2.10)$ implies that

where

Therefore,

Note that $\lim_{t\rightarrow +\infty}t^{-1}\int_0^t\sigma_i(s)dB_i(s) = 0(i = 1, 2, 3)$, we then deduce that if $\bar{b}_2 > 0$,

$\lim_{t\rightarrow +\infty}t^{-1}\ln L_4(t) = 0$, a.s.

This together with $(2.13)$, indicates

Using Itô's formula to $(2.3)$ deduces

In other words,

Clearly, for arbitrary $\varepsilon > 0$, there exists $T > 0$ such that, for $t > T$,

Then, we derive from $(2.14)$ that, for $t\geq T$,

Choose $\varepsilon$ be sufficiently small such that $0\leq2\varepsilon\leq\bar{b}_2$. Applying (i) and (ii) in Lemma 3 yields that

An application of the arbitrariness of $\varepsilon$, we can obtain that

Which implies that $\lim_{t\rightarrow +\infty}t^{-1}\ln n(t) = 0$ a.s. Thanks to $(2.7)$, one can derive that

The proof of (ⅰ) is completed.

The proof of (ⅱ) is similar to that of (ⅰ), so it is omitted here.

Theorem 1 For model $(1.3)$, the following conclusions hold:

(I) If $\bar{b}_1 < 0, \bar{b}_2 < 0$ and $\bar{b}_3 < 0$, then $y_1, y_2$ and $y_3$ become extinct, $i.e.$ $\lim_{t\rightarrow +\infty}y_i(t) = 0,$ a.s., $i = 1, 2, 3.$

(Ⅱ) If $\bar{b}_1 < 0, \bar{b}_2 > 0$ and $\bar{b}_3 < 0$, then $y_1$ and $y_3$ become extinct and $y_2$ is persistent in the mean, $i.e.$ $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_2(s)ds = h_2\bar{b}_2/f_2$, a.s.

(III) If $\bar{b}_1 < 0, \bar{b}_2 < 0$ and $\bar{b}_3 > 0$, then $y_1$ and $y_2$ become extinct and $y_3$ is persistent in the mean, $i.e.$ $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_3(s)ds = h_3\bar{b}_3/f_3$, a.s.

(IV) If $\bar{b}_1 > 0, \bar{b}_2 < 0$ and $\bar{b}_3 < 0$, then $y_2$ and $y_3$ become extinct and $y_1$ is persistent in the mean, $i.e.$ $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_1(s)ds = \bar{b}_1/a$, a.s.

(V) If $\bar{b}_1 < 0, \bar{b}_2 > 0$ and $\bar{b}_3 > 0$, then $y_1$ becomes extinct and $y_2$ and $y_3$ are persistent in the mean, $i.e.$ $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_2(s)ds = h_2\bar{b}_2/f_2, \lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_3(s)ds = h_3\bar{b}_3/f_3$, a.s.

(VI) If $\bar{b}_1 > 0, \bar{b}_2 > 0$ and $\bar{b}_3 < 0$, then $y_3$ becomes extinct, and

(i) If $\frac{\bar{b}_1}{c_2} < \frac{\bar{b}_2}{f_2}$, then $y_1$ becomes extinct and $y_2$ is persistent in the mean, $i.e.$ $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_2(s)ds = h_2\bar{b}_2/f_2$, a.s.

(ii) If $\frac{\bar{b}_1}{c_2} > \frac{\bar{b}_2}{f_2}$, then $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_1(s)ds = \bar{b}_1/a-c_2\bar{b}_2/af_2,$

$\lim_{t\rightarrow +\infty}t^{-1}\int_0^t\frac{y_2(s)}{h_2+y_1(s)}ds = \bar{b}_2/f_2$, a.s.

(VII) If $\bar{b}_1 > 0, \bar{b}_2 < 0$ and $\bar{b}_3 > 0$, then $y_2$ becomes extinct, and

(i) If $\frac{\bar{b}_1}{c_3} < \frac{\bar{b}_3}{f_3}$, then $y_1$ becomes extinct and $y_3$ is persistent in the mean, $i.e.$ $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_3(s)ds = h_3\bar{b}_3/f_3$, a.s.

(ii) If $\frac{\bar{b}_1}{c_3} > \frac{\bar{b}_3}{f_3}$, then $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_1(s)ds = \bar{b}_1/a-c_3\bar{b}_3/af_3,$

$\lim_{t\rightarrow +\infty}t^{-1}\int_0^t\frac{y_3(s)}{h_3+y_1(s)}ds = \bar{b}_3/f_3$, a.s.

(VIII) If $\bar{b}_1 > 0, \bar{b}_2 > 0$ and $\bar{b}_3 > 0$, then,

(i) If $\bar{b}_1 < \frac{c_2}{f_2}\bar{b}_2+\frac{c_3}{f_3}\bar{b}_3$, then $y_1$ becomes extinct and $y_2$ and $y_3$ are persistent in the mean, $i.e.$ $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_2(s)ds = h_2\bar{b}_2/f_2,$ $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_3(s)ds = h_3\bar{b}_3/f_3$, a.s.

(ii) If $\bar{b}_1 > \frac{c_2}{f_2}\bar{b}_2+\frac{c_3}{f_3}\bar{b}_3$, then $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_1(s)ds = \bar{b}_1/a-c_2\bar{b}_2/af_2-c_3\bar{b}_3/af_3,$ $\lim_{t\rightarrow +\infty}t^{-1}\int_0^t\frac{y_2(s)}{h_2+y_1(s)}ds = \bar{b}_2/f_2,$ $\lim_{t\rightarrow +\infty}t^{-1}\int_0^t\frac{y_3(s)}{h_3+y_1(s)}ds = \bar{b}_3/f_3$, a.s.

Proof. Applying Itô's formula to model $(1.3)$ gives

As a consequence,

First, let us prove (Ⅰ), we derive from $(2.19)$ that, for sufficiently large t,

Note that $\lim_{t\rightarrow +\infty}t^{-1}\int_0^t\sigma_1(s)dB_1(s) = 0$ and $\bar{b}_1+\varepsilon < 0$, thus $\lim_{t\rightarrow +\infty}y_1(t) = 0,$ a.s. In the same way, $\bar{b}_2 < 0$ means that $\lim_{t\rightarrow +\infty}y_2(t) = 0,$ a.s., $\bar{b}_3 < 0$ means that $\lim_{t\rightarrow +\infty}y_3(t) = 0,$ a.s.

(Ⅱ). Since $\bar{b}_1 < 0,$ $\bar{b}_3 < 0$, we then deduce from (Ⅰ) that $\lim_{t\rightarrow +\infty}y_1(t) = 0, \lim_{t\rightarrow +\infty}y_3(t) = 0,$ a.s. Then for sufficiently large t,

Making use of Lemma 3 to $(2.23)$ and $(2.24)$ respectively, we obtain

Therefore, from the arbitrariness of $\varepsilon$, we can derive that $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_2(s)ds = h_2\bar{b}_2/f_2,$ a.s.

The proof of (Ⅲ) and (Ⅳ) is similar to that of (Ⅱ), thus is omitted here.

Now let us prove (Ⅴ). According to $\bar{b}_1 < 0$, we have $\lim_{t\rightarrow +\infty}y_1(t) = 0$. The proof is similar to (Ⅱ), so here is omitted.

(Ⅵ). First of all, $\bar{b}_3 < 0$ implies that $\lim_{t\rightarrow +\infty}y_3(t) = 0,$ a.s.

(ⅰ) Compute that $(2.19)\times f_2-(2.20)\times c_2$, then we get

In view of Lemma 4, for any $\varepsilon > 0$, there exists a $T > 0$ such that for $t > T$, we can see that

As a consequence, for $t > T$,

Let $\varepsilon$ be sufficiently small such that $0 < \varepsilon < \frac{c_2\bar{b}_2-\bar{b}_1f_2}{f_2+c_2+1}$. Consequently, $\lim_{t\rightarrow +\infty}y_1(t) = 0$. The proof of $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_2(s)ds = \frac{h_2\bar{b}_2}{f_2}$ is similar to that of (Ⅱ) and here is omitted.

(ⅱ) It follows from $(2.20)$ that

Making use of Lemma 4 and $\lim_{t\rightarrow +\infty}t^{-1}\int_0^t\sigma_2(s)dB_2(s) = 0$, we can obtain

Then, for any $\varepsilon > 0$, there is $T > 0$ such that for $t > T$,

Substitute $(2.29)$ into $(2.19)$, and we can get that, for $t\geq T$,

Let $\varepsilon$ be sufficiently small such that $0 < \varepsilon < (\bar{b}_1-c_2\bar{b}_2/f_2)/2$. Thanks to Lemma 3, we have

We then derive from the arbitrariness of $\varepsilon$ that $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_1(s)ds = \bar{b}_1/a-c_2\bar{b}_2/(af_2)$.

The proof of (Ⅶ) is similar to that of (Ⅵ), so it is omitted here.

Now let us proof (Ⅷ).

(ⅰ) Compute that $(2.19)\times f_2f_3-(2.20)\times f_3c_2-(2.21)\times f_2c_3$,

In view of Lemma 4, for any $\varepsilon > 0$, there exists a $T > 0$ such that for $t > T$, we can see that

As a consequence, for $t > T$,

Let $\varepsilon$ be sufficiently small such that $0 < \varepsilon < \frac{c_3f_2\bar{b}_3+c_2f_3\bar{b}_2-f_2f_3\bar{b}_1}{1+f_2f_3+c_2f_3+c_3f_2}$, we have $\lim_{t\rightarrow +\infty}y_1(t) = 0$. According to the proof of (Ⅱ), we derive that

Similar to the proof of (Ⅱ), we can obtain

(ⅱ) It follows from $(2.20)$, $(2.21)$, Lemma 4 and $\lim_{t\rightarrow +\infty}t^{-1}\int_0^t\sigma_i(s)dB_i(s) = 0, i = 1, 2, 3$ that

Then, for any $\varepsilon > 0$, there is $T > 0$ such that for $t > T$,

Substitute $(2.31)$ into $(2.19)$, we can get that, for $t\geq T$,

Let $\varepsilon$ be sufficiently small such that $0 < \varepsilon < (\bar{b}_1-c_2\bar{b}_2/f_2-c_3\bar{b}_3/f_3)/2$. Thanks to Lemma 3, we have

We then derive from the arbitrariness of $\varepsilon$ that $\lim_{t\rightarrow +\infty}t^{-1}\int_0^ty_1(s)ds = \bar{b}_1/a-c_2\bar{b}_2/(af_2)-c_3\bar{b}_3/(af_3)$.

The proof of Theorem 1 is completed.

3.   Numerical simulations

Now we use the Milstein method offered in ^[36] 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 Figure 1-4, 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.5$, $f_3 = 0.47$, $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$, $h = 0.9$, $q = 0.5$. The only difference between conditions of Figure 1-4 is that the values of $\xi_1^2$, $\xi_2^2$ and $\xi_3^2$ are different.

(Ⅰ) In Figure 1, we choose $\xi_1^2 = 0.55$, $\xi_2^2 = 0.9$ and $\xi_3^2 = 0.5$. Then $\bar{b}_1 = -0.006 < 0$, $\bar{b}_2 = -0.2357 < 0$, $\bar{b}_3 = -0.0726 < 0$. According to (Ⅰ) in Theorem 1, $y_1$, $y_2$ and $y_3$ die out.

Figure 1 confirms these.

(Ⅱ) In Figure 2, we choose $\xi_1^2 = 0.55$, $\xi_2^2 = 0.16$ and $\xi_3^2 = 0.5$. Then $\bar{b}_1 = -0.006 < 0$, $\bar{b}_2 = 0.2929 > 0$, $\bar{b}_3 = -0.0726 < 0$. On the basis of (Ⅱ) in Theorem 1, $y_1$ and $y_3$ die out and

See Figure 2.

(Ⅲ) In Figure 3, we choose $\xi_1^2 = 0.15$, $\xi_2^2 = 0.9$ and $\xi_3^2 = 0.5$. Then $\bar{b}_1 = 0.2915 > 0$, $\bar{b}_2 = -0.2357 < 0$, $\bar{b}_3 = -0.0726 < 0$. On the basis of (Ⅳ) in Theorem 1, $y_2$ and $y_3$ die out and

See Figure 3.

(Ⅳ) In Figure 4, we choose $\xi_1^2 = 0.55$, $\xi_2^2 = 0.16$ and $\xi_3^2 = 0.09$. Then $\bar{b}_1 = -0.006 < 0$, $\bar{b}_2 = \bar{b}_2 = 0.2929 > 0$, $\bar{b}_3 = 0.1714 > 0$. On the basis of (Ⅴ) in Theorem 1, $y_1$ dies out and

Figure 4 confirms these.

In Figures 5 and 6, we choose $r_{10} = 0.6$, $r_{20} = 0.8$, $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.5$, $f_3 = 0.47$, $h_2 = h_3 = 1$, $\alpha_1 = 0.81$, $\alpha_2 = 0.65$, $\alpha_3 = 0.42$, $\gamma = 2$, $k_1 = 0.2$, $k_2 = 0.23$ $k_3 = 0.45$, $l_1 = 0.4$, $l_2 = 0.8$, $l_3 = 0.7$, $h = 0.9$, $q = 0.5$. The only difference between conditions of Figures 5 and 6 is that the values of $\xi_1^2$, $\xi_2^2$ and $\xi_3^2$ are different.

(Ⅴ) In Figure 5, we choose $\xi_1^2 = 0.5$, $\xi_2^2 = 0.06$ and $\xi_3^2 = 0.5$. Then $\bar{b}_1 = 0.3207 > 0$, $\bar{b}_2 = 0.7050 > 0$, $\bar{b}_3 = -0.1583 < 0$ and $\bar{b}_1 < c_2\bar{b}_2/f_2 = 0.5076$. On the basis of (Ⅵ) in Theorem 1, $y_1$ and $y_3$ die out and

See Figure 5.

(Ⅵ) In Figure 6, we choose $\xi_1^2 = 0.5$, $\xi_2^2 = 0.86$ and $\xi_3^2 = 0.5$. Then $\bar{b}_1 = 0.3207 > 0$, $\bar{b}_2 = 0.3974 > 0$, $\bar{b}_3 = -0.1583 < 0$ and $\bar{b}_1 > c_2\bar{b}_2/f_2 = 0.2861$. On the basis of (Ⅵ) in Theorem 1, $y_3$ dies out and

See Figure 6.

In Figures 7 and 8, we choose $r_{10} = 0.6$, $r_{20} = 0.6$, $r_{30} = 0.4$, $\tilde{r}_{10} = 0.3$, $\tilde{r}_{20} = 0.15$, $\tilde{r}_{30} = 0.3$, $r_{11} = r_{21} = r_{31} = 0.1$, $a = 0.25$, $c_2 = 0.36$, $c_3 = 0.4$, $f_2 = 0.75$, $f_3 = 0.77$, $h_2 = h_3 = 1$, $\alpha_1 = 0.61$, $\alpha_2 = 0.85$, $\alpha_3 = 0.89$, $\gamma = 2$, $k_1 = 0.32$, $k_2 = 0.05$ $k_3 = 0.25$, $l_1 = 0.9$, $l_2 = 0.7$, $l_3 = 0.8$, $h = 0.9$, $q = 0.6$. The only difference between conditions of Figures 7 and 8 is that the values of $\xi_1^2$, $\xi_2^2$ and $\xi_3^2$ are different.

(Ⅶ) In Figure 7, we choose $\xi_1^2 = 0.5$, $\xi_2^2 = 0.006$ and $\xi_3^2 = 0.009$. Then $\bar{b}_1 = 0.3832 > 0$, $\bar{b}_2 = 0.5959 > 0$, $\bar{b}_3 = 0.3871 > 0$ and $\bar{b}_1 < c_2\bar{b}_2/f_2+c_3\bar{b}_3/f_3 = 0.4871$. On the basis of (Ⅷ) in Theorem 1, $y_1$ dies out and

See Figure 7.

(Ⅷ) In Figure 8, we choose $\xi_1^2 = 0.0002$, $\xi_2^2 = 0.006$ and $\xi_3^2 = 0.009$. Then $\bar{b}_1 = 0.5881 > 0$, $\bar{b}_2 = 0.5959 > 0$, $\bar{b}_3 = 0.3871 > 0$ and $\bar{b}_1 > c_2\bar{b}_2/f_2+c_3\bar{b}_3/f_3 = 0.4871$. On the basis of (Ⅷ) in Theorem 1, we have

See Figure 8.

In Figures 9 and 10, 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.6$, $a = 0.5$, $c_2 = 0.2$, $c_3 = 0.2$, $f_2 = 0.95$, $f_3 = 0.97$, $h_2 = h_3 = 1$, $\alpha_1 = 0.81$, $\alpha_2 = 0.86$, $\alpha_3 = 0.8$, $k_1 = 0.11$, $k_2 = 0.12$ $k_3 = 0.125$, $l_1 = 0.6$, $l_2 = 0.7$, $l_3 = 0.7$, $h = 0.9$, $q = 0.12$, $\xi_1^2 = 0.02$, $\xi_2^2 = 0.006$, $\xi_3^2 = 0.009$. The only difference between conditions of Figures 9 and 10 is that the values of $\gamma$ are different.

(Ⅸ) In Figure 9, we choose $\gamma = 4$. Then $\bar{b}_1 = 0.1816 > 0$, $\bar{b}_2 = 0.0925 > 0$, $\bar{b}_3 = 0.1817 > 0$ and $\bar{b}_1 > c_2\bar{b}_2/f_2+c_3\bar{b}_3/f_3 = 0.0569$. On the basis of (Ⅷ) in Theorem 1, we have

See Figure 9.

In Figure 10, choose $\gamma = 0.23$. Then $\bar{b}_1 = -0.0187 < 0$, $\bar{b}_2 = -0.0011 < 0$ and $\bar{b}_3 = -0.0720 < 0$. On the basis of (Ⅰ) in Theorem 1, $y_1$, $y_2$ and $y_3$ die out.

See Figure 10.

By comparing Figure 9 with Figure 10, one can observe that with the decrease of toxicant impulsive period $\gamma$, 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 assume that the toxicants are released in regular pulses. Based on the classical deterministic predator-prey model with modified Leslie-Gower Holling-type Ⅱ schemes, we present a three-species predator prey stochastic model with modified Leslie-Gower Holling-type Ⅱ schemes, and use more appropriate methods to describe random perturbations in the environment. We obtain sharp sufficient conditions for persistence in the mean and extinction for each species of model $(1.4)$.

Theorem 1 has some interesting biological interpretations. By Theorem 1, we can see that each species is either extinct or persistent in the mean, relying on the sign of $\bar{b}_i (i = 1, 2, 3)$, $\bar{b}_1f_2-c_2\bar{b}_2$, $\bar{b}_1f_3-c_3\bar{b}_3$, and $f_2f_3\bar{b}_1-c_2f_3\bar{b}_2-c_3f_2\bar{b}_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,

Therefore, with the increase of $\alpha_i$ (respectively, $\xi^{2}_i$), species $y_{i}$ tends to be persistent (respectively, extinct), $i = 1, 2, 3$. Furthermore, with the increase of $\alpha_2$ or $\alpha_3$ (respectively, $\xi^{2}_2$ or $\xi^{2}_3$), the prey population $y_1$ tends to die out (respectively, be persistent) provided $\bar{b}_i > 0, i = 1, 2, 3.$

Some interesting topics remain to be solved. For example, it would be interesting to dissect other random noises such as the telephone noise (see ^[37]), the L$\rm \acute{e}$vy noise (see ^[38]) or reaction diffusion (see ^[39]) etc. We leave these questions for future research.

Acknowledgments

The work is supported by National Science Foundation of China(No.11672326) and Scientific Research Project of Tianjin Municipal Education Commission (No.2019KJ131).

Conflict of interest

The authors declare that they have no conflicts of interest.