1. Introduction

As far as we know, rate-dependent predator-prey system models have become the focus of mathematical ecology and have been extensively studied in recent years (see e.g., [1] [2] [3] ). The dynamic relationship between predator and prey is ubiquitous in ecology and mathematical ecology [4]. The relationship between two species is usually thought of as competition, predation and cooperation. Here is the Lotka-Volterra model of a predator-predator with intra-species competition:

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=a_{1}x-b_1{x}^2-c_{1}xy,\hfill \\ \frac{\text{d}y}{\text{d}t}=-a_{2}y-b_2{y}^2+c_{2}xy,\hfill \end{array}$ (1.1)

where $x\left(t\right)$ and $y\left(t\right)$ represent the population density of prey species and predators species at time t, parameters $a_1$, $a_2$, $c_1$ and $c_2$ are all positive constants describing the interaction of two species. The positive constants $b_1$, $b_2$ represent the effect of one species on the other. It is well known that the solution of system (1.1) is asymptotically stable.

In the study of biological phenomena, there are many factors affecting the dynamic properties of biological and mathematical models, and functional reaction is one of the common nonlinear factors [4]. System (1.1) assumes the prey biomass is enough and a individual predator consume the prey with functional response of type $c_{1}x$. When predators face an increase in the density of their local prey, they usually change their consumption rate. The concept of functional response was first proposed by Solomon and later discussed in detail by Holling. Holling [5] proposed three types of functional responses, i.e., Hollings type

I of $\alpha x$, Hollings type II of $\frac{\alpha x}{\beta +x}$ and Hollings type III of $\frac{\alpha x}{\beta +x^2}$. Proportional-dependent functional responses are a better description of how predators

must find food and therefore must share or compete for it. Based on Holling type II function, Arditi and Ginzburg [2] first proposed a rate-dependent functional response model:

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=a_{1}x-b_1{x}^2-\frac{c_{1}xy}{x+ey},\hfill \\ \frac{\text{d}y}{\text{d}t}=-a_{2}y+\frac{c_{2}xy}{x+ey},\hfill \end{array}$ (1.2)

where parameters $c_1,c_2,e$ are all positive constants, representing capturing rate, conversion rate and half capturing saturation constant, respectively.

These papers are all deterministic models, which do not consider the impact of environmental fluctuations, nor the impact of population randomness. Population dynamics in nature will inevitably be affected by environmental noise in the ecosystem. More recently, a number of authors have looked at stochastic predator-prey models with white noise and revealed how white noise affects population systems, such as [4] [6] [7] [8]. Considering that the environmental fluctuation is mainly manifested as the fluctuation of the internal growth rate of the predator population and the mortality rate of the predator population [9], they supposed parameters $a_1$ and $a_2$ were perturbed with

$a_1\to a_1+{\sigma }_{1}B\left(t\right)\mathrm{,}{a}_2\to a_2+{\sigma }_{2}B\left(t\right)\mathrm{,}$ (1.3)

where $B\left(t\right)$ is a standard Brownian process defined on a complete probability space $\left(\Omega \mathrm{,}\mathcal{F}\mathrm{,}\mathbb{P}\right)$ with a filtration ${\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}$ satisfying the usual conditions, ${\sigma }_1^2$ and ${\sigma }_2^2$ represent the intensities of the white noise.

Nguyen and Ta [7] introduced intra-specific competition into the stochastic rate-dependent model to obtain the model (4), and considered the corresponding non-autonomous stochastic system, and estimated the high population growth rate and exponential mortality rate. The stochastic predator-prey system with rate-dependent functional response can be expressed as follows:

