Dynamical Analysis of a Stochastic Predator-Prey System with Lévy Noise and Impulsive Toxicant Input

Xiaowan She, Yuanfu Shao, Taolin Zhang

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

DOI: 10.4236/jamp.2020.89149                

Abstract

This paper established a modified Leslie-Gower and Holling-type IV stochastic predator-prey model with Lévy noise and impulsive toxicant input. We study the stability in distribution of solutions by inequality techniques and ergodic method. By comparison method and Itô’s formula, we obtain the sufficient conditions for the survival of each species. Some numerical simulations are introduced to show the theoretical results.

Keywords

Lévy Noise, Impulsive Toxicant Input, Comparison Method, Extinction and Persistence

1. Introduction

With the development of modernization, pollution is also being produced. Air pollution, water pollution, noise pollution and other pollution affect the stability of the ecosystem. At the same time, environmental pollution affects the survival of the natural population and human life [1] [2] [3] [4] [5]. For example, the use of chemical pesticides has effectively controlled the pest problem in agriculture, but it is also widely regarded as one of the problems that have a negative impact on the environment and food safety [5]. The ecotoxicity produced by microplastics will be transferred and diffused to the entire aquatic environment, which affects the stability of the ecological environment [6]. These examples show that uncontrolled input of toxicant affects the balance of the ecosystem, and even leads to the extinction of populations. Therefore, environmental pollution will inevitably attract people’s great attention. Research on the survival of populations in polluted environments has become a hot spot [7] [8] [9] [10].

Zhang and Tan [11] considered a stochastic predator-prey system in a polluted environment with impulsive toxicant input and impulsive perturbations. They obtained a set of sufficient conditions for extinction, weak persistence in the mean and global attraction to any positive solution of the system. Lv, Meng and Wang [12] investigated an impulsive stochastic chemostat model with nonlinear perturbation in a polluted environment. They showed that both stochastic and impulsive toxicant inputs have great effects on the survival and extinction of the microorganism. Liu, Du and Deng [13] established a stochastic modified Leslie-Gower Holling-type II predator-prey model with impulsive toxicant input. They got the threshold between persistence in the mean and extinction for each population; then they concluded that the white noise is harmful to the sustainable growth of species.

Ecosystems may suffer sudden and catastrophic environmental disturbances, such as earthquakes, tsunamis, volcanoes, hurricanes or epidemics, etc. To explain these phenomena, Bao et al. [14] [15] considered a jump process into the stochastic Lotka-Volterra population systems and studied population dynamics of their systems at the first time. Zhao, Yuan and Zhang [16] established a stochastic competitive model with Lévy noise in an impulsive polluted environment. They showed that Lévy noise can significantly affect the persistence and extinction of each species. In this paper, we consider adding Lévy noise to the stochastic modified Leslie-Gower and Holling-type IV predator-prey system proposed by Xu et al. [17]. Then we get