$\{\begin{array}{l}\text{d}x=x\left[a_1-b_{1}x-\frac{c_{1}y}{x+ey}\right]\text{d}t+{\sigma }_{1}x\text{d}B\left(t\right),\hfill \\ \text{d}y=y\left[-a_2-b_{2}y+\frac{c_{2}x}{x+ey}\right]\text{d}t+{\sigma }_{2}y\text{d}B\left(t\right).\hfill \end{array}$ (1.4)

Furthermore, from a biological point of view, population dynamics may encounter sudden environmental disturbances that cannot be described by white noise, such as earthquakes, epidemics, floods and hurricanes. In this case, there are many references to stochastic differential equations with jumps (see e.g., [9] [10] [11] [12] [13] ), and Lévy jumping into a potential population system may be a reasonable method to describe these phenomena. Therefore, this paper considers the random ratio-dependent model with jumps:

$\{\begin{array}{l}\text{d}x\left(t\right)=x\left(t^{-}\right)\left\{\left[a_1-b_{1}x\left(t^{-}\right)-\frac{c_{1}y\left(t^{-}\right)}{x\left(t^{-}\right)+ey\left(t^{-}\right)}\right]\text{d}t+{\sigma }_1\text{d}B\left(t\right)+{\displaystyle {\int }_{\mathbb{Y}}}{\gamma }_1\left(u\right)\tilde{N}\left(\text{d}t\mathrm{,}\text{d}u\right)\right\}\mathrm{,}\\ \text{d}y\left(t\right)=y\left(t^{-}\right)\left\{\left[-a_2-b_{2}y\left(t^{-}\right)+\frac{c_{2}x\left(t^{-}\right)}{x\left(t^{-}\right)+ey\left(t^{-}\right)}\right]\text{d}t+{\sigma }_2\text{d}B\left(t\right)+{\displaystyle {\int }_{\mathbb{Y}}}{\gamma }_2\left(u\right)\tilde{N}\left(\text{d}t\mathrm{,}\text{d}u\right)\right\}\mathrm{,}\end{array}$ (1.5)

where $x\left(t^{-}\right)$ and $y\left(t^{-}\right)$ represent the left limit of $x\left(t\right)$ and $y\left(t\right)$, respectively: the parameters $a_l\left(i\right)$, $b_l\left(i\right)$, $c_l\left(i\right)$, $e\left(i\right)$ and ${\sigma }_l\left(i\right)\left(l=1,2\right)$ are all positive constant, N is a poisson counting measure with compensator $\tilde{N}$ and characteristic measure $\lambda$ on a measurable subset $\mathbb{Y}$ of $\left(\mathrm{0,}\infty \right)$ with $\lambda \left(\mathbb{Y}\right)<\infty$, and $\tilde{N}\left(\text{d}t\mathrm{,}\text{d}u\right)=N\left(\text{d}t\mathrm{,}\text{d}ut\right)-\lambda \left(\text{d}u\right)\text{d}t$, the function ${\gamma }_l\mathrm{<mo>\; :\; </mo>}\mathbb{Y}\times \left(\mathrm{0,}\infty \right)\to ℝ$ is bounded and continuous with respect to $\lambda$ and is $\mathfrak{B}\left(\mathbb{Y}\right)\times {\mathcal{F}}_t$ -measurable, $l=1,2$.

Now let’s go one step further and add another environmental noise that is Telegraph noise. Telegraph noise can be described as a random switch between two or more environmental states that differ in factors such as nutrition or rainfall [14] [15]. The stochastic differential equations driven by a continuous-time Markov chain have been used to model the population system [16] [17] [18] [19] with this type of noise. Suppose the Markov chain on the state space $\mathbb{S}=\left\{\mathrm{1,2,}\cdots \mathrm{,}S\right\}$ controls the switching between the environmental regimes. Then the prey-predator model with two types of noise can be described by the following stochastic differential equation

$\{\begin{array}{l}\text{d}x\left(t\right)=x\left(t^{-}\right)\left\{\left[a_1\left(r\left(t\right)\right)-b_1\left(r\left(t\right)\right)x\left(t^{-}\right)-\frac{c_1\left(r\left(t\right)\right)y\left(t^{-}\right)}{x\left(t^{-}\right)+e\left(r\left(t\right)\right)y\left(t^{-}\right)}\right]\text{d}t\right\}\\ +{\sigma }_1\left(r\left(t\right)\right)x\left(t^{-}\right)\text{d}B\left(t\right)+{\displaystyle {\int }_{\mathbb{Y}}}{\gamma }_1\left(u\mathrm{,}r\left(t\right)\right)x\left(t^{-}\right)\tilde{N}\left(\text{d}t\mathrm{,}\text{d}u\right)\mathrm{,}\\ \text{d}y\left(t\right)=y\left(t^{-}\right)\left\{\left[-a_2\left(r\left(t\right)\right)-b_2\left(r\left(t\right)\right)y\left(t^{-}\right)+\frac{c_2\left(r\left(t\right)\right)x\left(t^{-}\right)}{x\left(t^{-}\right)+e\left(r\left(t\right)\right)y\left(t^{-}\right)}\right]\text{d}t\right\}\\ +{\sigma }_2\left(r\left(t\right)\right)y\left(t^{-}\right)\text{d}B\left(t\right)+{\displaystyle {\int }_{\mathbb{Y}}}{\gamma }_1\left(u\mathrm{,}r\left(t\right)\right)y\left(t^{-}\right)\tilde{N}\left(\text{d}t\mathrm{,}\text{d}u\right)\mathrm{,}\end{array}$ (1.6)

with initial data

$X\left(0\right)=\left(x\left(0\right)\mathrm{,}y\left(0\right)\right)\in {ℝ}_{+}^2\mathrm{,}r\left(0\right)\in \mathbb{S}\mathrm{,}$ (1.7)

where $r\left(t\right)$ is a right-continuous Markov chain, taking values in $\mathbb{S}=\left\{\mathrm{1,2,}\cdots \mathrm{,}S\right\}$. System (1.6) will switch from one mode to another according to the law of Markov chain. If $r\left(0\right)=i_0$, then system (1.6) obeys

$\{\begin{array}{l}\text{d}x\left(t\right)=x\left(t^{-}\right)\{\left[a_1\left(i\right)-b_1\left(i\right)x\left(t^{-}\right)-\frac{c_1\left(i\right)y\left(t^{-}\right)}{x\left(t^{-}\right)+e\left(i\right)y\left(t^{-}\right)}\right]\text{d}t\\ \begin{array}{c}\\ \underset{}{}\end{array}+{\sigma }_1\left(i\right)\text{d}B\left(t\right)+{\displaystyle {\int }_{\mathbb{Y}}}{\gamma }_1\left(u\mathrm{,}i\right)\tilde{N}\left(\text{d}t\mathrm{,}\text{d}u\right)\}\\ \text{d}y\left(t\right)=y\left(t^{-}\right)\{\left[-a_2\left(i\right)-b_2\left(i\right)y\left(t^{-}\right)+\frac{c_2\left(i\right)x\left(t^{-}\right)}{x\left(t^{-}\right)+e\left(i\right)y\left(t^{-}\right)}\right]\text{d}t\\ \begin{array}{c}\\ \underset{}{}\end{array}+{\sigma }_2\left(i\right)\text{d}B\left(t\right)+{\displaystyle {\int }_{\mathbb{Y}}}{\gamma }_2\left(u\mathrm{,}i\right)\tilde{N}\left(\text{d}t\mathrm{,}\text{d}u\right)\}\end{array}$ (1.8)

with $i=i_0$ until time ${\tau }_1$ when Markov chain jumps to $i_1$ from $i_0$ ; the system will then obey (8) with $i=i_1$ from ${\tau }_1$ to ${\tau }_2$ when the Markov chain jumps to $i_2$ from $i_1$. The system will continue to switch as long as the Markov chain jumps. In other words, system (8) can be regarded as system (6) switching from one to another according to the law of the Markov chain. If the switching between environmental regimes disappears, then system (6) degenerates into system (8).

More recently, several authors have discussed the effects of environmental fluctuations on populations [16] [20] [21] [22] [23]. C. Ji et al. [21] studied a stochastic predator-prey system with white noise and concluded that the predator population and ratio-dependent functions were stable on average time. M. Ouyang, X. Li [22] studied stochastic predator-prey systems with Markovian switching, and obtained that the stable distribution of Markov chains was related to the parameters of subsystems, and the switching between subsystems made them neither permanent nor dissipative. L. Bai et al. [23] studied a stochastic predator-prey system with Lévy noise and revealed that Lévy noise and white noise can have an impact on biological systems. But previous results were insufficient to reveal the effects of Markovian switching and Lévy noise on prey and predators. Thus, these factors let do investigate the dynamics of the stochastic random predator-prey system described in system (1.8). Till now, the available tools for analyzing stochastic population models with Markovian switching are limited [24] [25]. The key methods used in this paper are m-matrix analysis (e.g., [18] ) and stochastic analysis of Lyapunov functions developed by Khasminskii (see e.g., [18] [26] ).

This paper is organized as follows: In Section 2, we give the global existence and positive properties of the solutions of system (1.8). In section 3 and section 4, we give sufficient conditions for the non-persistence, weak persistence and extinction and stochastic permanence respectively. In Section 5, we demonstrate that system (5.1) has a unique stationary distribution under some appropriate conditions when Lévy jumps are not present. Finally, we illustrate our main results with two examples.

2. Preliminaries

Throughout this paper, let $\left(\Omega \mathrm{,}\mathcal{F}\mathrm{,}{\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}\mathrm{,}P\right)$ be a complete probability space with a filtration ${\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}$ satisfying the usual conditions (i.e. it is right continuous and ${\mathcal{F}}_0$ contains all P-null sets). Let $B\left(t\right)$, $t\ge 0$, be a scalar standard Brownian motion defined on this probability space. We also denote by ${ℝ}_{+}^2$ the positive cone in ${ℝ}^2$, and denote by ${\overline{ℝ}}_{+}^2$ the nonnegative cone in ${ℝ}^2$.

Let $r\left(t\right)$ be a right-continuous Markov chain on the probability space, taking values in a finite state space $\mathbb{S}=\left\{\mathrm{1,2,}\cdots \mathrm{,}S\right\}$ with the generator $\Gamma ={\left({\gamma }_{ij}\right)}_{S\times S}$ given by

$\mathbb{P}\left(r\left(t+\zeta \right)=j|r\left(t\right)=i\right)=\{\begin{array}{l}{\gamma }_{ij}+o\left(\zeta \right),i\ne j;\hfill \\ 1+{\gamma }_{ij}\zeta +o\left(\zeta \right),i=j,\hfill \end{array}$ (2.1)

where $\zeta >0$. Here ${\gamma }_{ij}$ represents the transition rate from i to j and ${\gamma }_{ij}>0$, if $i\ne j$ while

${\gamma }_{ii}=-{\displaystyle {\sum }_{j=1,j\ne i}^S{\gamma }_{ij}},\forall i\in \mathbb{S}.$ (2.2)

We assume that the Markov chain $r(\cdot )$ is independent of the Brownian motion $B(\cdot )$ and is irreducible. Under this condition, the Markov chain has a unique stationary (probability) distribution $\phi =\left({\phi }_1\mathrm{,}{\phi }_2\mathrm{,}\cdots \mathrm{,}{\phi }_S\right)\in R^{1\times S}$ which can be determined by solving the following linear equation

$\phi \Gamma =0,$ (2.3)

subject to

${\displaystyle {\sum }_{i=1}^S{\phi }_i}=1\text{and}{\phi }_i>\mathrm{0,}\forall i\in \mathbb{S}\mathrm{.}$ (2.4)

Next, let $x\left(t\right)\in {ℝ}^2$ be a solution of the stochastic differential equation with regime-switching jump-diffusion processes taking the form

$\begin{array}{c}\text{d}x\left(t\right)=F\left(x\left(t^{-}\right)\mathrm{,}{t}^{-}\mathrm{,}r\left(t^{-}\right)\right)\text{d}t+G\left(x\left(t^{-}\right)\mathrm{,}{t}^{-}\mathrm{,}r\left(t^{-}\right)\right)\text{d}B\left(t\right)\\ +{\displaystyle {\int }_{\mathbb{Y}}}H\left(x\left(t^{-}\right)\mathrm{,}{t}^{-}\mathrm{,}r\left(t^{-}\right)\mathrm{,}u\right)\tilde{N}\left(\text{d}t\mathrm{,}\text{d}u\right)\mathrm{,}\end{array}$ (2.5)

where $F\mathrm{<mo>\; :\; </mo>}{ℝ}^2\times {ℝ}_{+}\times \mathbb{S}\to {ℝ}^2$, $G\mathrm{<mo>\; :\; </mo>}{ℝ}^2\times {ℝ}_{+}\times \mathbb{S}\to {ℝ}^2$ and $H\mathrm{<mo>\; :\; </mo>}{ℝ}^2\times {ℝ}_{+}\times \mathbb{S}\times \mathbb{Y}\to {ℝ}^2$ are measurable functions. Moreover, let $C^2\left({ℝ}^2\times {ℝ}_{+}\times \mathbb{S}\mathrm{<mo>\; ;\; </mo>}ℝ\right)$ denote the family of all real-valued function $V\left(x\mathrm{,}t\mathrm{,}i\right)$ on ${ℝ}^2\times {ℝ}_{+}\times \mathbb{S}$ which are continuously twice differentiable in x. For each $V\in C^{\mathrm{2,1}}\left({ℝ}^2\times {ℝ}_{+}\times \mathbb{S}\mathrm{<mo>\; ;\; </mo>}{ℝ}_{+}\right)$, we define an operator $LV$ by

$\begin{array}{c}LV\left(x\mathrm{,}t\mathrm{,}i\right)=V_t\left(x\mathrm{,}t\mathrm{,}i\right)+V_x\left(x\mathrm{,}t\mathrm{,}i\right)F\left(x\mathrm{,}t\mathrm{,}i\right)\\ +\frac{1}{2}trace\left[G^{\text{T}}\left(x\mathrm{,}t\mathrm{,}i\right)V_{xx}\left(x\mathrm{,}t\mathrm{,}i\right)G\left(x\mathrm{,}t\mathrm{,}i\right)\right]\end{array}$

$\begin{array}{c}+{\displaystyle {\int }_{\mathbb{Y}}}\left\{V\left(x+H\left(x\mathrm{,}t\mathrm{,}i\mathrm{,}v\right)\mathrm{,}t\mathrm{,}i\right)-V\left(x\mathrm{,}t\mathrm{,}i\right)-V_x\left(x\mathrm{,}t\mathrm{,}i\right)H\left(x\mathrm{,}t\mathrm{,}i\mathrm{,}v\right)\right\}\lambda \left(\text{d}u\right)\\ +{\displaystyle \underset{j=1}{\overset{N}{\sum }}}{q}_{ij}V\left(x\mathrm{,}t\mathrm{,}j\right)\mathrm{,}\end{array}$

where

$\begin{array}{l}{V}_t\left(x,t,i\right)=\frac{\partial V_x\left(x,t,i\right)}{\partial t},\\ V_x\left(x,t,i\right)=\left(\frac{\partial V_x\left(x,t,i\right)}{\partial x_1},\frac{\partial V_x\left(x,t,i\right)}{\partial x_2}\right),\\ V_{xx}\left(x,t,i\right)={\left(\frac{{\partial }^2{V}_x\left(x,t,i\right)}{\partial x_i\partial x_j}\right)}_{2\times 2}.\end{array}$

Then the generalized Itô’s formula with jumps is given by

$\begin{array}{c}\text{d}V\left(x\mathrm{,}t\mathrm{,}i\right)=LV\left(x\mathrm{,}t\mathrm{,}i\right)\text{d}t+V_x\left(x\mathrm{,}t\mathrm{,}i\right)G\left(x\mathrm{,}t\mathrm{,}i\right)\text{d}B\left(t\right)\\ +{\displaystyle {\int }_{\mathbb{Y}}}\left\{V\left(x+H\left(x\mathrm{,}t\mathrm{,}i\mathrm{,}{u}_i\right)\mathrm{,}t\mathrm{,}i\right)-V\left(x\mathrm{,}t\mathrm{,}i\right)\right\}\tilde{N}\left(\text{d}t\mathrm{,}\text{d}u\right)\mathrm{.}\end{array}$

For convenience and simplicity in the following discussion, define

$\{\begin{array}{l}\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 }\langle f\left(t\right)\rangle ,{\langle f\left(t\right)\rangle }_{*}=\underset{t\to +\infty }{\lim \inf }\langle f\left(t\right)\rangle ,\hfill \\ B_1\left(i\right)=a_1\left(i\right)-\frac{{\sigma }_1^2\left(i\right)}{2}-{\displaystyle {\int }_{\mathbb{Y}}}\left[{\gamma }_1\left(\mu ,i\right)-\ln \left(1+{\gamma }_1\left(\mu ,i\right)\right)\right]\lambda \left(\text{d}\mu \right),\hfill \\ B_2\left(i\right)=c_2\left(i\right)-a_2\left(i\right)-\frac{{\sigma }_2^2\left(i\right)}{2}-{\displaystyle {\int }_{\mathbb{Y}}}\left[{\gamma }_2\left(\mu ,i\right)-\ln \left(1+{\gamma }_2\left(\mu ,i\right)\right)\right]\lambda \left(\text{d}\mu \right),\hfill \\ B\left(i\right)={\min }_{j=1,2}\left(B_j\left(i\right)\right),\mathcal{B}={\displaystyle {\sum }_{i=1}^{\mathbb{S}}{\phi }_{i}B}\left(i\right),{\mathcal{B}}_j={\displaystyle {\sum }_{i=1}^{\mathbb{S}}{\phi }_i{B}_j}\left(i\right),\hfill \\ {\mathcal{C}}_j={\displaystyle {\sum }_{i=1}^{\mathbb{S}}{\phi }_i{c}_j}\left(i\right),\sigma ={\max }_{i\in \mathbb{S},j=1,2}\left\{\left|{\sigma }_j\left(i\right)\right|\right\},\hfill \\ \bar{b}\left(i\right)=\max \left\{b\left(i\right)\right\},\underset{\_}{b}\left(i\right)=\min \left\{b\left(i\right)\right\},\hfill \\ {\varphi }_1\left(i\right)=a_1\left(i\right)-\frac{1}{2}{\sigma }_1^2\left(i\right),{\varphi }_2\left(i\right)=-a_2\left(i\right)-\frac{1}{2}{\sigma }_2^2\left(i\right).\hfill \end{array}$ (2.6)

In this paper, we impose the following assumptions:

Assumption 1. There is a constant $c>0$ such that ${\displaystyle {\int }_{\mathbb{Y}}}{\left[\ln \left(1+{\gamma }_2\left(r\left(t\right)\mathrm{,}u\right)\right)\right]}^2\lambda \left(\text{d}u\right)\le c$.

Assumption 2. $a_j\left(i\right)$, $b_j\left(i\right)$, $e\left(i\right)$, $c_j\left(i\right)$ are all positive constant, and there exists ${\gamma }_j^{*}\ge {\gamma }_{j*}\left(i\right)>-1$ such that ${\gamma }_{j\mathrm{<mo>\; *\; </mo>}}\left(i\right)\le {\gamma }_j\left(u\mathrm{,}i\right)\le {\gamma }_j^{\mathrm{<mo>\; *\; </mo>}}\left(i\right)\left(u\in \mathbb{Y}\right)$, $\forall i\in \mathbb{S}$, $j=1,2$.

Assumption 3. For some $i\in \mathbb{S}$, $a_j\left(i\right)>0$, $b_j\left(i\right)>0$, $e\left(i\right)>0$, $c_j\left(i\right)>0$ and there exists ${\gamma }_j^{*}\ge {\gamma }_{j*}\left(i\right)>-1$ such that ${\gamma }_{j\mathrm{<mo>\; *\; </mo>}}\left(i\right)\le {\gamma }_j\left(u\mathrm{,}i\right)\le {\gamma }_j^{\mathrm{<mo>\; *\; </mo>}}\left(i\right)\left(u\in \mathbb{Y}\right)$, $j=1,2$.

Assumption 4. For some $j\in \mathbb{S}$, ${\gamma }_{ij}>0$, $\forall i\ne j$.

Assumption 5. $e\left(i\right)>0$, $\forall i\in \mathbb{S}$, and ${\displaystyle {\sum }_{i=1}^{\mathbb{S}}}{\phi }_i\beta \left(i\right)>0$, where

$\beta \left(i\right)=B\left(i\right)-\frac{c_1\left(i\right)}{e\left(i\right)}-\frac{{\left[{\sigma }_1\left(i\right)-{\sigma }_2\left(i\right)\right]}^2}{8}.$

Lemma 2.1. If Assumption 1 holds, then for any given initial data $\left(x\left(0\right)\mathrm{,}y\left(0\right)\right)\in {ℝ}_{+}^2$ , $r\left(0\right)\in \mathbb{S}$, system (6) has a unique solution $\left(x\left(t\right)\mathrm{,}y\left(t\right)\right)$ on $t\ge 0$ a.s..

Proof. The proof is standard (see e.g., [27] [28] ) and hence is omitted.

3. The Persistence and Extinction of Populations

Definition 3.1. Let $x\left(t\right)$ be a solution of system (1.6). Then

1) the species $x\left(t\right)$ is said to be extinct, if ${\lim }_{t\to \infty }x\left(t\right)=0$, a.s.;

2) the species $x\left(t\right)$ is said to be weakly persistent in the mean a.s., if $\lim {\sup }_{t\to +\infty }\frac{1}{t}{\displaystyle {\int }_0^t}x\left(t\right)\text{d}s>0$, a.s.;

3) the species $x\left(t\right)$ is said to non-persistence in mean a.s., if ${\lim }_{t\to +\infty }\frac{1}{t}{\displaystyle {\int }_0^t}x\left(t\right)\text{d}s=0$, a.s.

Lemma 3.1. (see Lipster [29] ). Suppose that $M\left(t\right)\mathrm{,}t\ge 0$ is a local martingale with $M\left(0\right)=0$. Then

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

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,$

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

Lemma 3.2. If Assumption 1 holds, then for any given initial value $\left(x\left(0\right)\mathrm{,}y\left(0\right)\mathrm{,}r\left(0\right)\right)\in {ℝ}_{+}^2\times \mathbb{S}$ , the solutio $\left(x\left(t\right)\mathrm{,}y\left(t\right)\right)$ to system (1.6) has the following property

${\langle \frac{\ln x\left(t\right)}{t}\rangle }^{\mathrm{<mo>\; *\; </mo>}}\le \mathrm{0,}{\langle \frac{\ln y\left(t\right)}{t}\rangle }^{\mathrm{<mo>\; *\; </mo>}}\le 0.$ (3.1)

Proof. Applying the generalized Itô’s formula to ${\text{e}}^t\ln x$ and ${\text{e}}^t\ln y$ leads to

$\begin{array}{l}{\text{e}}^t\ln x\left(t\right)-\ln x\left(0\right)\\ ={\displaystyle {\int }_0^t}{\text{e}}^s[\ln x(s)+a_1\left(r\left(s\right)\right)-\frac{{\sigma }_1^2\left(r\left(s\right)\right)}{2}-\frac{c_1\left(r\left(s\right)\right)y\left(s\right)}{x\left(s\right)+e\left(r\left(s\right)\right)y\left(s\right)}\\ -{\displaystyle {\int }_Y}\left[{\gamma }_1\left(u,r\left(s\right)\right)-\ln \left(1+{\gamma }_1\left(u,r\left(s\right)\right)\right)\right]\lambda \left(\text{d}u\right)-b_1\left(r\left(s\right)\right)x\left(s\right)]\text{d}s+N_{1j}\left(t\right),\end{array}$

$\begin{array}{l}{\text{e}}^t\ln y\left(t\right)-\ln y\left(0\right)\\ ={\displaystyle {\int }_0^t}{\text{e}}^s[\ln y\left(s\right)-a_2\left(r\left(s\right)\right)-\frac{{\sigma }_2^2\left(r\left(s\right)\right)}{2}+\frac{c_2\left(r\left(s\right)\right)x\left(s\right)}{x\left(s\right)+e\left(r\left(s\right)\right)y\left(s\right)}\\ -{\displaystyle {\int }_Y}\left[{\gamma }_2\left(u,r\left(s\right)\right)-\ln \left(1+{\gamma }_2\left(u,r\left(s\right)\right)\right)\right]\lambda \left(\text{d}u\right)-b_2\left(r\left(s\right)\right)y\left(s\right)]\text{d}s+N_{2j}\left(t\right),\end{array}$ (3.2)

where

$\{\begin{array}{l}{N}_{j1}\left(t\right):={\displaystyle {\int }_0^t}{\text{e}}^s{\sigma }_j\left(r\left(s\right)\right)\text{d}B\left(s\right),\\ N_{j2}\left(t\right):={\displaystyle {\int }_0^t}{\displaystyle {\int }_{\mathbb{Y}}}{\text{e}}^s\ln \left(1+{\gamma }_j\left(u,r\left(t\right)\right)\right)\tilde{N}\left(\text{d}s,\text{d}u\right),j=1,2,\end{array}$

are local martingale with the quadratic variations

$\{\begin{array}{l}\langle N_{j1}\rangle \left(t\right):={\displaystyle {\int }_0^t}{\text{e}}^s{\sigma }_j\left(r\left(s\right)\right)\text{d}B\left(s\right)\le \frac{{\sigma }_1^2}{2}\left({\text{e}}^{2t}-1\right)\text{a}\text{.s}.\\ \langle N_{j2}\rangle \left(t\right):={\displaystyle {\int }_0^t}{\displaystyle {\int }_{\mathbb{Y}}}{\text{e}}^s\ln \left(1+{\gamma }_j\left(u,r\left(t\right)\right)\right)\tilde{N}\left(\text{d}s,\text{d}u\right)\le \frac{c}{2}\left({\text{e}}^{2t}-1\right)\text{a}\text{.s}.\end{array}$

In view of the exponential martingale inequality, for any positive number $\alpha ,\beta ,T$,

$P\left\{\underset{0\le t\le T}{\sup }\left[N_j\left(t\right)-\frac{\alpha }{2}\langle N_j\rangle \left(t\right)\right]>\beta \right\}\le {\text{e}}^{-\alpha \beta }.$

Let $T=\gamma k$, $\alpha ={\text{e}}^{-\gamma k}$, $\beta =\xi {\text{e}}^{\gamma k}\ln k$, we can get

$P\left\{\underset{0\le t\le \gamma k}{\sup }\left[N_j\left(t\right)-\frac{{\text{e}}^{-\gamma k}}{2}\langle N_j\rangle \left(t\right)\right]>\xi {\text{e}}^{\gamma k}\ln k\right\}\le \frac{1}{k^{\xi }}\left(j=1,2\right),$

where $\gamma >1$ and $\xi >1$. Applying the Borel-Cantelli Lemma, we can obtain that for all $\omega \in \Omega$, there is a random integer $k_0=k_0\left(\omega \right)$ such that for each $k>k_0\left(\omega \right)$,

$\underset{0\le t\le \gamma k}{\sup }\left[N_j\left(t\right)-\frac{{\text{e}}^{-\gamma k}}{2}\langle N_j\rangle \left(t\right)\right]\le \xi {\text{e}}^{\gamma k}\ln k.$

Namely, we have shown

$N_j\left(t\right)\le \frac{1}{2}{\text{e}}^{-\gamma k}\langle N_j\left(t\right)\rangle +\xi {\text{e}}^{\gamma k}\ln k,0\le t\le \gamma k\left(j=1,2\right).$

Substituting the above inequalities into (3.2) leads to

$\begin{array}{l}{\text{e}}^t\ln x\left(t\right)-\ln x\left(0\right)\\ \le {\displaystyle {\int }_0^t}{\text{e}}^s[\ln x\left(s\right)+a_1\left(r\left(s\right)\right)-\frac{{\sigma }_1^2\left(r\left(s\right)\right)}{2}-\frac{c_1\left(r\left(s\right)\right)y\left(s\right)}{x\left(s\right)+e\left(r\left(s\right)\right)y\left(s\right)}\\ -{\displaystyle {\int }_{\mathbb{Y}}}\left[{\gamma }_1\left(u,r\left(s\right)\right)-\ln \left(1+{\gamma }_1\left(u,r\left(s\right)\right)\right)\right]\lambda \left(\text{d}u\right)-b_1\left(r\left(s\right)\right)x\left(s\right)]\text{d}s\\ +\frac{1}{2}{\text{e}}^{-\gamma k}{\displaystyle {\int }_0^t}{\text{e}}^{2s}{\sigma }_1^2\left(r\left(s\right)\right)\text{d}s+\frac{1}{2}c{\text{e}}^{-\gamma k}{\displaystyle {\int }_0^t}{\text{e}}^{2s}\text{d}s+2\xi {\text{e}}^{\gamma k}\ln k\end{array}$

$\begin{array}{l}={\displaystyle {\int }_0^t}{\text{e}}^s[\ln x\left(s\right)+a_1\left(r\left(s\right)\right)-\frac{{\sigma }_1^2\left(r\left(s\right)\right)}{2}-\frac{c_1\left(r\left(s\right)\right)y\left(s\right)}{x\left(s\right)+e\left(r\left(s\right)\right)y\left(s\right)}\\ -{\displaystyle {\int }_{\mathbb{Y}}}\left[{\gamma }_1\left(u,r\left(s\right)\right)-\ln \left(1+{\gamma }_1\left(u,r\left(s\right)\right)\right)\right]\lambda \left(\text{d}u\right)+\frac{1}{2}{\text{e}}^{s-\gamma k}{\sigma }_1^2\left(r\left(s\right)\right)\\ -\frac{1}{2}c\left(1-{\text{e}}^{s-\gamma k}\right)-b_1\left(r\left(s\right)\right)x\left(s\right)]\text{d}s+2\xi {\text{e}}^{\gamma k}\ln k\\ \le {\displaystyle {\int }_0^t}{\text{e}}^s[\ln x\left(s\right)+{\bar{a}}_1\left(i\right)+{\displaystyle {\int }_{\mathbb{Y}}}\left[\left|{\bar{\gamma }}_1\left(u\right)\right|+\left|\ln \left(1+{\bar{\gamma }}_1\left(u\right)\right)\right|\right]\lambda \left(\text{d}u\right)\\ -\left(1-{\text{e}}^{s-\gamma k}\right)\frac{1}{2}{\sigma }_1^2\left(i\right)-\frac{1}{2}c\left(1-{\text{e}}^{s-\gamma k}\right)]\text{d}s+2\xi {\text{e}}^{\gamma k}\ln k,\end{array}$

$\begin{array}{l}{\text{e}}^t\ln y\left(t\right)-\ln y\left(0\right)\\ \le {\displaystyle {\int }_0^t}{\text{e}}^s[\ln y\left(s\right)-a_2\left(r\left(s\right)\right)-\frac{{\sigma }_2^2\left(r\left(s\right)\right)}{2}+\frac{c_2\left(r\left(s\right)\right)x\left(s\right)}{x\left(s\right)+e\left(r\left(s\right)\right)y(\; s\; )}\end{array}$

$\begin{array}{l}-{\displaystyle {\int }_{\mathbb{Y}}}\left[{\gamma }_2\left(u,r\left(s\right)\right)-\ln \left(1+{\gamma }_2\left(u,r\left(s\right)\right)\right)\right]\lambda \left(\text{d}u\right)-b_2\left(r\left(s\right)\right)y\left(s\right)]\text{d}s\\ +\frac{1}{2}{\text{e}}^{-\gamma k}{\displaystyle {\int }_0^t}{\text{e}}^{2s}{\sigma }_2^2\left(r\left(s\right)\right)\text{d}s+\frac{1}{2}c{\text{e}}^{-\gamma k}{\displaystyle {\int }_0^t}{\text{e}}^{2s}\text{d}s+2\xi {\text{e}}^{\gamma k}\ln k\end{array}$

$\begin{array}{l}={\displaystyle {\int }_0^t}{\text{e}}^s[\ln y\left(s\right)-a_2\left(r\left(s\right)\right)-\frac{{\sigma }_2^2\left(r\left(s\right)\right)}{2}+\frac{c_2\left(r\left(s\right)\right)x\left(s\right)}{x\left(s\right)+e\left(r\left(s\right)\right)y\left(s\right)}\\ -{\displaystyle {\int }_{\mathbb{Y}}}\left[{\gamma }_2\left(u,r\left(s\right)\right)-\ln \left(1+{\gamma }_2\left(u,r\left(s\right)\right)\right)\right]\lambda \left(\text{d}u\right)+\frac{1}{2}{\text{e}}^{s-\gamma k}{\sigma }_2^2\left(r\left(s\right)\right)\\ -\frac{1}{2}c\left(1-{\text{e}}^{s-\gamma k}\right)-b_2\left(r\left(s\right)\right)y\left(s\right)]\text{d}s+2\xi {\text{e}}^{\gamma k}\ln k\\ \le {\displaystyle {\int }_0^t}{\text{e}}^s[\ln y\left(s\right)-a_2\left(i\right)+{\displaystyle {\int }_{\mathbb{Y}}}\left[\left|{\bar{\gamma }}_2\left(u\right)\right|+\left|\ln \left(1+{\bar{\gamma }}_2\left(u\right)\right)\right|\right]\lambda \left(\text{d}u\right)\\ +{\bar{c}}_2\left(i\right)-\frac{1}{2}\left(1-{\text{e}}^{s-\gamma k}\right){\sigma }_1^2\left(i\right)-\frac{1}{2}c\left(1-{\text{e}}^{s-\gamma k}\right)]\text{d}s+2\xi {\text{e}}^{\gamma k}\ln k.\end{array}$

It is easy to find that for any $0\le s\le \gamma k$ and $x>0,y>0$, there exists two constants $C_1$ and $C_2$ which are independent of k such that

$\begin{array}{l}\ln x\left(s\right)+{\bar{a}}_1\left(s\right)+{\displaystyle {\int }_{\mathbb{Y}}}\left(\left|{\bar{\gamma }}_1\left(u\right)\right|+\left|\ln \left(1+{\bar{\gamma }}_1\left(u\right)\right)\right|\right)\lambda \left(\text{d}u\right)\\ -\frac{1}{2}{\sigma }_1^2\left(1-{\text{e}}^{s-\gamma k}\right)-\frac{1}{2}c\left(1-{\text{e}}^{s-\gamma k}\right)\le C_1,\end{array}$

$\begin{array}{l}\ln y\left(s\right)-a_2\left(s\right)+{\bar{c}}_2\left(s\right)+{\displaystyle {\int }_{\mathbb{Y}}}\left(\left|{\bar{\gamma }}_2\left(u\right)\right|+\left|\ln \left(1+{\bar{\gamma }}_2\left(u\right)\right)\right|\right)\lambda \left(\text{d}u\right)\\ -\frac{1}{2}{\sigma }_2^2\left(1-{\text{e}}^{s-\gamma k}\right)-\frac{1}{2}c\left(1-{\text{e}}^{s-\gamma k}\right)\le C_2\mathrm{.}\end{array}$

In other words, for any $0\le t\le \gamma k,k>k_0$, we obtain

$\ln x\left(t\right)\le {\text{e}}^{-t}\ln x\left(0\right)+C_1\left(1-{\text{e}}^{-t}\right)+2\theta {\text{e}}^{-t}{\text{e}}^{\gamma k}\ln k\mathrm{,}$

$\ln y\left(t\right)\le {\text{e}}^{-t}\ln y\left(0\right)+C_2\left(1-{\text{e}}^{-t}\right)+2\theta {\text{e}}^{-t}{\text{e}}^{\gamma k}\ln k\mathrm{.}$

Letting $t\to \infty$, we get

${\langle \frac{\ln x\left(t\right)}{t}\rangle }^{\mathrm{<mo>\; *\; </mo>}}\le \mathrm{0,}{\langle \frac{\ln y\left(t\right)}{t}\rangle }^{\mathrm{<mo>\; *\; </mo>}}\le 0.$

This proof is completed.

Theorem 3.1. Let Assumption 1 hold, if ${\mathcal{B}}_1<0$ and ${\mathcal{B}}_2<0$ , then species $x\left(t\right)$ and $y\left(t\right)$ of system (1.6) will tend to be extinct.

Proof. By the generalized Itô’s formula, we derive from (1.6) that

$\begin{array}{l}\text{d}\ln x\left(t\right)=[a_1\left(r\left(t\right)\right)-\frac{{\sigma }_1^2\left(r\left(t\right)\right)}{2}-{\displaystyle {\int }_{\mathbb{Y}}}\left[{\gamma }_1\left(\mu ,r\left(t\right)\right)-\ln \left(1+{\gamma }_1\left(\mu ,r\left(t\right)\right)\right)\right]\lambda \left(\text{d}\mu \right)\\ -\frac{c_1\left(r\left(t\right)\right)y\left(t\right)}{x\left(t\right)+e\left(r\left(t\right)\right)y\left(t\right)}-b_1\left(r\left(t\right)\right)x\left(t\right)]\text{d}t\\ +{\sigma }_1\left(r\left(t\right)\right)\text{d}B\left(t\right)+{\displaystyle {\int }_{\mathbb{Y}}}\ln \left[1+{\gamma }_1\left(\mu ,r\left(t\right)\right)\right]\tilde{N}\left(\text{d}t,\text{d}\mu \right),\\ \text{d}\ln y\left(t\right)=[-a_2\left(r\left(t\right)\right)-\frac{{\sigma }_2^2\left(r\left(t\right)\right)}{2}-{\displaystyle {\int }_{\mathbb{Y}}}\left[{\gamma }_2\left(\mu ,r\left(t\right)\right)-\ln \left(1+{\gamma }_2\left(\mu ,r\left(t\right)\right)\right)\right]\lambda \left(\text{d}\mu \right)\\ +\frac{c_2\left(r\left(t\right)\right)x\left(t\right)}{x\left(t\right)+e\left(r\left(t\right)\right)y\left(t\right)}-b_2\left(r\left(t\right)\right)y\left(t\right)]\text{d}t\\ +{\sigma }_2\left(r\left(t\right)\right)\text{d}B\left(t\right)+{\displaystyle {\int }_{\mathbb{Y}}}\ln \left[1+{\gamma }_2\left(\mu ,r\left(t\right)\right)\right]\tilde{N}\left(\text{d}t,\text{d}\mu \right).\end{array}$ (3.3)

Then

$\begin{array}{l}\ln x\left(t\right)-\ln x\left(0\right)\\ ={\displaystyle {\int }_0^t}[\ln x\left(s\right)+a_1\left(r\left(s\right)\right)-\frac{{\sigma }_1^2\left(r\left(s\right)\right)}{2}-{\displaystyle {\int }_{\mathbb{Y}}}\left[{\gamma }_1\left(u\mathrm{,}r\left(s\right)\right)-\ln \left(1+{\gamma }_1\left(u\mathrm{,}r\left(s\right)\right)\right)\right]\lambda \left(\text{d}u\right)\\ -\frac{c_1\left(r\left(s\right)\right)y\left(s\right)}{x\left(s\right)+e\left(r\left(s\right)\right)y\left(s\right)}-b_1\left(r\left(s\right)\right)x\left(s\right)]\text{d}s+M_{1j}\left(t\right)\mathrm{,}\\ \ln y\left(t\right)-\ln y\left(0\right)\\ ={\displaystyle {\int }_0^t}[\ln y\left(s\right)-a_2\left(r\left(s\right)\right)-\frac{{\sigma }_2^2\left(r\left(s\right)\right)}{2}-{\displaystyle {\int }_{\mathbb{Y}}}\left[{\gamma }_2\left(u\mathrm{,}r\left(s\right)\right)-\ln \left(1+{\gamma }_2\left(u\mathrm{,}r\left(s\right)\right)\right)\right]\lambda \left(\text{d}u\right)\\ +\frac{c_2\left(r\left(s\right)\right)x\left(s\right)}{x\left(s\right)+e\left(r\left(s\right)\right)y\left(s\right)}\text{d}s-b_2\left(r\left(s\right)\right)y\left(s\right)]+M_{2j}\left(t\right)\mathrm{,}\end{array}$ (3.4)

where

$\{\begin{array}{l}{M}_{j1}\left(t\right):={\displaystyle {\int }_0^t}{\sigma }_j\left(r\left(s\right)\right)\text{d}B\left(s\right),\\ M_{j2}\left(t\right):={\displaystyle {\int }_0^t}{\displaystyle {\int }_{\mathbb{Y}}}\ln \left(1+{\gamma }_j\left(u,r\left(t\right)\right)\right)\tilde{N}\left(\text{d}s,\text{d}u\right),\end{array}$