$\{\begin{array}{l}\text{d}x\left(t\right)=x\left(t^{-}\right)\left(a_1-\eta x\left(t\right)-\frac{cy\left(t\right)}{{\theta }_1+nx\left(t\right)+x^2\left(t\right)}\right)\text{d}t\\ +{\sigma }_{1}x\left(t^{-}\right)\text{d}{B}_1\left(t\right)+x\left(t^{-}\right){\displaystyle {\int }_{\mathbb{Z}}}{\gamma }_1\left(u\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\\ \text{d}y\left(t\right)=y\left(t^{-}\right)\left(a_2-\frac{fy\left(t\right)}{{\theta }_2+x\left(t\right)}\right)\text{d}t+{\sigma }_{2}y\left(t^{-}\right)\text{d}{B}_2\left(t\right)\\ +y\left(t^{-}\right){\displaystyle {\int }_{\mathbb{Z}}}{\gamma }_2\left(u\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\end{array}$ (1)

where $x\left(t^{-}\right)$ and $y\left(t^{-}\right)$ are the left limit of prey populations $x\left(t\right)$ and predator populations $y\left(t\right)$ respectively. $a_i$ represents the intrinsic growth rate of the $i\text{th}$ population in a non-polluting environment. $\eta$ is the intensity of competition among individuals of $x\left(t\right)$. $\frac{fy\left(t\right)}{{\theta }_2+x\left(t\right)}$ is the modified Leslie-Gower term, which states that the number of predators has fallen due to the shortage of the most important food. $\frac{cy\left(t\right)}{{\theta }_1+nx\left(t\right)+x^2\left(t\right)}$ is Holling-type IV functional response, which refers to the change in the density of the prey that each predator is attached to per unit time. $a_1\mathrm{,}{a}_2\mathrm{,}\eta \mathrm{,}c\mathrm{,}f\mathrm{,}{\theta }_1\mathrm{,}{\theta }_2$, and $n$ are positive constants, where $\frac{a_1}{\eta }$ indicates the maximum endurance of the environment without predators. N is a real-valued Poisson counting measure with characteristic measure $\lambda$ on a measurable subset $\mathbb{Z}$ of $\left(\mathrm{0,}+\infty \right)$ with $\lambda \left(\mathbb{Z}\right)<+\infty$, $\tilde{N}\left(\text{d}t,\text{d}u\right)=N\left(\text{d}t,\text{d}u\right)-\lambda \left(\text{d}u\right)\text{d}t$. ${\gamma }_i\left(u\right)$ is the jump-diffusion coefficient.

According to the actual situation, we consider the impact of environmental pollution on the system (1). Let $C_1\left(t\right)$ and $C_2\left(t\right)$ be the concentration of toxicant in the prey organism and predator organism at time t respectively. Suppose that the growth rate $a_i$ is an affine function of $C_i\left(t\right)$ [13], the parameter $a_{i1}$ represents the dose response rate of the ith population to the concentration of the organismal toxicant:

$a_i\to a_i-a_{i1}{C}_i\left(t\right),i=1,2.$ (2)

Suppose $C_e\left(t\right)$ denotes the concentration of toxicant in the environment at time t. $k_i$ is the organism’s net uptake rate of environmental toxicant. $l_i$ and $m_i$ represent the net ingestion rate and the depuration rate respectively. h represents the rate of toxin loss in the environment due to evaporation or other reasons. Assuming that external toxins affect the entire predator-prey system by impulsive toxicant input, let $\tau$ and q represent the period and the amount of impulsive toxicant input each time respectively. So we can get a stochastic modified Leslie-Gower Holling-type IV predator-prey model with Lévy noise in impulsive toxicant input environments:

$\{\begin{array}{l}\text{d}x\left(t\right)=x\left(t^{-}\right)\left(a_1-a_{11}{C}_1\left(t\right)-\eta x\left(t\right)-\frac{cy\left(t\right)}{{\theta }_1+nx\left(t\right)+x^2\left(t\right)}\right)\text{d}t\\ +{\sigma }_{1}x\left(t^{-}\right)\text{d}{B}_1\left(t\right)+x\left(t^{-}\right){\displaystyle {\int }_{\mathbb{Z}}}{\gamma }_1\left(u\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\\ \text{d}y\left(t\right)=y\left(t^{-}\right)\left(a_2-a_{21}{C}_2\left(t\right)-\frac{fy\left(t\right)}{{\theta }_2+x\left(t\right)}\right)\text{d}t+{\sigma }_{2}y\left(t^{-}\right)\text{d}{B}_2\left(t\right)\\ +y\left(t^{-}\right){\displaystyle {\int }_{\mathbb{Z}}}{\gamma }_2\left(u\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\\ \frac{\text{d}{C}_1\left(t\right)}{\text{d}t}=k_1{C}_e\left(t\right)-\left(l_1+m_1\right)C_1\left(t\right),\\ \frac{\text{d}{C}_2\left(t\right)}{\text{d}t}=k_2{C}_e\left(t\right)-\left(l_2+m_2\right)C_2\left(t\right),\\ \frac{\text{d}{C}_e\left(t\right)}{\text{d}t}=-h{C}_e\left(t\right),t\ne n\tau ,\\ \begin{array}{l}\Delta x\left(t\right)=0,\Delta y\left(t\right)=0,\\ \Delta C_1\left(t\right)=\Delta C_2\left(t\right)=0,\Delta C_e\left(t\right)=q,\end{array}\}t=n\tau ,n\in N^{+}.\end{array}$ (3)

Toxicants affect the system (3) by impulsive input, and the system also contains Lévy noise. There are few studies on the impact of this type of model on system dynamics, so it is of great significance. We first turn the system (3) into an impulseless system through approximate solving methods. Then we can use the ergodic method to prove the distribution stability of the system. We also get the extinction and persistence of the population by use of the comparison theorem and some inequality techniques.

The organization of this paper is as follows. In Section 2, we provide preparations for the proof and calculation of the system (3). Section 3 discusses the stability of the distribution of the impulseless system (6). Then in Section 4, the threshold between persistence in the mean and extinction for each species is established. We introduce some numerical simulations to support the theory in Section 5. The final section concludes this paper.

2. Preliminaries

For the sake of convenience, we define the following notations:

$R_{+}^2=\left\{\zeta \in R^2|{\zeta }_i>0,i=1,2\right\};$

$\langle f\left(t\right)\rangle =t^{-1}{\displaystyle {\int }_0^t}f\left(s\right)\text{d}s,{\langle f\left(t\right)\rangle }^{*}=\underset{t\to +\infty }{\lim \sup }f\left(t\right),{\langle f\left(t\right)\rangle }_{*}=\underset{t\to +\infty }{\lim \inf }f\left(t\right);$

${\beta }_i\left(t\right)=a_i-a_{i1}{\tilde{C}}_i\left(t\right)-\frac{1}{2}{\sigma }_i^2+{\displaystyle {\int }_{\mathbb{Z}}}\left(\ln \left(1+{\gamma }_i\left(u\right)\right)-{\gamma }_i\left(u\right)\right)\lambda \text{d}u,i=1,2;$

${\bar{\beta }}_i=a_i-a_{i1}{C}_i^e-\frac{1}{2}{\sigma }_i^2+{\displaystyle {\int }_{\mathbb{Z}}}\left(\ln \left(1+{\gamma }_i\left(u\right)\right)-{\gamma }_i\left(u\right)\right)\lambda \text{d}u,i=1,2;$

$k_i\left(t\right)={\displaystyle {\int }_0^t}{\displaystyle {\int }_{\mathbb{Z}}}ln\left(1+{\gamma }_i\left(u\right)\right)\tilde{N}\left(\text{d}t,\text{d}u\right)\mathrm{,}i=\mathrm{1,2.}$

Moreover, as a standing hypothesis throughout this paper, we assume that $B_1\left(t\right)\mathrm{,}{B}_2\left(t\right)$ and N are independent. We also suppose that $1+{\gamma }_i\left(u\right)>0,u\in \mathbb{Z},i=1,2$. In order to facilitate the search when using the formula later, we suppose

Assumption 1. There is a constant $c>0$ such that

${\displaystyle {\int }_{\mathbb{Z}}}{\left[ln\left(1+{\gamma }_i\left(u\right)\right)\right]}^2\lambda \left(\text{d}u\right)<c\mathrm{,}i=\mathrm{1,2.}$

Then we put forward some necessary lemmas to prepare for the main results later.

Lemma 1 [18]. Consider the following subsystem of the system (3),

$\{\begin{array}{l}\begin{array}{l}\frac{\text{d}{C}_1\left(t\right)}{\text{d}t}=k_1{C}_e\left(t\right)-\left(l_1+m_1\right)C_1\left(t\right),\\ \frac{\text{d}{C}_2\left(t\right)}{\text{d}t}=k_2{C}_e\left(t\right)-\left(l_2+m_2\right)C_2\left(t\right),\\ \frac{\text{d}{C}_e\left(t\right)}{\text{d}t}=-h{C}_e\left(t\right),\end{array}\}t\ne n\tau \\ \Delta C_1\left(t\right)=\Delta C_2\left(t\right)=0,\Delta C_e\left(t\right)=q,t=n\tau ,n\in N^{+}.\end{array}$ (4)

Subsystem (4) has a unique $\tau$ -periodic solution ${\left({\tilde{C}}_1\left(t\right)\mathrm{,}{\tilde{C}}_2\left(t\right)\mathrm{,}{\tilde{C}}_e\left(t\right)\right)}^{\text{T}}$, which satisfies

${\langle {\tilde{C}}_i\left(t\right)\rangle }^{*}={\langle {\tilde{C}}_i\left(t\right)\rangle }_{*}=\underset{t\to +\infty }{\lim }\langle {\tilde{C}}_i\left(t\right)\rangle =\frac{k_{i}q}{h\left(l_i+m_i\right)\tau }=:C_i^e,i=1,2.$ (5)

Remark. $C_1\left(t\right)\mathrm{,}{C}_2\left(t\right)$ and $C_e\left(t\right)$ denote the concentrations of toxicant. According to their practical significance, we get $0\le C_1\left(t\right)\mathrm{,}{C}_2\left(t\right)\mathrm{,}{C}_e\left(t\right)\le 1$ must hold for all $t\ge 0$. From Lemma 1, it requires the following constraints [19]:

$k_i\le l_i+m_i,q\le 1-{\text{e}}^{-h\tau },i=1,2.$

In the following, we apply Lemma 1 to the system (3). Therefore, we only need to consider the following system

$\{\begin{array}{l}\text{d}x\left(t\right)=x\left(t^{-}\right)\left(a_1-a_{11}{\tilde{C}}_1\left(t\right)-\eta x\left(t\right)-\frac{cy\left(t\right)}{{\theta }_1+nx\left(t\right)+x^2\left(t\right)}\right)\text{d}t\\ +{\sigma }_{1}x\left(t^{-}\right)\text{d}{B}_1\left(t\right)+x\left(t^{-}\right){\displaystyle {\int }_{\mathbb{Z}}}{\gamma }_1\left(u\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\\ \text{d}y\left(t\right)=y\left(t^{-}\right)\left(a_2-a_{21}{\tilde{C}}_2\left(t\right)-\frac{fy\left(t\right)}{{\theta }_2+x\left(t\right)}\right)\text{d}t+{\sigma }_{2}y\left(t^{-}\right)\text{d}{B}_2\left(t\right)\\ +y\left(t^{-}\right){\displaystyle {\int }_{\mathbb{Z}}}{\gamma }_2\left(u\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\end{array}$ (6)

with initial value $x\left(0\right)>0,y\left(0\right)>0$.

Lemma 2 [20]. Suppose that $M\left(t\right)\mathrm{,}t\ge 0$, is a local martingale vanishing at time zero. Then

$\underset{t\to +\infty }{\lim }{\rho }_M\left(t\right)<+\infty \Rightarrow \underset{t\to +\infty }{\lim }\frac{M\left(t\right)}{t}=0a.s.,$ (7)

where

${\rho }_M\left(t\right)={\displaystyle {\int }_0^t}\frac{\text{d}\langle M,M\rangle \left(s\right)}{{\left(1+s\right)}^2},t\ge 0,$ (8)

and $\langle M,M\rangle \left(t\right)$ is Meyer’s angle bracket process.

Lemma 3 [21]. Suppose that population $Z\left(t\right)\in C\left[\Omega \times R_{+}\mathrm{,}{R}_{+}\right]$.

(i) If there exist some constants $T>0$, ${\lambda }_0>0$, $\lambda$, ${\sigma }_i$, and ${\lambda }_i$ such that, for all $t\ge T$,

$\ln Z\left(t\right)\le \lambda t-{\lambda }_0{\displaystyle {\int }_0^{t}z\left(s\right)\text{d}s}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\sigma }_i{B}_i\left(t\right)}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_i{\displaystyle {\int }_0^t{\displaystyle {\int }_{\mathbb{Z}}\ln \left(1+{\gamma }_i\left(v\right)\right)\tilde{N}\left(\text{d}s,\text{d}v\right)}}},a.s.,$

then

$(\begin{array}{ll}{\langle Z\rangle }^{*}\le \lambda /{\lambda }_{0}a.s.\hfill & \text{if}\lambda \ge 0,\hfill \\ \underset{t\to \infty }{\lim }Z\left(t\right)=0a.s.\hfill & \text{if}\lambda <0.\hfill \end{array}$

(ii) If there exist some constants $T>0$, ${\lambda }_0>0$, $\lambda$, ${\sigma }_i$, and ${\lambda }_i$ such that, for all $t\ge T$,

$\ln Z\left(t\right)\ge \lambda t-{\lambda }_0{\displaystyle {\int }_0^{t}z\left(s\right)\text{d}s}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\sigma }_i{B}_i\left(t\right)}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\lambda }_i{\displaystyle {\int }_0^t{\displaystyle {\int }_{\mathbb{Z}}\ln \left(1+{\gamma }_i\left(v\right)\right)\tilde{N}\left(\text{d}s,\text{d}v\right)}}},a.s.,$

then

${\langle Z\rangle }_{\mathrm{<mo>\; *\; </mo>}}\ge \lambda /{\lambda }_{0}a\mathrm{.}s\mathrm{..}$

Definition 1 [16] [21] [22] [23]. Let $X\left(t\right)={\left(x\left(t\right)\mathrm{,}y\left(t\right)\right)}^{\text{T}}\in R_{+}^2$ be a solution of system (6). Then

(a) the population $X\left(t\right)$ is said to go to extinction if ${\lim }_{t\to \infty }X\left(t\right)=0$ ;

(b) the population $X\left(t\right)$ is said to be stable in mean if ${\lim }_{t\to \infty }\langle X\left(t\right)\rangle =K$ a.s., where K is a constant;

(c) the population $X\left(t\right)$ is said to be stochastic strong persistence in mean if ${\langle X\left(t\right)\rangle }_{*}>0$ a.s.;

(d) the population $X\left(t\right)$ is said to be stochastic weak persistence in mean if ${\langle X\left(t\right)\rangle }^{*}>0$ a.s..

Lemma 4. For any given initial value ${\left(x\left(0\right)\mathrm{,}y\left(0\right)\right)}^{\text{T}}\in R_{+}^2$, system (6) has a unique global positive solution ${\left(x\left(t\right)\mathrm{,}y\left(t\right)\right)}^{\text{T}}$ on $t\ge 0$ and the solution will remain in $R_{+}^2$ almost surely.

Proof. To begin with, let us consider the following equations

$\{\begin{array}{l}\text{d}u\left(t\right)=\left(a_1-a_{11}{\tilde{C}}_1\left(t\right)-\frac{1}{2}{\sigma }_1^2-\eta {\text{e}}^{u\left(t\right)}-\frac{c{\text{e}}^{v\left(t\right)}}{{\theta }_1+n{\text{e}}^{u\left(t\right)}+{\text{e}}^{2u\left(t\right)}}\right)\text{d}t\\ +{\sigma }_1\text{d}{B}_1\left(t\right)+{\displaystyle {\int }_{\mathbb{Z}}}\ln \left(1+{\gamma }_1\left(u\right)\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\\ \text{d}v\left(t\right)=\left(a_2-a_{21}{\tilde{C}}_2\left(t\right)-\frac{1}{2}{\sigma }_2^2-\frac{f{\text{e}}^{v\left(t\right)}}{{\theta }_2+{\text{e}}^{u\left(t\right)}}\right)\text{d}t\\ +{\sigma }_2\text{d}{B}_2\left(t\right)+{\displaystyle {\int }_{\mathbb{Z}}}\ln \left(1+{\gamma }_2\left(u\right)\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\end{array}$ (9)

with initial value $u\left(0\right)=\ln x\left(0\right),v\left(0\right)=\ln y\left(0\right)$. Clearly, the coefficient of system (6) satisfy the local Lipschitz condition, then there is a unique local solution ${\left(u\left(t\right)\mathrm{,}v\left(t\right)\right)}^{\text{T}}$ on $\left[\mathrm{0,}{\tau }_e\right)$, where ${\tau }_e$ is the explosion time. According to Itô’s formula, ${\left(x\left(t\right)={\text{e}}^{u\left(t\right)},y\left(t\right)={\text{e}}^{v\left(t\right)}\right)}^{\text{T}}$ is the unique positive local solution to system (6) with initial value $x\left(0\right)>0,y\left(0\right)>0$. Now let us prove ${\tau }_e=+\infty ,a.s.$ to show this solution is global. Consider the following four auxiliary equations

$\begin{array}{l}\text{d}\Phi \left(t\right)=\Phi \left(t^{-}\right)\left(a_1-a_{11}{\tilde{C}}_1\left(t\right)-\eta \Phi \left(t\right)\right)\text{d}t+{\sigma }_1\Phi \left(t^{-}\right)\text{d}{B}_1\left(t\right)\\ +\Phi \left(t^{-}\right){\displaystyle {\int }_{\mathbb{Z}}}{\gamma }_1\left(u\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\Phi \left(0\right)=x\left(0\right);\end{array}$ (10)

$\begin{array}{l}\text{d}\varphi \left(t\right)=\varphi \left(t^{-}\right)\left(a_1-a_{11}{\tilde{C}}_1\left(t\right)-\eta \varphi \left(t\right)-\frac{c\Psi \left(t\right)}{{\theta }_1}\right)\text{d}t\\ +\varphi \left(t^{-}\right)\left({\sigma }_1\text{d}{B}_1\left(t\right)+{\displaystyle {\int }_{\mathbb{Z}}}{\gamma }_1\left(u\right)\tilde{N}\left(\text{d}t,\text{d}u\right)\right),\varphi \left(0\right)=x\left(0\right);\end{array}$ (11)

$\begin{array}{l}\text{d}\psi \left(t\right)=\psi \left(t^{-}\right)\left(a_2-a_{21}{\tilde{C}}_2\left(t\right)-\frac{f\psi \left(t\right)}{{\theta }_2}\right)\text{d}t+{\sigma }_2\psi \left(t^{-}\right)\text{d}{B}_2\left(t\right)\\ +\psi \left(t^{-}\right){\displaystyle {\int }_{\mathbb{Z}}}{\gamma }_2\left(u\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\psi \left(0\right)=y\left(0\right);\end{array}$ (12)

$\begin{array}{l}\text{d}\Psi \left(t\right)=\Psi \left(t^{-}\right)\left(a_2-a_{21}{\tilde{C}}_2\left(t\right)-\frac{f\Psi \left(t\right)}{{\theta }_2+\Phi \left(t\right)}\right)\text{d}t+{\sigma }_2\Psi \left(t^{-}\right)\text{d}{B}_2\left(t\right)\\ +\Psi \left(t^{-}\right){\displaystyle {\int }_{\mathbb{Z}}}{\gamma }_2\left(u\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\Psi \left(0\right)=y\left(0\right);\end{array}$ (13)

By the comparison theorem for stochastic differential equations [24], we have

$\varphi \left(t\right)\le x\left(t\right)\le \Phi \left(t\right),\psi \left(t\right)\le y\left(t\right)\le \Psi \left(t\right)a.s.,$ (14)

where $t\in \left[\mathrm{0,}{\tau }_e\right)$. According to Lemma 4.2 in [15], we have

$\Phi \left(t\right)=\frac{\exp \left\{{\displaystyle {\int }_0^t}{\beta }_1\left(s\right)\text{d}s+{\sigma }_1{B}_1\left(t\right)+k_1\left(t\right)\right\}}{\frac{1}{x\left(0\right)}+\eta {\displaystyle {\int }_0^t}\exp \left\{{\displaystyle {\int }_0^s}{\beta }_1\left(v\right)\text{d}v+{\sigma }_1{B}_1\left(s\right)+k_1\left(s\right)\right\}\text{d}s};$ (15)

$\varphi \left(t\right)=\frac{\exp \left\{{\displaystyle {\int }_0^t}\left[{\beta }_1\left(s\right)-\frac{c\Psi \left(s\right)}{{\theta }_1}\right]\text{d}s+{\sigma }_1{B}_1\left(t\right)+k_1\left(t\right)\right\}}{\frac{1}{x\left(0\right)}+\eta {\displaystyle {\int }_0^t}\exp \left\{{\displaystyle {\int }_0^s}\left[{\beta }_1\left(v\right)-\frac{c\Psi \left(v\right)}{{\theta }_1}\right]\text{d}v+{\sigma }_1{B}_1\left(s\right)+k_1\left(s\right)\right\}\text{d}s};$ (15)

$\psi \left(t\right)=\frac{\exp \left\{{\displaystyle {\int }_0^t}{\beta }_2\left(s\right)\text{d}s+{\sigma }_2{B}_2\left(t\right)+k_2\left(t\right)\right\}}{\frac{1}{y\left(0\right)}+\frac{f}{{\theta }_2}{\displaystyle {\int }_0^t}\exp \left\{{\displaystyle {\int }_0^s}{\beta }_2\left(v\right)\text{d}v+{\sigma }_2{B}_2\left(s\right)+k_2\left(s\right)\right\}\text{d}s};$ (17)

$\Psi \left(t\right)=\frac{\exp \left\{{\displaystyle {\int }_0^t}{\beta }_2\left(s\right)\text{d}s+{\sigma }_2{B}_2\left(t\right)+k_2\left(t\right)\right\}}{\frac{1}{y\left(0\right)}+{\displaystyle {\int }_0^t}\frac{f}{{\theta }_2+\Phi \left(s\right)}\exp \left\{{\displaystyle {\int }_0^s}{\beta }_2\left(v\right)\text{d}v+{\sigma }_2{B}_2\left(s\right)+k_2\left(s\right)\right\}\text{d}s}.$ (18)

Noting that $\Phi \left(t\right)>0,\varphi \left(t\right)>0,\psi \left(t\right)>0$ and $\Psi \left(t\right)>0$ are existent on $t\ge 0$, then we obtain ${\tau }_e=+\infty$ (Theorem 2.1 in [15] ).

3. Stability in Distribution

Lemma 5. Suppose that ${\left(x\left(t\right)\mathrm{,}y\left(t\right)\right)}^{\text{T}}$ is the positive solution of system (6) with any initial value ${\left(x\left(0\right)\mathrm{,}y\left(0\right)\right)}^{\text{T}}$, then for any $\iota >0$, there exists a positive constant K such that

$\underset{t\to +\infty }{\lim \sup }\mathbb{E}\left|x^{\iota }\left(t\right)\right|\le K,\underset{t\to +\infty }{\lim \sup }\mathbb{E}\left|y^{\iota }\left(t\right)\right|\le K.$

Proof. The proofs are very standard. Detailed proofs can refer to [15] [25] and hence are omitted here.

Lemma 6. If ${\bar{\beta }}_1>0,{\bar{\beta }}_2>0,{\left(n-c\right)}^2\le 4\left({\theta }_1-c{\theta }_2\right),{\bar{\beta }}_1>\frac{{\bar{\beta }}_2}{f}$ and $f>1$, then system (6) is asymptotically stable in distribution, i.e., when $t\to +\infty$, there is a unique probability measure $\nu (\cdot )$ such that the transition density $p\left(t\mathrm{,}\xi \mathrm{,}\cdot \right)$ of ${\left(x\left(t\right)\mathrm{,}y\left(t\right)\right)}^{\text{T}}$ converges weakly to $\nu (\cdot )$ with any given initial value $\xi \left(\theta \right)\in R_{+}^2$.

Proof. Let ${\left(x\left(t\mathrm{<mo>\; ;\; </mo>}\xi \right)\mathrm{,}y\left(t\mathrm{<mo>\; ;\; </mo>}\xi \right)\right)}^{\text{T}}$ and ${\left(x\left(t\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)\mathrm{,}y\left(t\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)\right)}^{\text{T}}$ be two solutions of system (6) with the same initial value $\xi \left(\theta \right)\in R_{+}^2$ and $\bar{\xi }\left(\theta \right)\in R_{+}^2$ respectively. Then we define

$V\left(t\right)=\left|lnx\left(t\mathrm{<mo>\; ;\; </mo>}\xi \right)-lnx\left(t\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)\right|+\left|lny\left(t\mathrm{<mo>\; ;\; </mo>}\xi \right)-lny\left(t\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)\right|\mathrm{.}$

According to the Itô’s formula with noise, we have

$\begin{array}{l}\text{d}\ln x\left(t;\xi \right)=\left({\beta }_1\left(t\right)-\eta x\left(t\right)-\frac{cy\left(t\right)}{{\theta }_1+nx\left(t\right)+x^2\left(t\right)}\right)\text{d}t\\ +{\sigma }_1\text{d}{B}_1\left(t\right)+{\displaystyle {\int }_{\mathbb{Z}}}\ln \left(1+{\gamma }_1\left(u\right)\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\\ \text{d}\ln x\left(t;\bar{\xi }\right)=\left({\bar{\beta }}_1\left(t\right)-\eta x\left(t\right)-\frac{cy\left(t\right)}{{\theta }_1+nx\left(t\right)+x^2\left(t\right)}\right)\text{d}t\\ +{\sigma }_1\text{d}{B}_1\left(t\right)+{\displaystyle {\int }_{\mathbb{Z}}}\ln \left(1+{\gamma }_1\left(u\right)\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\end{array}$

$\begin{array}{l}\text{d}\ln y\left(t;\xi \right)=\left({\beta }_2\left(t\right)-\frac{fy\left(t\right)}{{\theta }_2+x\left(t\right)}\right)\text{d}t+{\sigma }_2\text{d}{B}_2\left(t\right)+{\displaystyle {\int }_{\mathbb{Z}}}\ln \left(1+{\gamma }_2\left(u\right)\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\\ \text{d}\ln y\left(t;\bar{\xi }\right)=\left({\bar{\beta }}_2(\left(t\right)-\frac{fy\left(t\right)}{{\theta }_2+x\left(t\right)}\right)\text{d}t+{\sigma }_2\text{d}{B}_2\left(t\right)+{\displaystyle {\int }_{\mathbb{Z}}}\ln \left(1+{\gamma }_2\left(u\right)\right)\tilde{N}\left(\text{d}t,\text{d}u\right).\end{array}$

Hence,

$\begin{array}{c}{D}^{+}V\left(t\right)=-a_{11}\mathrm{sgn}\left(x\left(t;\xi \right)-x\left(t;\bar{\xi }\right)\right)\left({\tilde{C}}_1\left(t\right)-C_1^e\right)\text{d}t\\ -\eta \mathrm{sgn}\left(x\left(t;\xi \right)-x\left(t;\bar{\xi }\right)\right)\left(x\left(t;\xi \right)-x\left(t;\bar{\xi }\right)\right)\text{d}t\\ -c\mathrm{sgn}\left(x\left(t;\xi \right)-x\left(t;\bar{\xi }\right)\right)\\ \times \left(\frac{y\left(t;\xi \right)}{{\theta }_1+nx\left(t;\xi \right)+x^2\left(t;\xi \right)}-\frac{y\left(t;\bar{\xi }\right)}{{\theta }_1+nx\left(t;\bar{\xi }\right)+x^2\left(t;\bar{\xi }\right)}\right)\text{d}t\\ -a_{21}\mathrm{sgn}\left(y\left(t;\xi \right)-y\left(t;\bar{\xi }\right)\right)\left({\tilde{C}}_2\left(t\right)-C_2^e\right)\text{d}t\end{array}$

$\begin{array}{l}-f\mathrm{sgn}\left(y\left(t;\xi \right)-y\left(t;\bar{\xi }\right)\right)\left(\frac{y\left(t;\xi \right)}{{\theta }_2+x\left(t;\xi \right)}-\frac{y\left(t;\bar{\xi }\right)}{{\theta }_2+x\left(t;\bar{\xi }\right)}\right)\text{d}t\\ \le -a_{11}\mathrm{sgn}\left(x\left(t;\xi \right)-x\left(t;\bar{\xi }\right)\right)\left({\tilde{C}}_1\left(t\right)-C_1^e\right)\text{d}t\\ -\eta \mathrm{sgn}\left(x\left(t;\xi \right)-x\left(t;\bar{\xi }\right)\right)\left(x\left(t;\xi \right)-x\left(t;\bar{\xi }\right)\right)\text{d}t\\ +c\mathrm{sgn}\left(x\left(t;\xi \right)-x\left(t;\bar{\xi }\right)\right)\\ \times \left(\frac{y\left(t;\xi \right)}{{\theta }_1+nx\left(t\mathrm{<mo>\; ;\; </mo>}\xi \right)+x^2\left(t\mathrm{<mo>\; ;\; </mo>}\xi \right)}-\frac{y\left(t\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)}{{\theta }_1+nx\left(t\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)+x^2\left(t\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)}\right)\text{d}t\end{array}$

$\begin{array}{l}-a_{21}\mathrm{sgn}\left(y\left(t;\xi \right)-y\left(t;\bar{\xi }\right)\right)\left({\tilde{C}}_2\left(t\right)-C_2^e\right)\text{d}t\\ -cf\mathrm{sgn}\left(y\left(t;\xi \right)-y\left(t;\bar{\xi }\right)\right)\\ \times \left(\frac{y\left(t;\xi \right)}{{\theta }_1+nx\left(t;\xi \right)+x^2\left(t;\xi \right)}-\frac{y\left(t\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)}{{\theta }_1+nx\left(t\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)+x^2\left(t\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)}\right)\text{d}t\\ \le -{\displaystyle \underset{i\mathrm{<mo>\; =\; </mo>1}}{\overset{2}{\sum }}}{a}_{i1}\left|{\tilde{C}}_i\left(t\right)-C_i^e\right|\text{d}t-\eta \left|x\left(t\mathrm{<mo>\; ;\; </mo>}\xi \right)-x\left(t\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)\right|\text{d}t\\ -\frac{c\left(f-1\right)}{{\theta }_1}\left|y\left(t;\xi \right)-y\left(t;\bar{\xi }\right)\right|\text{d}t\mathrm{.}\end{array}$

Therefore,

$\begin{array}{c}\mathbb{E}\left(V\left(t\right)\right)\le V\left(0\right)-{\displaystyle \underset{i=1}{\overset{2}{\sum }}}{a}_{i1}{\displaystyle {\int }_0^t}\mathbb{E}\left|C_i\left(s\right)-C_i^e\right|\text{d}s-\eta {\displaystyle {\int }_0^t}\mathbb{E}\left|x\left(s;\xi \right)-x\left(s;\bar{\xi }\right)\right|\text{d}s\\ -\frac{c\left(f-1\right)}{{\theta }_1}{\displaystyle {\int }_0^t}\mathbb{E}\left|y\left(s\mathrm{<mo>\; ;\; </mo>}\xi \right)-y\left(s\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)\right|\text{d}s\mathrm{.}\end{array}$

Note that $V\left(t\right)\ge 0$, we have

$\begin{array}{l}{\displaystyle \underset{i=1}{\overset{2}{\sum }}}{a}_{i1}{\displaystyle {\int }_0^t}\mathbb{E}\left|{\tilde{C}}_i\left(s\right)-C_i^e\right|\text{d}s+\eta {\displaystyle {\int }_0^t}\mathbb{E}\left|x\left(s\mathrm{<mo>\; ;\; </mo>}\xi \right)-x\left(s\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)\right|\text{d}s\\ +\frac{c\left(f-1\right)}{{\theta }_1}{\displaystyle {\int }_0^t}\mathbb{E}\left|y\left(s\mathrm{<mo>\; ;\; </mo>}\xi \right)-y\left(s\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)\right|\text{d}s\le V\left(0\right)<+\infty \mathrm{.}\end{array}$

Hence,

$\mathbb{E}\left|x\left(t\mathrm{<mo>\; ;\; </mo>}\xi \right)-x\left(t\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)\right|\in L^1\left[\mathrm{0,}+\infty \right)\mathrm{,}$

$\mathbb{E}\left|y\left(t\mathrm{<mo>\; ;\; </mo>}\xi \right)-y\left(t\mathrm{<mo>\; ;\; </mo>}\bar{\xi }\right)\right|\in L^1\left[\mathrm{0,}+\infty \right)\mathrm{.}$

According to the first equation of system (6), we get

$\mathbb{E}\left(x\left(t\right)\right)\le x\left(0\right)+{\displaystyle {\int }_0^t}\left[\left(a_1-a_{11}{\tilde{C}}_1\left(t\right)\right)\mathbb{E}\left(x\left(s\right)\right)-\eta \mathbb{E}\left(x^2\left(s\right)\right)+\frac{c}{{\theta }_1}\mathbb{E}\left(x\left(s\right)y\left(s\right)\right)\right]\text{d}s\mathrm{.}$

Therefore, $\mathbb{E}\left(x\left(t\right)\right)$ is continuously differentiable function. By Lemma 5, it gives that

$\frac{\text{d}\mathbb{E}\left(x\left(t\right)\right)}{\text{d}t}\le a_1\mathbb{E}\left(x\left(t\right)\right)+\frac{c}{2{\theta }_1}\mathbb{E}\left(x^2\left(t\right)+y^2\left(t\right)\right)\le K_1\mathrm{,}$ (19)

where $K_1$ is a positive constant. Similarly, by the second equation of system (6), it can

$\frac{\text{d}\mathbb{E}\left(y\left(t\right)\right)}{\text{d}t}\le a_2\mathbb{E}\left(y\left(t\right)\right)\le K_2\mathrm{,}$ (20)

where $K_2$ is a positive constant.

Based on Lemma 5, we can get that $\mathbb{E}\left(x\left(t\right)\right)$ and $\mathbb{E}\left(y\left(t\right)\right)$ are uniformly continuous function through (19) and (20). By the Barbalat’s conclusion of [26], we can observe that

$\underset{t\to +\infty }{\lim }\mathbb{E}\left|x\left(t;\xi \right)-x\left(t;\bar{\xi }\right)\right|=0a.s.\text{and}\underset{t\to +\infty }{\lim }\mathbb{E}\left|y\left(t;\xi \right)-y\left(t;\bar{\xi }\right)\right|=0a.s.$ (21)

Suppose that $p\left(t\mathrm{,}\xi \mathrm{,}\text{d}z\right)$ represents the transition probability density of the process ${\left(x\left(t\right)\mathrm{,}y\left(t\right)\right)}^{\text{T}}$ and $P\left(t\mathrm{,}\xi \mathrm{,}\mathcal{A}\right)$ denotes the probability of $z\left(t\mathrm{,}\xi \right)\in \mathcal{A}$ with initial value $\xi \left(\theta \right)\in R_{+}^2$. By Chebyshev’s inequality [27] and Lemma 5, the family of $p\left(t\mathrm{,}\xi \mathrm{,}\text{d}z\right)$ is tight. So we can obtained that a compact subset $\mathcal{B}\in R_{+}^2$ such that $P\left(t\mathrm{,}\xi \mathrm{,}\mathcal{B}\right)\ge 1-\epsilon$ for any given $\epsilon >0$.

Let $\mathcal{P}\left(R_{+}^2\right)$ be the probability measures on $R_{+}^2$. For any given two measures $P_1\mathrm{,}{P}_2\in \mathcal{P}$, we define the metric

$d_L\left(P_1\mathrm{,}{P}_2\right)=\underset{g\in L}{sup}\left|{\displaystyle {\int }_{R_{+}^2}}g\left(s\right)P_1\left(\text{d}s\right)-{\displaystyle {\int }_{R_{+}^2}}g\left(s\right)P_2\left(\text{d}s\right)\right|\mathrm{,}$

where

$L=\left\{g\mathrm{<mo>\; :\; </mo>}{R}_{+}^2\to R|\left|g\left(s_1\right)-g\left(s_2\right)\right|\le \left|s_1-s_2\right|\mathrm{,}\left|g(\cdot )\le 1\right|\right\}\mathrm{.}$

For any $g\in L$ and $t,s>0$, we can get

$\begin{array}{l}\left|\mathbb{E}g\left(x\left(t+s;\xi \right)\right)-\mathbb{E}g\left(x\left(t;\xi \right)\right)\right|\\ =\left|\mathbb{E}\left[\mathbb{E}\left(g\left(x\left(t+s;\xi \right)\right)|{\mathcal{F}}_s\right)\right]-\mathbb{E}g\left(x\left(t;\xi \right)\right)\right|\\ =\left|{\displaystyle {\int }_{R_{+}^2}}\mathbb{E}g\left(x\left(t;z\right)\right)p\left(s,\xi ,\text{d}z\right)-\mathbb{E}g\left(x\left(t;\xi \right)\right)\right|\\ \le {\displaystyle {\int }_{R_{+}^2}}\left|\mathbb{E}g\left(x\left(t;z\right)\right)-\mathbb{E}g\left(x\left(t;\xi \right)\right)\right|p\left(s,\xi ,\text{d}z\right)\\ \le 2p\left(s,\xi ,U_K^c\right)+{\displaystyle {\int }_{U_K}}\left|\mathbb{E}g\left(x\left(t;z\right)\right)-\mathbb{E}g\left(x\left(t;\xi \right)\right)\right|p\left(s,\xi ,\text{d}z\right),\end{array}$

where $U_K=\left\{s\in R_{+}^2\mathrm{<mo>\; :\; </mo>}\left|s\right|\le K\right\}$, $U_K^c$ is a complementary set of $U_K$. Because the family of $p\left(t\mathrm{,}\xi \mathrm{,}\text{d}z\right)$ is tight, so there exists a sufficiently large K such that $p\left(s,\xi ,U_K^c\right)<\frac{\epsilon }{4}$ for any given $s\ge 0$. By (21), there exists a $T>0$ such that

$\underset{g\in L}{sup}\left|\mathbb{E}g\left(x\left(t\mathrm{<mo>\; ;\; </mo>}z\right)\right)-\mathbb{E}g\left(x\left(t\mathrm{<mo>\; ;\; </mo>}\xi \right)\right)\right|\le \frac{\epsilon }{2}\mathrm{,}$

holds for $t\ge T$, which yields $\left|\mathbb{E}g\left(x\left(t+s\mathrm{<mo>\; ;\; </mo>}\xi \right)\right)-\mathbb{E}g\left(x\left(t\mathrm{<mo>\; ;\; </mo>}\xi \right)\right)\right|\le \epsilon$.

It follows from the arbitrariness of $\epsilon$ that for $t\ge T$, we have

$\underset{g\in L}{sup}\left|\mathbb{E}g\left(x\left(t+s\mathrm{<mo>\; ;\; </mo>}\xi \right)\right)-\mathbb{E}g\left(x\left(t\mathrm{<mo>\; ;\; </mo>}\xi \right)\right)\right|\le \epsilon \mathrm{.}$

Similarly, we get

$\underset{g\in L}{sup}\left|\mathbb{E}g\left(y\left(t+s\mathrm{<mo>\; ;\; </mo>}\xi \right)\right)-\mathbb{E}g\left(y\left(t\mathrm{<mo>\; ;\; </mo>}\xi \right)\right)\right|\le \epsilon \mathrm{.}$

Hence,

$d_L\left(p\left(t+s\mathrm{,}\xi \mathrm{,}\cdot \right)\mathrm{,}p\left(t\mathrm{,}\xi \mathrm{,}\cdot \right)\right)\le \epsilon \mathrm{,}$

holds for $\forall t\ge T$ and $\forall s>0$.

That is to say, for any $\xi \left(\theta \right)\in R_{+}^2$, the transition probability $\left\{p\left(t,\xi ,\cdot \right):t>0\right\}$ is Cauchy in $\mathcal{P}$ with metric $d_L$. So $\left\{p\left(t,0,\cdot \right):t>0\right\}$ is Cauchy in $\mathcal{P}$ with metric $d_L$. In a word, there has a unique $\nu (\cdot )\in \mathcal{P}\left(R_{+}^2\right)$ such that $\underset{t\to 0}{lim}{d}_L\left(p\left(t\mathrm{,0,}\cdot \right)\mathrm{,}\nu (\cdot )\right)=0$. From (21), we get

$\underset{t\to 0}{lim}{d}_L\left(p\left(t\mathrm{,}\xi \mathrm{,}\cdot \right)\mathrm{,}p\left(t\mathrm{,0,}\cdot \right)\right)=0.$

Using triangle inequalities, we have

$\begin{array}{l}\underset{t\to 0}{lim}{d}_L\left(p\left(t\mathrm{,}\xi \mathrm{,}\cdot \right)\mathrm{,}\nu (\cdot )\right)\\ \le \underset{t\to 0}{lim}{d}_L\left(p\left(t\mathrm{,}\xi \mathrm{,}\cdot \right)\mathrm{,}p\left(t\mathrm{,0,}\cdot \right)\right)+\underset{t\to 0}{lim}{d}_L\left(p\left(t\mathrm{,0,}\cdot \right)\mathrm{,}\nu (\cdot )\right)=0.\end{array}$

This completes the proof of Lemma 6.

4. Extinction and Persistence

Lemma 7. If ${\bar{\beta }}_1>0$ and ${\bar{\beta }}_2>0$, then

$\underset{t\to +\infty }{lim}{t}^{-1}lny\left(t\right)=0a\mathrm{.}s\mathrm{.}$ (22)

Proof. By Lemma 1, we have

$\underset{t\to +\infty }{\lim }{t}^{-1}{\displaystyle {\int }_0^t}{\beta }_i\left(s\right)\text{d}s={\bar{\beta }}_i,i=1,2.$ (23)

Then for arbitrary $0<\epsilon <{\bar{\beta }}_i$, there is a T such that for $t\ge T$,

$\left({\bar{\beta }}_i-\epsilon \right)t\le {\displaystyle {\int }_0^t}{\beta }_i\left(s\right)\text{d}s\le \left({\bar{\beta }}_i+\epsilon \right)t\mathrm{,}i=\mathrm{1,2.}$ (24)

Set $T_1=\frac{l{n}_2}{{\bar{\beta }}_1-\epsilon }+T$, therefore for $t\ge T_1$,

$\frac{1}{2}exp\left\{\left({\bar{\beta }}_1-\epsilon \right)t\right\}\ge exp\left\{\left({\bar{\beta }}_1-\epsilon \right)T\right\}\mathrm{.}$

Consequently, for $t\ge T_1$, by (10), we have

$\begin{array}{l}\Phi \left(t\right)=\frac{\exp \left\{{\displaystyle {\int }_0^t}{\beta }_1\left(s\right)\text{d}s+{\sigma }_1{B}_1\left(t\right)+k_1\left(t\right)\right\}}{\frac{1}{x\left(0\right)}+\eta {\displaystyle {\int }_0^t}\exp \left\{{\displaystyle {\int }_0^s}{\beta }_1\left(v\right)\text{d}v+{\sigma }_1{B}_1\left(s\right)+k_1\left(s\right)\right\}\text{d}s}\\ \le \frac{\exp \left\{{\displaystyle {\int }_0^t}{\beta }_1\left(s\right)\text{d}s+{\sigma }_1{B}_1\left(t\right)+k_1\left(t\right)\right\}}{\eta {\displaystyle {\int }_0^t}\exp \left\{{\displaystyle {\int }_0^s}{\beta }_1\left(v\right)\text{d}v+{\sigma }_1{B}_1\left(s\right)+k_1\left(s\right)\right\}\text{d}s}\\ \le \frac{\exp \left\{{\displaystyle {\int }_0^t}{\beta }_1\left(s\right)\text{d}s+{\sigma }_1{B}_1\left(t\right)+k_1\left(t\right)\right\}}{\eta \exp \left\{\underset{0\le s\le t}{\min }{\sigma }_1{B}_1\left(s\right)+\underset{0\le s\le t}{\min }{k}_1\left(s\right)\right\}{\displaystyle {\int }_0^t}\exp \left\{{\displaystyle {\int }_0^s}{\beta }_1\left(v\right)\text{d}v\right\}\text{d}s}\\ \le \frac{\exp \left\{\left({\bar{\beta }}_1+\epsilon \right)t+{\sigma }_1{B}_1\left(t\right)+k_1\left(t\right)\right\}}{\eta \exp \left\{\underset{0\le s\le t}{\min }{\sigma }_1{B}_1\left(s\right)+\underset{0\le s\le t}{\min }{k}_1\left(s\right)\right\}{\displaystyle {\int }_{T_1}^t}\exp \left\{\left({\bar{\beta }}_1-\epsilon \right)s\right\}\text{d}s}\end{array}$

$\begin{array}{l}=\frac{{\bar{\beta }}_1-\epsilon }{\eta }\frac{\exp \left\{\left({\bar{\beta }}_1+\epsilon \right)t+{\sigma }_1{B}_1\left(t\right)+k_1\left(t\right)\right\}}{\exp \left\{\underset{0\le s\le t}{\min }{\sigma }_1{B}_1\left(s\right)+\underset{0\le s\le t}{\min }{k}_1\left(s\right)\right\}\left\{\exp \left[\left({\bar{\beta }}_1-\epsilon \right)t\right]-\exp \left[\left({\bar{\beta }}_1-\epsilon \right)T_1\right]\right\}}\\ \le \frac{2\left({\bar{\beta }}_1-\epsilon \right)}{\eta }\frac{\exp \left\{\left({\bar{\beta }}_1+\epsilon \right)t+{\sigma }_1{B}_1\left(t\right)+k_1\left(t\right)\right\}}{\exp \left\{\underset{0\le s\le t}{\min }{\sigma }_1{B}_1\left(s\right)+\underset{0\le s\le t}{\min }{k}_1\left(s\right)\right\}\exp \left\{\left({\bar{\beta }}_1-\epsilon \right)t\right\}}\\ \le \frac{2\left({\bar{\beta }}_1-\epsilon \right)}{\eta }\exp \left(2\epsilon t\right)\exp \left\{\left|{\sigma }_1\right|\left(B_1\left(t\right)-\underset{0\le s\le t}{\min }{B}_1\left(s\right)\right)+k_1\left(t\right)-\underset{0\le s\le t}{\min }{k}_1\left(s\right)\right\}.\end{array}$

Clearly,

$\exp \left(2\epsilon t\right)\exp \left\{\left|{\sigma }_1\right|\left(B_1\left(t\right)-\underset{0\le s\le t}{\min }{B}_1\left(s\right)\right)+k_1\left(t\right)-\underset{0\le s\le t}{\min }{k}_1\left(s\right)\right\}>1.$

Therefore,

$\begin{array}{l}{\displaystyle {\int }_{T_1}^t}\frac{f}{{\theta }_2+\Phi \left(s\right)}\exp \left\{{\displaystyle {\int }_0^s}{\beta }_2\left(v\right)\text{d}v+{\sigma }_2{B}_2\left(s\right)+k_2\left(s\right)\right\}\text{d}s\\ \ge {\displaystyle {\int }_{T_1}^t}\frac{f\exp \left\{\left({\bar{\beta }}_2-\epsilon \right)s+{\sigma }_2{B}_2\left(s\right)+k_2\left(s\right)\right\}}{{\theta }_2+\frac{2\left({\bar{\beta }}_1-\epsilon \right)}{\eta }\exp \left(2\epsilon s\right)\exp \left\{\left|{\sigma }_1\right|\left(B_1\left(s\right)-\underset{0\le u\le s}{\min }{B}_1\left(u\right)\right)+k_1\left(s\right)-\underset{0\le u\le s}{\min }{k}_1\left(u\right)\right\}}\text{d}s\\ \ge {\displaystyle {\int }_{T_1}^t}\frac{f\exp \left\{\left({\bar{\beta }}_2-\epsilon \right)s+{\sigma }_2{B}_2\left(s\right)+k_2\left(s\right)\right\}}{\left[{\theta }_2+\frac{2\left({\bar{\beta }}_1-\epsilon \right)}{\eta }\right]\exp \left(2\epsilon s\right)\exp \left\{\left|{\sigma }_1\right|\left(B_1\left(s\right)-\underset{0\le u\le s}{\min }{B}_1\left(u\right)\right)+k_1\left(s\right)-\underset{0\le u\le s}{\min }{k}_1\left(u\right)\right\}}\text{d}s\\ =\frac{\eta f}{\eta {\theta }_2+2\left({\bar{\beta }}_1-\epsilon \right)}{\displaystyle {\int }_{T_1}^t}\exp \left\{-\left|{\sigma }_1\right|\left(B_1\left(s\right)-\underset{0\le u\le s}{\min }{B}_1\left(u\right)\right)-k_1\left(s\right)+\underset{0\le u\le s}{\min }{k}_1\left(u\right)\right\}\\ \times \exp \left\{\left({\bar{\beta }}_2-3\epsilon \right)s+{\sigma }_2{B}_2\left(s\right)+k_2\left(s\right)\right\}\text{d}s\end{array}$

$\begin{array}{l}\ge \frac{\eta f}{\eta {\theta }_2+2\left({\bar{\beta }}_1-\epsilon \right)}\times \exp \left\{\left|{\sigma }_1\right|\underset{0\le s\le t}{\min }{B}_1\left(s\right)-\left|{\sigma }_1\right|\underset{T_1\le s\le t}{\max }{B}_1\left(s\right)+\underset{0\le s\le t}{\min }{k}_1\left(s\right)-\underset{T_1\le s\le t}{\max }{k}_1\left(s\right)\right\}\\ \times \exp \left\{\underset{T_1\le s\le t}{\min }{\sigma }_2{B}_2\left(s\right)+\underset{T_1\le s\le t}{\min }{k}_2\left(t\right)\right\}{\displaystyle {\int }_{T_1}^t}\exp \left\{\left({\bar{\beta }}_2-3\epsilon \right)s\right\}\text{d}s\end{array}$

$=:R_1\left(t\right)\left\{\exp \left[\left({\bar{\beta }}_2-3\epsilon \right)t\right]-\exp \left[\left({\bar{\beta }}_2-3\epsilon \right)T_1\right]\right\},$

where

$\begin{array}{c}{R}_1\left(t\right)=\frac{\eta f}{\left[\eta {\theta }_2+2\left({\bar{\beta }}_1-\epsilon \right)\right]\left({\bar{\beta }}_2-3\epsilon \right)}\\ \times \exp \left\{\left|{\sigma }_1\right|\underset{0\le s\le t}{\min }{B}_1\left(s\right)-\left|{\sigma }_1\right|\underset{T_1\le s\le t}{\max }{B}_1\left(s\right)+\underset{0\le s\le t}{\min }{k}_1\left(s\right)-\underset{T_1\le s\le t}{\max }{k}_1\left(s\right)\right\}\\ \times \exp \left\{\underset{T_1\le s\le t}{\min }{\sigma }_2{B}_2\left(s\right)+\underset{T_1\le s\le t}{\min }{k}_2\left(t\right)\right\}.\end{array}$

Substituting the above inequality into (18), we can get

$\begin{array}{c}\frac{1}{\Psi \left(t\right)}\ge \exp \left\{-{\displaystyle {\int }_{T_1}^t}{\beta }_2\left(s\right)\text{d}s-{\sigma }_2\left(B_2\left(t\right)-B_2\left(T_1\right)\right)-\left(k_2\left(t\right)-k_2\left(T_1\right)\right)\right\}\\ \times \left[\frac{1}{y\left(0\right)}+R_1\left(t\right)\left\{\exp \left[\left({\bar{\beta }}_2-3\epsilon \right)t\right]-\exp \left[\left({\bar{\beta }}_2-3\epsilon \right)T_1\right]\right\}\right]\\ \ge \exp \left\{\text{}-\left({\bar{\beta }}_2+\epsilon \right)t+B_2\left(T_1\right)+k_2\left(T_1\right)\right\}\exp \left[\left({\bar{\beta }}_2-3\epsilon \right)t\right]\\ \times R_1\left(t\right)\left(1-\exp \left\{-\left({\bar{\beta }}_2-3\epsilon \right)\left(t-T_1\right)\right\}\right)\exp \left[-\underset{T_1\le s\le t}{\max }{\sigma }_2{B}_2\left(s\right)-\underset{T_1\le s\le t}{\max }{k}_2\left(s\right)\right]\\ =:R_2\left(t\right)\times R_3\left(t\right),\end{array}$

where

$R_2\left(t\right)=\exp \left\{-4\epsilon t+B_2\left(T_1\right)+k_2\left(T_1\right)\right\}\left(1-\exp \left\{-\left({\bar{\beta }}_2-3\epsilon \right)\left(t-T_1\right)\right\}\right),$

$R_3\left(t\right)=R_1\left(t\right)\exp \left[-\underset{T_1\le s\le t}{\max }{\sigma }_2{B}_2\left(s\right)-\underset{T_1\le s\le t}{\max }{k}_2\left(s\right)\right].$

Thus

$t^{-1}\ln \Psi \left(t\right)\le -t^{-1}\ln R_2\left(t\right)-t^{-1}\ln R_3\left(t\right).$ (25)

Note that ${\lim }_{t\to +\infty }{B}_i\left(t\right)/t=0a.s.,i=1,2$, under Assumptions 1,

$\begin{array}{c}\langle k_i\left(t\right),k_i\left(t\right)\rangle \left(t\right)={\displaystyle {\int }_0^t}{\displaystyle {\int }_{\left(t\right)}}{\left[\ln \left(1+{\gamma }_i\left(u\right)\right)\right]}^2\lambda \left(\text{d}u\right)\text{d}s\\ \le t{\displaystyle {\int }_{\mathbb{Z}}}{\left[\ln \left(1+{\gamma }_i\left(u\right)\right)\right]}^2\lambda \left(\text{d}u\right)<tc,i=1,2.\end{array}$

Applying Lemma 2, we can obtain that

$\underset{t\to +\infty }{\lim }{k}_i\left(t\right)/t=0a.s.,i=1,2.$ (26)

Then, it follows from ${\bar{\beta }}_2>0$ and for arbitrary $\epsilon >0$ that

$\underset{t\to +\infty }{\lim }{t}^{-1}\ln R_2\left(t\right)=0,\underset{t\to +\infty }{\lim }{t}^{-1}\ln R_3\left(t\right)=0a.s\mathrm{..}$

By substituting the above identities into (25) results in

$\underset{t\to +\infty }{\lim \sup }{t}^{-1}\ln y\left(t\right)\le \underset{t\to +\infty }{\lim \sup }{t}^{-1}\ln \Psi \left(t\right)\le 0a.s\mathrm{..}$

And then let us prove $\underset{t\to +\infty }{\lim \inf }{t}^{-1}\ln y\left(t\right)\ge 0a.s.$. Applying Itô’s formula to (12) gives

$\begin{array}{c}\text{d}\ln \psi \left(t\right)=\left[a_2-a_{21}{\tilde{C}}_2\left(t\right)-\frac{1}{2}{\sigma }_2^2-\frac{f\psi \left(t\right)}{{\theta }_2}\right]\text{d}t\\ +{\displaystyle {\int }_{\mathbb{Z}}}\left[\ln \left(\psi \left(t^{-}\right)+{\gamma }_2\left(u\right)\psi \left(t^{-}\right)\right)-\ln \psi \left(t^{-}\right)\right]\tilde{N}\left(\text{d}t,\text{d}u\right)+{\sigma }_2\text{d}{B}_2\left(t\right)\\ +{\displaystyle {\int }_{\mathbb{Z}}}\left[\ln \left(\psi \left(t\right)+{\gamma }_2\left(u\right)\psi \left(t\right)\right)-\ln \psi \left(t\right)-\frac{1}{\psi \left(t\right)}{\gamma }_2\left(u\right)\psi \left(t\right)\right]\lambda \left(\text{d}u\right)\text{d}t\\ =\left[{\beta }_2\left(t\right)-\frac{f\psi \left(t\right)}{{\theta }_2}\right]\text{d}t+{\sigma }_2\text{d}{B}_2\left(t\right)+{\displaystyle {\int }_{\mathbb{Z}}}\ln \left(1+{\gamma }_2\left(u\right)\right)\tilde{N}\left(\text{d}t,\text{d}u\right).\end{array}$

That is to say

$\frac{1}{t}\ln \psi \left(t\right)=\frac{1}{t}\ln \psi \left(0\right)+\frac{1}{t}{\displaystyle {\int }_0^t}{\beta }_2\left(s\right)\text{d}s-\frac{f}{t{\theta }_2}{\displaystyle {\int }_0^t}\psi \left(s\right)\text{d}s+\frac{1}{t}{\sigma }_2{B}_2\left(t\right)+\frac{1}{t}{k}_2\left(t\right).$ (27)

For arbitrary $\epsilon >0$, there exists $T_2>0$ such that for $t\ge T_2$,

$t^{-1}\ln \psi \left(t\right)\le {\bar{\beta }}_2+2\epsilon -\frac{f}{{\theta }_2}{t}^{-1}{\displaystyle {\int }_0^t}\psi \left(s\right)\text{d}s+t^{-1}{\sigma }_2{B}_2\left(t\right)+t^{-1}{k}_2\left(t\right),$

$t^{-1}\ln \psi \left(t\right)\ge {\bar{\beta }}_2-2\epsilon -\frac{f}{{\theta }_2}{t}^{-1}{\displaystyle {\int }_0^t}\psi \left(s\right)\text{d}s+t^{-1}{\sigma }_2{B}_2\left(t\right)+t^{-1}{k}_2\left(t\right).$

Let $\epsilon$ be sufficiently small such that ${\bar{\beta }}_2-2\epsilon >0$, then applying (i) and (ii) in Lemma 3 to the above two inequalities respectively, we have

$\frac{\left({\bar{\beta }}_2-2\epsilon \right){\theta }_2}{f}\le {\langle \psi \left(t\right)\rangle }_{*}\le {\langle \psi \left(t\right)\rangle }^{*}\le \frac{\left({\bar{\beta }}_2+2\epsilon \right){\theta }_2}{f}a.s\mathrm{..}$

According the arbitrariness of $\epsilon$ gives that

$\underset{t\to +\infty }{\lim }{t}^{-1}{\displaystyle {\int }_0^t}\psi \left(s\right)\text{d}s=\underset{t\to +\infty }{\lim }\langle \psi \left(t\right)\rangle =\frac{{\bar{\beta }}_2{\theta }_2}{f}a.s\mathrm{..}$

Substituting this equation into (27), and then noting that ${\lim }_{t\to \text{}+\text{}\infty }{t}^{\text{}-\text{}1}\ln \psi \left(0\right)=0$, ${\lim }_{t\to \text{}+\text{}\infty }{B}_2\left(t\right)/t=0$ and ${\lim }_{t\to \text{}+\text{}\infty }{k}_2\left(t\right)/t=0a.s.$, we can derive that ${\lim }_{t\to \text{}+\text{}\infty }{t}^{-1}\ln \psi \left(t\right)=0a.s.$. Thus by (14), we obtain

$\underset{t\to +\infty }{\lim \inf }{t}^{-1}\ln y\left(t\right)\ge \underset{t\to +\infty }{\lim }{t}^{-1}\ln \psi \left(t\right)=0a.s\mathrm{..}$

This completes the proof.

Next we will discuss the ecological dynamics of the system (3) or system (6).

Theorem 1. Consider system (6), we have the following valid statements

(i) If ${\bar{\beta }}_1<0$ and ${\bar{\beta }}_2<0$, then both $x\left(t\right)$ and $y\left(t\right)$ are extinct, i.e.

$\underset{t\to +\infty }{\lim }x\left(t\right)=0a.s.\text{and}\underset{t\to +\infty }{\lim }y\left(t\right)=0a.s.;$

(ii) If ${\bar{\beta }}_1<0$ and ${\bar{\beta }}_2>0$, then $x\left(t\right)$ is extinct and $y\left(t\right)$ is stable in mean almost surely, i.e.

$\underset{t\to +\infty }{lim}{t}^{-1}{\displaystyle {\int }_0^t}y\left(s\right)\text{d}s=\frac{{\bar{\beta }}_2{\theta }_2}{f}a\mathrm{.}s\mathrm{.\; <mo>\; ;\; </mo>}$

(iii) If ${\bar{\beta }}_1>0$ and ${\bar{\beta }}_2<0$, then $x\left(t\right)$ is stable in mean almost surely and $y\left(t\right)$ is extinct, i.e.

$\underset{t\to +\infty }{\lim }{t}^{-1}{\displaystyle {\int }_0^t}x\left(s\right)\text{d}s=\frac{{\bar{\beta }}_1}{\eta }a.s.;$

(iv) If ${\bar{\beta }}_1>0,{\bar{\beta }}_2>0,{\left(n-c\right)}^2\le 4\left({\theta }_1-c{\theta }_2\right)$ and ${\bar{\beta }}_1>\frac{{\bar{\beta }}_2}{f}$ then $x\left(t\right)$ and $y\left(t\right)$ are both stochastic strong persistence in mean almost surely.

Proof. According to Itô’s formula, we can obtain that

$\begin{array}{l}\text{d}\ln x\left(t\right)=\left({\beta }_1\left(t\right)-\eta x\left(t\right)-\text{}\frac{cy\left(t\right)}{{\theta }_1+nx\left(t\right)+x^2\left(t\right)}\right)\text{d}t+{\sigma }_1\text{d}{B}_1\left(t\right)\\ +{\displaystyle {\int }_{\mathbb{Z}}}\ln \left(1+{\gamma }_1\left(u\right)\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\end{array}$

$\text{d}\ln y\left(t\right)=\left({\beta }_2\left(t\right)-\frac{fy\left(t\right)}{{\theta }_2+x\left(t\right)}\right)\text{d}t+{\sigma }_2\text{d}{B}_2\left(t\right)+{\displaystyle {\int }_{\mathbb{Z}}}\ln \left(1+{\gamma }_2\left(u\right)\right)\tilde{N}\left(\text{d}t,\text{d}u\right).$

In other words, we have

$\begin{array}{l}\ln x\left(t\right)-\ln x\left(0\right)={\displaystyle {\int }_0^t}{\beta }_1\left(s\right)\text{d}s-\eta {\displaystyle {\int }_0^t}x\left(s\right)\text{d}s-c{\displaystyle {\int }_0^t}\frac{y\left(s\right)}{{\theta }_1+nx\left(s\right)+x^2\left(s\right)}\text{d}s\\ +{\sigma }_1{B}_1\left(t\right)+k_1\left(t\right),\end{array}$ (28)

$\ln y\left(t\right)-\ln y\left(0\right)={\displaystyle {\int }_0^t}{\beta }_2\left(s\right)\text{d}s-f{\displaystyle {\int }_0^t}\frac{y\left(s\right)}{{\theta }_2+x\left(s\right)}\text{d}s+{\sigma }_2{B}_2\left(t\right)+k_2\left(t\right).$ (29)

Now let us prove (i). It follows from (24) and (28), for sufficiently large t, that

$t^{-1}ln\frac{x\left(t\right)}{x\left(0\right)}\le {\bar{\beta }}_1+\epsilon +t^{-1}{\sigma }_1{B}_1\left(t\right)+t^{-1}{k}_1\left(t\right)\mathrm{,}$

where $\epsilon$ is sufficiently small such that ${\bar{\beta }}_1+\epsilon <0$. Noting that ${\lim }_{t\to \text{}+\text{}\infty }{B}_1\left(t\right)/t=0$ and ${\lim }_{t\to \text{}+\text{}\infty }{k}_1\left(t\right)/t=0a\mathrm{.}s\mathrm{.}$ and, hence ${\lim }_{t\to \text{}+\text{}\infty }x\left(t\right)=0a\mathrm{.}s\mathrm{.}$. In the same way, if ${\bar{\beta }}_2<0$, and according to (29), ${\lim }_{t\to \text{}+\text{}\infty }y\left(t\right)=0a.s.$.

(ii) Since ${\bar{\beta }}_1<0$, thus (i) implies ${\lim }_{t\to \text{}+\text{}\infty }x\left(t\right)=0a.s.$. By (24) and (29), for sufficiently large t, we have

$\ln y\left(t\right)-\ln y\left(0\right)\le \left({\bar{\beta }}_2+\epsilon \right)t-\frac{f}{{\theta }_2+\epsilon }{\displaystyle {\int }_0^t}y\left(s\right)+{\sigma }_2{B}_2\left(t\right)+k_2\left(t\right),$ (30)

$\ln y\left(t\right)-\ln y\left(0\right)\ge \left({\bar{\beta }}_2-\epsilon \right)t-\frac{f}{{\theta }_2-\epsilon }{\displaystyle {\int }_0^t}y\left(s\right)+{\sigma }_2{B}_2\left(t\right)+k_2\left(t\right).$ (31)

If ${\bar{\beta }}_2>0$, then there exist arbitrarily sufficiently small $\epsilon >0$ and $T_3>0$, for all $t\ge T_3$, by making use of (i) and (ii) in Lemma 3 to (30) and (31) respectively, we have

$\frac{\left({\bar{\beta }}_2-\epsilon \right)\left({\theta }_2-\epsilon \right)}{f}\le {\langle y\left(t\right)\rangle }_{\mathrm{<mo>\; *\; </mo>}}\le {\langle y\left(t\right)\rangle }^{\mathrm{<mo>\; *\; </mo>}}\le \frac{\left({\bar{\beta }}_2+\epsilon \right)\left({\theta }_2+\epsilon \right)}{f}a\mathrm{.}s\mathrm{..}$

According to the arbitrariness of $\epsilon$, the above inequality gives

$\underset{t\to +\infty }{\lim }{t}^{-1}{\displaystyle {\int }_0^t}y\left(s\right)\text{d}s=\frac{{\bar{\beta }}_2{\theta }_2}{f}a.s\mathrm{..}$

(iii) Since ${\bar{\beta }}_2<0$, thus (i) implies ${\lim }_{t\to \text{}+\text{}\infty }y\left(t\right)=0a.s.$, By (24) and (28), for sufficiently large t, we have

$\ln x\left(t\right)-\ln x\left(0\right)\le \left({\bar{\beta }}_1+\epsilon \right)t-\eta {\displaystyle {\int }_0^t}x\left(s\right)\text{d}s+{\sigma }_1{B}_1\left(t\right)+k_1\left(t\right),$ (32)

$\ln x\left(t\right)-\ln x\left(0\right)\ge \left({\bar{\beta }}_1-\epsilon \right)t-\eta {\displaystyle {\int }_0^t}x\left(s\right)\text{d}s-\frac{c}{{\theta }_1}{\displaystyle {\int }_0^t}y\left(s\right)+{\sigma }_1{B}_1\left(t\right)+k_1\left(t\right).$ (33)

If ${\bar{\beta }}_1>0$, then there exist arbitrarily sufficiently small $\epsilon >0$ and $T_4>0$, for all $t\ge T_4$, by making use of (i) and (ii) in Lemma 3 to (32) and (33) respectively, we have

$\frac{{\bar{\beta }}_1-\epsilon }{\eta }\le {\langle x\left(t\right)\rangle }_{\mathrm{<mo>\; *\; </mo>}}\le {\langle x\left(t\right)\rangle }^{\mathrm{<mo>\; *\; </mo>}}\le \frac{{\bar{\beta }}_1+\epsilon }{\eta }a\mathrm{.}s\mathrm{..}$

According to the arbitrariness of $\epsilon$, the above inequality gives

$\underset{t\to +\infty }{lim}{t}^{-1}{\displaystyle {\int }_0^t}x\left(s\right)\text{d}s=\frac{{\bar{\beta }}_1}{\eta }a\mathrm{.}s\mathrm{..}$

(iv) From (29), we get

$\begin{array}{l}{t}^{-1}\left(\ln y\left(t\right)-\ln y\left(0\right)\right)\\ =t^{-1}{\displaystyle {\int }_0^t}{\beta }_2\left(s\right)\text{d}s-t^{-1}f{\displaystyle {\int }_0^t}\frac{y\left(s\right)}{{\theta }_2+x\left(s\right)}\text{d}s+t^{-1}{\sigma }_2{B}_2\left(t\right)+t^{-1}{k}_2\left(t\right).\end{array}$

Then, we have

$\begin{array}{l}{t}^{\text{}-\text{}1}f{\displaystyle {\int }_0^t}\frac{y\left(s\right)}{{\theta }_2+x\left(s\right)}\text{d}s\\ =-t^{\text{}-\text{}1}\left(\ln y\left(t\right)-\ln y\left(0\right)\right)+t^{\text{}-\text{}1}{\displaystyle {\int }_0^t}{\beta }_2\left(s\right)\text{d}s+t^{\text{}-\text{}1}{\sigma }_2{B}_2\left(t\right)+t^{\text{}-\text{}1}{k}_2\left(t\right),\end{array}$

noting that ${\lim }_{t\to \text{}+\text{}\infty }{t}^{-1}\ln y\left(0\right)=0$, ${\lim }_{t\to \text{}+\text{}\infty }{B}_2\left(t\right)/t=0$ and ${\lim }_{t\to \text{}+\text{}\infty }{k}_2\left(t\right)/t=0a\mathrm{.}s\mathrm{.}$ and (23), for sufficiently large t, we can derive that

$\underset{t\to +\infty }{lim}{t}^{-1}f{\displaystyle {\int }_0^t}\frac{y\left(s\right)}{{\theta }_2+x\left(s\right)}\text{d}s=\frac{{\bar{\beta }}_2}{f}a\mathrm{.}s\mathrm{..}$ (34)

Moreover, follows from (24) and (29), for sufficiently large t, we have

$\ln y\left(t\right)-\ln y\left(0\right)\ge \left({\bar{\beta }}_2-\epsilon \right)t-\frac{f}{{\theta }_2}{\displaystyle {\int }_0^t}y\left(s\right)\text{d}s+{\sigma }_2{B}_2\left(t\right)+k_2\left(t\right),$

then there exist arbitrarily sufficiently small $\epsilon >0$ and $T_5>0$, for all $t\ge T_5$, by making use of (ii) in Lemma 3, we get

${\langle y\left(t\right)\rangle }_{\mathrm{<mo>\; *\; </mo>}}\ge \frac{{\bar{\beta }}_2{\theta }_2}{f}a\mathrm{.}s\mathrm{..}$

According to the given condition ${\left(n-c\right)}^2\le 4\left({\theta }_1-c{\theta }_2\right)$, we get

$\begin{array}{c}\text{d}\ln x\left(t\right)=\left({\beta }_1\left(t\right)-\eta x\left(t\right)-\frac{cy\left(t\right)}{{\theta }_1+nx\left(t\right)+x^2\left(t\right)}\right)\text{d}t+{\sigma }_1\text{d}{B}_1\left(t\right)\\ +{\displaystyle {\int }_{\mathbb{Z}}}\ln \left(1+{\gamma }_1\left(u\right)\right)\tilde{N}\left(\text{d}t,\text{d}u\right),\end{array}$ (35)

$\begin{array}{c}\ge \left({\beta }_1\left(t\right)-\eta x\left(t\right)-\frac{y\left(t\right)}{{\theta }_2+x\left(t\right)}\right)\text{d}t+{\sigma }_1\text{d}{B}_1\left(t\right)\\ +{\displaystyle {\int }_{\mathbb{Z}}}\ln \left(1+{\gamma }_1\left(u\right)\right)\tilde{N}\left(\text{d}t,\text{d}u\right).\end{array}$

Integrating the both side of (35), then by (34) and for sufficiently large t, we have

$lnx\left(t\right)-lnx\left(0\right)\ge \left({\bar{\beta }}_1-\frac{{\bar{\beta }}_2}{f}-2\epsilon \right)t-\eta {\displaystyle {\int }_0^t}x\left(s\right)\text{d}s+{\sigma }_1{B}_1\left(t\right)+k_1\left(t\right)\mathrm{.}$ (36)

There exist arbitrarily sufficiently small $\epsilon >0$ and $T_5>0$, for all $t\ge T_5$, by making use of (i) and (ii) in Lemma 3 to (33) and (36) respectively, we get

$\frac{{\bar{\beta }}_1-\frac{{\bar{\beta }}_2}{f}}{\eta }\le {\langle x\left(t\right)\rangle }_{\mathrm{<mo>\; *\; </mo>}}\le {\langle x\left(t\right)\rangle }^{\mathrm{<mo>\; *\; </mo>}}\le \frac{{\bar{\beta }}_1}{\eta }a\mathrm{.}s\mathrm{..}$

5. Numerical Simulations

In this section, we apply Split-step Backward Euler method [28] [29] [30] to prove our theoretical results.

(1) We assume the parameters $a_1=0.8$, $a_2=0.6$, $a_{11}=a_{21}=0.6$, $\eta =0.3$, $c=0.2$, $f=1.1$, $k_1=0.4$, $k_2=0.5$, $l_1=0.5$, $l_2=0.4$, $m_1=0.5$, $m_2=0.6$, $h=0.8$, $q=0.2$, $\tau =1$, ${\theta }_1={\theta }_2=1$, $n=1$, ${\sigma }_1={\sigma }_2=2$, ${\gamma }_1\left(u\right)\equiv 2$, ${\gamma }_2\left(u\right)\equiv 2$, then ${\bar{\beta }}_1=-0.2414<0$, ${\bar{\beta }}_2=-0.4564<0$.

We observe that two species will go to extinction from Figure 1, and the result of (i) in Theorem 1 are shown.

(2) Let ${\gamma }_1\left(u\right)\equiv 2$, ${\gamma }_2\left(u\right)\equiv 2$, then other conditions remain unchanged, we have ${\bar{\beta }}_1=-0.2414<0$, ${\bar{\beta }}_2=0.4273>0$.

We observe that $x\left(t\right)$ will go to extinction, and $y\left(t\right)$ will be stable in mean from Figure 2, where

$\underset{t\to +\infty }{lim}{t}^{-1}{\displaystyle {\int }_0^t}y\left(s\right)\text{d}s=\frac{{\bar{\beta }}_2{\theta }_2}{f}=\mathrm{0.3885,}$

and then the result of (ii) in Theorem 1 are shown.

(3) We set ${\gamma }_1\left(u\right)\equiv 0.1$, ${\gamma }_2\left(u\right)\equiv 1.5$, then other conditions remain unchanged, we have ${\bar{\beta }}_1=0.6553>0$, ${\bar{\beta }}_2=0.1387<0$.

From Figure 3, we get that $x\left(t\right)$ will be stable in mean and $y\left(t\right)$ will be extinct, and then the result of (iii) in Theorem 1 are shown, where

$\underset{t\to +\infty }{\lim }{t}^{-1}{\displaystyle {\int }_0^t}x\left(s\right)\text{d}s=\frac{{\bar{\beta }}_1}{\eta }=\mathrm{2.1844.}$

(4) We set ${\gamma }_1\left(u\right)\equiv 0.03$, ${\gamma }_2\left(u\right)\equiv 0.01$, then other conditions remain unchanged, we have ${\bar{\beta }}_1=0.6596>0$, ${\bar{\beta }}_2=0.4450>0$, then ${\left(n-c\right)}^2-4\left({\theta }_1-c{\theta }_2\right)=-2.5600<0$, $0.6596={\bar{\beta }}_1>\frac{{\bar{\beta }}_2}{f}=0.4045$.

From Figure 4, we get that both $x\left(t\right)$ and $y\left(t\right)$ will be stochastic strong

Figure 1. This figure is time series graph of $x\left(t\right)$ and $y\left(t\right)$. We choose ${\gamma }_1\left(u\right)\equiv 2$, ${\gamma }_2\left(u\right)\equiv 2$, step size $\Delta t=0.03$, initial value $x\left(0\right)=0.6$, $y\left(0\right)=0.6$, $C_1\left(0\right)=0.2$, $C_2\left(0\right)=0.3$, and $C_e\left(0\right)=0.65$, $\mathbb{Z}=\left(0,+\infty \right)$, $\lambda \left(\mathbb{Z}\right)=1$.

Figure 2. This figure is time series of $x\left(t\right)$ and $y\left(t\right)$ with ${\gamma }_1\left(u\right)\equiv 2$, ${\gamma }_2\left(u\right)\equiv 0.2$, step size $\Delta t=0.03$, initial value $x\left(0\right)=0.6$, $y\left(0\right)=0.6$, $C_1\left(0\right)=0.2$, $C_2\left(0\right)=0.3$, and $C_e\left(0\right)=0.65$, $\mathbb{Z}=\left(0,+\infty \right)$, $\lambda \left(\mathbb{Z}\right)=1$.

Figure 3. This figure is time series of $x\left(t\right)$ and $y\left(t\right)$ with ${\gamma }_1\left(u\right)\equiv 0.1$, ${\gamma }_2\left(u\right)\equiv 1.5$, step size $\Delta t=0.03$, initial value $x\left(0\right)=0.6$, $y\left(0\right)=0.6$, $C_1\left(0\right)=0.2$, $C_2\left(0\right)=0.3$, and $C_e\left(0\right)=0.65$, $\mathbb{Z}=\left(0,+\infty \right)$, $\lambda \left(\mathbb{Z}\right)=1$.

Figure 4. The above figure is time series graph of $x\left(t\right)$ and $y\left(t\right)$. The under figure is the 2D phase diagram of $x\left(t\right)$ and $y\left(t\right)$. We choose ${\gamma }_1\left(u\right)\equiv 0.03$, ${\gamma }_2\left(u\right)\equiv 0.01$, step size $\Delta t=0.03$, initial value $x\left(0\right)=0.6$, $y\left(0\right)=0.6$, $C_1\left(0\right)=0.2$, $C_2\left(0\right)=0.3$, and $C_e\left(0\right)=0.65$, $\mathbb{Z}=\left(0,+\infty \right)$, $\lambda \left(\mathbb{Z}\right)=1$.

persistent in mean. The 2D phase diagram of $x\left(t\right)$ and $y\left(t\right)$ mean that two species are in a predator-prey relationship. The result of (iv) in Theorem 1 are shown, where

${\langle y\left(t\right)\rangle }_{*}\ge \frac{{\bar{\beta }}_2{\theta }_2}{f}=0.4045,$

$0.8502=\frac{{\bar{\beta }}_1-\frac{{\bar{\beta }}_2}{f}}{\eta }\le {\langle x\left(t\right)\rangle }_{*}\le {\langle x\left(t\right)\rangle }^{*}\le \frac{{\bar{\beta }}_1}{\eta }=\mathrm{2.1985.}$

6. Discussion and Conclusions

In this paper, we add Lévy noise to the stochastic modified Leslie-Gower and Holling-type IV predator-prey model, and assume that the toxicants are added in periodic pulses in the model. We show that the model has a unique global solution and study the stability in distribution of solutions. We get the thresholds ${\bar{\beta }}_i$ to determine extinction and persistent in mean of two species; thus sufficient and necessary conditions are established for the extinction and persistent in mean of two species.

From the Theorem 1 and the numerical simulation results in Figures 1-4, we can see that Lévy noise has a strong effect on the system (3). At the same time, through the expression of the thresholds ${\bar{\beta }}_i$ and changing the parameter value multiple times, it shows that the line shape in the numerical simulation is undulating, because white noise can reflect that the model is affected by the environment. We also know that the value of $n\mathrm{,}c\mathrm{,}{\theta }_1\mathrm{,}{\theta }_2$ and f will affect the survival dynamics of the species from (iv) in Theorem 1. The expression of ${\bar{\beta }}_i$ also reflects that the toxicants and population’s own performance also more or less affects the survival dynamics of the species.

Indeed, when the population encounters sudden environmental disturbances, such as tsunamis, earthquakes, etc., the survival environment of the population is threatened. Ecological stability is bound to be affected because they can’t adapt to this sudden environmental fluctuation in a short time. Lévy jump has a great impact on the survival of species. With the rapid development of modern industrial technology, pollution has been increased as well. Impulsive toxicant will inevitably have a certain impact to species’ living environment and their own growth.

This article has practical significance for the survival analysis of a stochastic modified Leslie-Gower and Holling-type IV predator-prey model with Lévy noise in impulsive toxicant input environments. But considering that some more complex systems will be more in line with the actual situation, for example, during the rainy season and the dry season, the growth rate and mortality rate of the species are different, so we can consider adding the regime switching to the system (3). In the next research work, we can try to consider the influence of continuous-time Markov chain on the system.

Acknowledgements

The authors would like to thank the referees for their useful suggestions which have significantly improved the paper.

Conflicts of Interest

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

References

[1]	Sarker, A., et al. (2020) A Review of Microplastics Pollution in the Soil and Terrestrial Ecosystems: A Global and Bangladesh Perspective. Science of the Total Environment, 733, Article ID: 139296. https://doi.org/10.1016/j.scitotenv.2020.139296
[2]	Kumar, V., Parihar, R.D., Sharma, A., et al. (2019) Global Evaluation of Heavy Metal Content in Surface Water Bodies: A Meta-Analysis Using Heavy Metal Pollution Indices and Multivariate Statistical Analyses. Chemosphere, 236, Article ID: 124364. https://doi.org/10.1016/j.chemosphere.2019.124364
[3]	Jiang, R., Wang, M.E., Chen, W.P., et al. (2019) Ecological Risk of Combined Pollution on Soil Ecosystem Functions: Insight from the Functional Sensitivity and Stability. Environmental Pollution, 255, Article ID: 113184. https://doi.org/10.1016/j.envpol.2019.113184
[4]	Zhai, T.L., Wang, J., Fang, Y., Qin, Y., Huang, L.Y. and Chen, Y. (2020) Assessing Ecological Risks Caused by Human Activities in Rapid Urbanization Coastal Areas: Towards an Integrated Approach to Determining Key Areas of Terrestrial-Oceanic Ecosystems Preservation and Restoration. Science of the Total Environment, 708, Article ID: 135153. https://doi.org/10.1016/j.scitotenv.2019.135153
[5]	Wan, N.F., Ji, X.Y., Jiang, J.X., Zhang, Y.M., Liang, J.H. and Li, B. (2015) An Ecological Indicator to Evaluate the Effect of Chemical Insecticide Pollution Management on Complex Ecosystems. Ecological Indicators, 53, 11-17. https://doi.org/10.1016/j.ecolind.2015.01.014
[6]	Ma, H., Pu, S.Y., Liu, S.B., Bai, Y.C., Mandal, S.and Xing, B.S. (2020) Microplastics in Aquatic Environments: Toxicity to Trigger Ecological Consequences. Environmental Pollution, 261, Article ID: 114089. https://doi.org/10.1016/j.envpol.2020.114089
[7]	Jiang, K., Merrill, R., You, D.M., Pan, P. and Li, Z.D. (2019) Optimal Control for Transboundary Pollution under Ecological Compensation: A Stochastic Differential Game Approach. Journal of Cleaner Production, 241, Article ID: 118391. https://doi.org/10.1016/j.jclepro.2019.118391
[8]	He, J. and Wang, K. (2009) The Survival Analysis for a Population in a Polluted Environment. Nonlinear Analysis: Real World Applications, 10, 1555-1571. https://doi.org/10.1016/j.nonrwa.2008.01.027
[9]	Liu, M. and Wang, K. (2009) Survival Analysis of Stochastic Single-Species Population Models in Polluted Environments. Ecological Modelling, 220, 1347-1357. https://doi.org/10.1016/j.ecolmodel.2009.03.001
[10]	Liu, M. (2015) Survival Analysis of a Cooperation System with Random Perturbations in a Polluted Environment. Nonlinear Analysis: Hybrid Systems, 18, 100-116. https://doi.org/10.1016/j.nahs.2015.06.005
[11]	Zhang, S. and Tan, D. (2015) Dynamics of a Stochastic Predator-Prey System in a Polluted Environment with Pulse Toxicant Input and Impulsive Perturbations. Applied Mathematical Modelling, 39, 6319-6331. https://doi.org/10.1016/j.apm.2014.12.020
[12]	Lv, X.J., Meng, X.Z. and Wang, X.Z. (2018) Extinction and Stationary Distribution of an Impulsive Stochastic Chemostat Model with Nonlinear Perturbation. Chaos, Solitons & Fractals, 110, 273-279. https://doi.org/10.1016/j.chaos.2018.03.038
[13]	Liu, M., Du, C.X. and Deng, M.L. (2018) Persistence and Extinction of a Modified Leslie-Gower Holling-Type II Stochastic Predator? Cprey Model with Impulsive Toxicant Input in Polluted Environments. Nonlinear Analysis: Hybrid Systems, 27, 177-190. https://doi.org/10.1016/j.nahs.2017.08.001
[14]	Bao, J.H. and Yuan, C.G. (2012) Stochastic Population Dynamics Driven by Lévy Noise. Journal of Mathematical Analysis and Applications, 391, 363-375. https://doi.org/10.1016/j.jmaa.2012.02.043
[15]	Bao, J.H., Mao, X.R., Yin, G. and Yuan, C.G. (2011) Competitive Lotka-Volterra Population Dynamics with Jumps. Nonlinear Analysis: Theory, Methods & Applications, 74, 6601-6616. https://doi.org/10.1016/j.na.2011.06.043
[16]	Zhao, Y., Yuan, S.L. and Zhang, Q.M. (2016) The Effect of Lévy Noise on the Survival of a Stochastic Competitive Model in an Impulsive Polluted Environment. Applied Mathematical Modelling, 40, 7583-7600. https://doi.org/10.1016/j.apm.2016.01.056
[17]	Xu, D.S., Liu, M. and Xu, X.F. (2020) Analysis of a Stochastic Predator-Prey System with Modified Leslie-Gower and Holling-Type IV Schemes. Physica A: Statistical Mechanics and its Applications, 537, Article ID: 122761.
[18]	Yang, X.F., Jin, Z. and Xue, Y.K. (2007) Weak Average Persistence and Extinction of a Predator-Prey System in a Polluted Environment with Impulsive Toxicant Input. Chaos, Solitons & Fractals, 31, 726-735. https://doi.org/10.1016/j.chaos.2005.10.042
[19]	Liu, B., Chen, L.S. and Zhang, Y.J. (2003) The Effects of Impulsive Toxicant Input on a Population in a Polluted Environment. Journal of Biological Systems, 11, 265-274. https://doi.org/10.1142/S0218339003000907
[20]	Liptser, R.S. (1980) A Strong Law of Large Numbers for Local Martingales. Stochastics, 3, 217-228. https://doi.org/10.1080/17442508008833146
[21]	Liu, M. and Wang, K. (2014) Stochastic Lotka-Volterra Systems with Lévy Noise. Journal of Mathematical Analysis and Applications, 410, 750-763. https://doi.org/10.1016/j.jmaa.2013.07.078
[22]	Liu, M. and Wang, K. (2012) Global Asymptotic Stability of a Stochastic Lotka-Volterra Model with Infinite Delays. Communications in Nonlinear Science and Numerical Simulation, 17, 3115-3123. https://doi.org/10.1016/j.cnsns.2011.09.021
[23]	Chen, L. and Chen, J. (1993) Nonlinear Biological Dynamical System. Science Press, Beijing.
[24]	Ikeda, N. and Watanabe, S. (2014) Stochastic Differential Equations and Diffusion Processes. Elsevier, Amsterdam.
[25]	Li, X.Y. and Mao, X.R. (2009) Population Dynamical Behavior of Non-Autonomous Lotka-Volterra Competitive System with Random Perturbation. American Institute of Mathematical Sciences, 24, 523-545. https://doi.org/10.3934/dcds.2009.24.523
[26]	Barbalat, I. (1959) Systems d’equations differentielles d’osci d’oscillations nonlinéaires. Revue Roumaine de Mathématique Pures et Appliquées, 4, 267-270.
[27]	Kot, M. (2001) Elements of Mathematical Biology, Mathematics in Science and Engineering. Cambridge University Press, Cambridge.
[28]	Bastani, A.F. and Tahmasebi, M. (2012) Strong Convergence of Split-Step Backward Euler Method for Stochastic Differential Equations with Non-Smooth Drift. Journal of Computational and Applied Mathematics, 236, 1903-1918. https://doi.org/10.1016/j.cam.2011.10.023
[29]	Zhang, Q.M., Liu, Y.T. and Li, X.N. (2014) Strong Convergence of Split-Step Backward Euler Method for Stochastic Age-Dependent Capital System with Markovian Switching. Applied Mathematics and Computation, 235, 439-453. https://doi.org/10.1016/j.amc.2013.12.189
[30]	Tan, J.G. and Wang, H.L. (2010) Convergence and Stability of the Split-Step Backward Euler Method for Linear Stochastic Delay Integro-Differential Equations. Mathematical and Computer Modelling, 51, 504-515. https://doi.org/10.1016/j.mcm.2009.11.020