Analysis of a Stochastic Ratio-Dependent Predator-Prey System with Markovian Switching and Lévy Jumps

Xuegui Zhang, Yuanfu Shao^*, Taolin Zhang

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

DOI: 10.4236/jamp.2020.811195                

Abstract

In this paper, the dynamics of a stochastic ratio-dependent predator-prey system with markovian switching and Lévy noise is studied. Firstly, we show the existence condition of global positive solution under the given positive initial value. Secondly, sufficient conditions for system extinction and persistence are obtained through some assumptions. Then, the sufficient conditions of stochastically persistence are obtained by combining stochastic analysis technique and M-matrix analysis. In addition, under appropriate conditions, we demonstrate the existence of a unique stationary distribution for a system without Lévy jumps. Finally, the empirical and Mlistein methods are used to verify the theoretical results through numerical simulation.

Keywords

Ratio-Dependent Predator-Prey System, Lévy, Stochastically Permanence, Ergodicity

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

are local martingale with the quadratic variations

$\{\begin{array}{l}\langle M_{j1}\rangle \left(t\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle {\int }_0^t}{\sigma }_j\left(r\left(s\right)\right)\text{d}B\left(s\right)\le {\sigma }_1^{2}t\text{a}\text{.s}\hfill \\ \langle M_{j2}\rangle \left(t\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle {\int }_0^t}{\displaystyle {\int }_{\text{Y}}}\ln \left(1+{\gamma }_j\left(u\mathrm{,}r\left(t\right)\right)\right)\tilde{N}\left(\text{d}s\mathrm{,}\text{d}u\right)\le ct\text{a}\text{.s}\text{.}\hfill \end{array}$

It follows from Lemma 3.1 that

$\underset{x\to +\infty }{\lim }\frac{M_j\left(t\right)}{t}=\mathrm{0,}j=\mathrm{1,2,}\text{a}\text{.s}\mathrm{.}$ (3.5)

Dividing by t on both sides of (18) and then taking the superior limit, we obtain

${\langle \frac{\ln x\left(t\right)}{t}\rangle }^{*}\le {\langle B_1\left(r\left(t\right)\right)\rangle }^{*}={\mathcal{B}}_1<0,$ (3.6)

${\langle \frac{\ln y\left(t\right)}{t}\rangle }^{*}\le {\langle B_2\left(r\left(t\right)\right)\rangle }^{*}={\mathcal{B}}_2<0,$ (3.7)

that is

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

This proof is completed.

Theorem 3.2. 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 be non-persistent in the mean.

Proof. Due to the fact that ${{\displaystyle \lim \sup }}_{t\to \infty }\frac{1}{t}{\displaystyle {\int }_0^t}{B}_1\left(r\left(t\right)\right)\text{d}s={\mathcal{B}}_1$ and (3.5), for any given $\epsilon >0$, there is a positive constant T such that

$\frac{{\displaystyle {\int }_0^t}{B}_1\left(r\left(t\right)\right)\text{d}s}{t}\le {\mathcal{B}}_1+\frac{\epsilon }{2}=\frac{\epsilon }{2},t>T.$

Then for any $\epsilon >0$ and sufficiently large $t>T$, we have

$\begin{array}{c}\ln x\left(t\right)-\ln x\left(0\right)\le {\displaystyle {\int }_0^t}\left[B_1\left(r\left(s\right)\right)-b_1\left(r\left(s\right)\right)x\left(s\right)\right]\text{d}s+M_{1j}\left(t\right)\\ \le \epsilon t-b_1{\displaystyle {\int }_0^t}x\left(s\right)\text{d}s.\end{array}$

Let $u\left(t\right)={\displaystyle {\int }_0^t}x\left(s\right)\text{d}s$ for $t>T$, then one gets

$\ln \left(\frac{\text{d}u}{\text{d}t}\right)\le \epsilon t-b_{1}u\left(t\right)+\ln x\left(0\right),$

which shows that

${\text{e}}^{b_{1}u\left(t\right)}\text{d}u\le x\left(0\right){\text{e}}^{\epsilon t}\text{d}t,t>T.$

Integrating the above inequality from T to t results in

$\frac{1}{{\hat{b}}_1}\left({\text{e}}^{b_{1}u\left(t\right)}-{\text{e}}^{b_{1}u\left(T\right)}\right)\le \frac{x\left(0\right)}{\epsilon }\left({\text{e}}^{\epsilon t}-{\text{e}}^{\epsilon T}\right),$

which implies

${\text{e}}^{b_{1}u\left(t\right)}<{\text{e}}^{b_{1}u\left(T\right)}+\frac{x\left(0\right)b_1}{\epsilon }{\text{e}}^{\epsilon t}-\frac{x\left(0\right)b_1}{\epsilon }{\text{e}}^{\epsilon T}.$

Taking the logarithm on both sides, we have

$u\left(t\right)<\frac{1}{b_1}\ln \left[{\text{e}}^{b_{1}u\left(T\right)}+\frac{x\left(0\right)b_1}{\epsilon }{\text{e}}^{\epsilon t}-\frac{x\left(0\right)b_1}{\epsilon }{\text{e}}^{\epsilon T}\right].$

Then letting $t\to \infty$, and making use of the L’Hospital’s rule results in

${\langle x\left(t\right)\rangle }^{*}=\underset{t\to \infty }{\lim \sup }\frac{1}{t}{\displaystyle {\int }_0^t}x\left(s\right)\text{d}s\le \frac{1}{b_1}\underset{t\to \infty }{\lim \sup }\frac{1}{t}\ln b_{1}x\left(0\right){\epsilon }^{-1}{\text{e}}^{\epsilon t}=\frac{\epsilon }{b_1},\text{a}\text{.s}\text{.}$

By the arbitrariness of $\epsilon$, we can get

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

which is the desired assertion for $x\left(t\right)$. The corresponding result of $y\left(t\right)$ can also be proved by the same method, so it is omitted. Thus the proof of Theorem 3.2 is completed.

Theorem 3.3. Let Assumption 1 hold, if ${\mathcal{B}}_1>{\mathcal{C}}_1$ , then species $x\left(t\right)$ of system (1.6) will be weakly persistent in the mean.

Proof. By the generalized Itô’s formula, and dividing by t on both sides, we have

$\begin{array}{l}\frac{\ln x\left(t\right)-\ln x\left(0\right)}{t}\\ =\frac{1}{t}{\displaystyle {\int }_0^t}\left[B_1\left(r\left(s\right)\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(\xi \left(s\right)\right)x\left(s\right)\right]\text{d}s+\frac{M_{j1}\left(t\right)}{t}\\ \ge \frac{1}{t}{\displaystyle {\int }_0^t}\left[B_1\left(r\left(s\right)\right)-\left(c_1\left(r\left(s\right)\right)-b_1\left(r\left(s\right)\right)\right)x\left(s\right)\right]\text{d}s+\frac{M_{j1}\left(t\right)}{t},\end{array}$

Taking the superior limit and combining with (3.1), (3.5), we conclude

$\underset{t\to \infty }{\lim \sup }\frac{\ln x\left(t\right)}{t}\ge {\langle B_1\left(r\left(t\right)\right)\rangle }^{*}-{\langle c_1\left(r\left(t\right)\right)\rangle }_{*}={\mathcal{B}}_1-{\mathcal{C}}_1>0,\text{a}\text{.s}.$

The proof is completed.

Corollary 3.1. According to Theorem 3.1, the solutions of subsystems (1.8) can be obtained with the following properties:

1) If $B_1\left(i\right)<0$, $B_2\left(i\right)<0$, then both $x\left(t\right)$ and $y\left(t\right)$ tend to extinction a.s., i.e., ${\lim }_{t\to +\infty }x\left(t\right)=0$ and ${\lim }_{t\to +\infty }y\left(t\right)=0$ ;

2) If $B_1\left(i\right)=0$, $B_2\left(i\right)=0$, then both $x\left(t\right)$ and $y\left(t\right)$ will be non-persistent in the mean a.s.;

3) If $B_1\left(i\right)>c_1\left(i\right)$, then $x\left(t\right)$ is weakly persistent in the mean a.s.

4. Stochastic Permanence

In the study of population dynamics, stochastic permanence is one of the most important properties. We first introduce the definition of stochastic permanence [6], which is widely used in the field of population dynamics (see e.g., [22] [28] ). In this section, we shall discuss this property.

Definition 4.1. System (1.6) is stochastically permanent if for any $ϵ>0$ , there exist ${\delta }_{*}={\delta }_{*}\left(ϵ\right)>0$ and ${\delta }^{*}={\delta }^{*}\left(ϵ\right)>0$ such that

$\underset{t\to \infty }{\lim \inf }P\left\{x\left(t\right)\ge {\delta }_{*}\right\}\ge 1-\epsilon ,\underset{t\to \infty }{\lim \inf }P\left\{y\left(t\right)\ge {\delta }_{*}\right\}\ge 1-\epsilon ,$ (4.1)

$\underset{t\to \infty }{\lim \inf }P\left\{x\left(t\right)\le {\delta }^{\mathrm{<mo>\; *\; </mo>}}\right\}\ge 1-\epsilon \mathrm{,}\underset{t\to \infty }{\lim \inf }P\left\{y\left(t\right)\le {\delta }^{\mathrm{<mo>\; *\; </mo>}}\right\}\ge 1-\epsilon \mathrm{.}$ (4.2)

Lemma 4.1. Suppose Assumption 2 hold, let $\left(x\left(t\right)\mathrm{,}y\left(t\right)\right)$ be the solution to system (1.6) with initial value $\left(\left(x\left(0\right)\mathrm{,}y\left(0\right)\right)\mathrm{,}r\left(0\right)\right)\in {ℝ}_{+}^2\times \mathbb{S}$ , then for any $\left({\theta }_1\mathrm{,}{\theta }_2\right)\in {ℝ}_{+}^2$ , there exists $K\left({\theta }_1,{\theta }_2\right)>0$ such that

$\underset{t\to \infty }{\lim \sup }E\left[x^{{\theta }_1}\left(t\right)+y^{{\theta }_2}\left(t\right)\right]\le K\left({\theta }_1,{\theta }_2\right).$ (4.3)

Proof. Define $W\left(x,y,i\right)=x^{{\theta }_1}+y^{{\theta }_2}$. By the generalized Itô’s formula, we deduce

$\begin{array}{c}L\left[W\left(x,y,i\right)\right]=\frac{1}{2}{\theta }_1\left({\theta }_1-1\right){\sigma }_1^2\left(i\right)x^{{\theta }_1}+\frac{1}{2}{\theta }_2\left({\theta }_2-1\right){\sigma }_2^2\left(i\right)y^{{\theta }_2}\\ +{\theta }_1{x}^{{\theta }_1}\left[a_1\left(i\right)-b_1\left(i\right)x-\frac{c_1\left(i\right)y}{x+e\left(i\right)y}\right]\\ +{\theta }_2{y}^{{\theta }_2}\left[-a_2\left(i\right)-b_2\left(i\right)y+\frac{c_2\left(i\right)x}{x+e\left(i\right)y}\right]\\ +x^{{\theta }_1}{\displaystyle {\int }_{\mathbb{Y}}}\left[{\left(1+{\gamma }_1\left(t,u\right)\right)}^{{\theta }_1}-1-{\theta }_1{\gamma }_1\left(t,u\right)\right]\lambda \left(\text{d}u\right)\\ +y^{{\theta }_2}{\displaystyle {\int }_{\mathbb{Y}}}\left[{\left(1+{\gamma }_2\left(t,u\right)\right)}^{{\theta }_2}-1-{\theta }_2{\gamma }_2\left(t,u\right)\right]\lambda \left(\text{d}u\right),\end{array}$ (4.4)

dropping t from $x\left(t\right)$ and $y\left(t\right)$. We notice that both coefficients of the higher $x^{{\theta }_1+1}$, $y^{{\theta }_2+1}$ are negative. From (4.4) and ${\underset{\_}{b}}_j(i)>0,i=1,2$, we derive that there exists a positive constant $K\left({\theta }_1\mathrm{,}{\theta }_2\right)$ such that

$\begin{array}{l}L\left[W\left(x,y,i\right)\right]+W\left(x,y,i\right)\\ =x^{{\theta }_1}\{1+\frac{1}{2}{\theta }_1\left({\theta }_1-1\right){\sigma }_1^2\left(i\right)+{\theta }_1\left[a_1\left(i\right)-b_1\left(i\right)x-\frac{c_1\left(i\right)y}{x+e\left(i\right)y}\right]\end{array}$

$\begin{array}{l}+{\displaystyle {\int }_{\mathbb{Y}}}\left[{\left(1+{\gamma }_1\left(t,u\right)\right)}^{{\theta }_1}-1-{\theta }_1{\gamma }_1\left(t,u\right)\right]\lambda \left(\text{d}u\right)\}\\ +y^{{\theta }_2}\{1+\frac{1}{2}{\theta }_2\left({\theta }_2-1\right){\sigma }_2^2\left(i\right)+{\theta }_2\left[-a_2\left(i\right)-b_2\left(i\right)y+\frac{c_2\left(i\right)x}{x+e\left(i\right)y}\right]\\ +{\displaystyle {\int }_{\mathbb{Y}}}\left[{\left(1+{\gamma }_2\left(t,u\right)\right)}^{{\theta }_2}-1-{\theta }_2{\gamma }_2\left(t,u\right)\right]\lambda \left(\text{d}u\right)\}\end{array}$

$\begin{array}{l}\le x^{{\theta }_1}\{1+\frac{1}{2}{\theta }_1\left({\theta }_1-1\right){\sigma }_1^2\left(i\right)+{\theta }_1\left[a_1\left(i\right)-b_1\left(i\right)x\right]\\ \begin{array}{c}\\ \end{array}+{\displaystyle {\int }_{\mathbb{Y}}}\left[{\left(1+{\gamma }_1\left(t,u\right)\right)}^{{\theta }_1}-1-{\theta }_1{\gamma }_1\left(t,u\right)\right]\lambda \left(\text{d}u\right)\}\\ +y^{{\theta }_2}\{1+\frac{1}{2}{\theta }_2\left({\theta }_2-1\right){\sigma }_2^2\left(i\right)+{\theta }_2\left[-a_2\left(i\right)-b_2\left(i\right)y+c_2\left(i\right)\right]\\ \begin{array}{c}\\ \end{array}+{\displaystyle {\int }_{\mathbb{Y}}}\left[{\left(1+{\gamma }_2\left(t,u\right)\right)}^{{\theta }_2}-1-{\theta }_2{\gamma }_2\left(t,u\right)\right]\lambda \left(\text{d}u\right)\}\\ \le {\displaystyle \underset{i=1}{\overset{2}{\sum }}}{K}_i\left({\theta }_i\right)=:K\left({\theta }_1,{\theta }_2\right).\end{array}$ (4.5)

According to (4.5) and the generalized Itô’s formula, we obtain

$L\left[{\text{e}}^{t}W\left(x,y,i\right)\right]={\text{e}}^{t}W\left(x,y,i\right)+{\text{e}}^{t}L\left[W\left(x,y,i\right)\right]\le {\text{e}}^{t}K\left({\theta }_1,{\theta }_2\right).$

Integrating $d\left({\text{e}}^{t}W\left(x\left(t\right)\mathrm{,}y\left(t\right)\mathrm{,}r\left(t\right)\right)\right)$ from 0 to t and taking expectations of both sides, we obtain that

${\text{e}}^{t}E\left[x^{{\theta }_1}\left(t\right)+y^{{\theta }_2}\left(t\right)\right]\le x^{{\theta }_1}\left(0\right)+y^{{\theta }_2}\left(0\right)+{\text{e}}^{t}K\left({\theta }_1\mathrm{,}{\theta }_2\right)\mathrm{,}$

which implies the required assertion (4.3).

To state our main result, we give some notations. Let C be vector or matrix. By $C\gg 0$ we mean all elements are positive. Let

$Z^{N\times N}=\left\{C={\left(c_{ij}\right)}_{S\times S}:c_{ij}\le 0,\forall i.j\in \mathbb{S},i\ne j\right\}.$

Lemma 4.2. (see Lemma 5.3 in [18] ). If $C=\left(c_{ij}\right)\in Z^{S\times S}$ has all of its row sums positive, that is

${\displaystyle \underset{j=1}{\overset{S}{\sum }}}{c}_{ij}>0$ for all $1\le i\le S$.

Then $\det A>0$.

Lemma 4.3. (see Lemma 5.3 in [18] ). If $C\in Z^{S\times S}$ , then the following statements are equivalent:

1) C is nonsingular M-matrix.

2) All of the principal minors of are positive; that is

$\left|\begin{array}{cccc}{c}_{11}& c_{12}& \cdots & c_{1k}\\ c_{21}& c_{22}& \cdots & c_{2k}\\ ⋮& ⋮& \ddots & ⋮\\ c_{k1}& c_{k2}& \cdots & c_{kk}\end{array}\right|>0$ for every $k\in \mathbb{S}$.

3) C is semi-positive; that is, there exists $x\gg 0$ in ${ℝ}^N$ such that $Cx\gg 0$.

The proof of the stochastic permanence is rather technical, so we prepare several useful lemmas.

Lemma 4.4. Assumptions 3 and 4 imply that there is a constant ${\theta }_0>0$ such that for $0<\theta <{\theta }_0$ , the matrix

$G\left(\theta \right)=\text{diag}\left({\xi }_1\left(\theta \right)\mathrm{,}{\xi }_2\left(\theta \right)\mathrm{,}\cdots \mathrm{,}{\xi }_S\left(\theta \right)\right)-\Gamma \mathrm{,}$

is a nonsingular M-matrix, where

$\begin{array}{c}{\xi }_i\left(\theta \right)=\theta \beta \left(i\right)-\frac{1}{2}{\theta }^2{\sigma }^2-\theta {\displaystyle {\int }_{\mathbb{Y}}}\ln \left(1+\underset{j=1,2}{\max }\left\{{\gamma }_j^{*}\left(i\right)\right\}\right)\lambda \left(\text{d}u\right)\\ -{\displaystyle {\int }_{\mathbb{Y}}}\left[{\left(1+\underset{j=1,2}{\min }\left\{{\gamma }_{j\mathrm{<mo>\; *\; </mo>}}\left(i\right)\right\}\right)}^{-\theta }-1\right]\lambda \left(\text{d}u\right)\mathrm{,}\end{array}$

and $\beta \left(i\right)$ is defined in Assumption 5.

Proof. The proof is rather standard and hence is omitted (see [15] ).

Lemma 4.5. If $e\left(i\right)>0$ and there exists a constant $\theta >0$ such that $G\left(\theta \right)$ is a nonsingular M-matrix, then the solution $\left(x\left(t\right)\mathrm{,}y\left(t\right)\right)$ of system (1.6) with initial value $\left(\left(x\left(0\right)\mathrm{,}y\left(0\right)\right)\mathrm{,}r\left(0\right)\right)\in {ℝ}_{+}^2\times \mathbb{S}$ has the property that

$\underset{t\to \infty }{\lim \sup }E\left(\frac{1}{x^{\theta }\left(t\right)}\right)\le H_1,\underset{t\to \infty }{\lim \sup }E\left(\frac{1}{y^{\theta }\left(t\right)}\right)\le H_2.$

Proof. Define

$\begin{array}{l}{U}_1=\frac{1}{x},U_2=\frac{1}{y},U=U_1+U_2.\hfill \end{array}$ (4.6)

Applying the generalized Itô’s formula, we have

$\begin{array}{c}\text{d}{U}_1=U_1\{-a_1\left(r\left(t\right)\right)+{\sigma }_1^2\left(r\left(t\right)\right)+{\displaystyle {\int }_{\mathbb{Y}}}\frac{{\gamma }_1^2\left(u\mathrm{,}r\left(t\right)\right)}{1+{\gamma }_1\left(u\mathrm{,}r\left(t\right)\right)}\lambda \left(\text{d}u\right)\\ +\frac{b_1\left(r\left(t\right)\right)}{U_1}+\frac{c_1\left(r\left(t\right)\right)U_1}{U_2+e\left(r\left(t\right)\right)U_1}\}\text{d}t-{\sigma }_1\left(r\left(t\right)\right)U_1\text{d}B\left(t\right)\\ -{\displaystyle {\int }_{\mathbb{Y}}}\frac{{\gamma }_1\left(u\mathrm{,}r\left(t\right)\right)}{1+{\gamma }_1\left(u\mathrm{,}r\left(t\right)\right)}{U}_1\tilde{N}\left(\text{d}t\mathrm{,}\text{d}u\right)\mathrm{.}\end{array}$

$\begin{array}{c}\text{d}{U}_2=U_2\{a_2\left(r\left(t\right)\right)+{\sigma }_1^2\left(r\left(t\right)\right)+{\displaystyle {\int }_{\mathbb{Y}}}\frac{{\gamma }_2^2\left(u\mathrm{,}r\left(t\right)\right)}{1+{\gamma }_2\left(u\mathrm{,}r\left(t\right)\right)}\lambda \left(\text{d}u\right)\\ +\frac{b_2\left(r\left(t\right)\right)}{U_2}-\frac{c_2\left(r\left(t\right)\right)U_2}{U_2+e\left(r\left(t\right)\right)U_1}\}\text{d}t-{\sigma }_2\left(r\left(t\right)\right)U_2\text{d}B\left(t\right)\\ -{\displaystyle {\int }_{\mathbb{Y}}}\frac{{\gamma }_2\left(u\mathrm{,}r\left(t\right)\right)}{1+{\gamma }_2\left(u\mathrm{,}r\left(t\right)\right)}{U}_2\tilde{N}\left(\text{d}t\mathrm{,}\text{d}u\right)\mathrm{.}\end{array}$

Noting that $G\left(\theta \right)$ is a nonsingular M-matrix, thank to Lemma 4.3, there exists $P={\left(p_1,p_2,\cdots ,p_S\right)}^{\text{T}}\gg 0$ such that $G\left(\theta \right)P\gg 0$, namely, for every $i\in \mathbb{S}$, ${\xi }_i\left(\theta \right)p_i-{\displaystyle \underset{j=1}{\overset{S}{\sum }}}{\gamma }_{ij}{p}_j>0$. So there exists $\kappa >0$ such that

${\xi }_i\left(\theta \right)p_i-{\displaystyle \underset{j=1}{\overset{S}{\sum }}}{\gamma }_{ij}{p}_j-\kappa p_i>0,\forall i\in \mathbb{S}.$

With the help of Itô’s formula, we have

$L\left[{\text{e}}^{\kappa t}{p}_i{\left(1+U\right)}^{\theta }\right]={\text{e}}^{\kappa t}\left\{\kappa p_i{\left(1+U\right)}^{\theta }+L\left[p_i{\left(1+U\right)}^{\theta }\right]\right\},$ (4.7)

where

$\begin{array}{l}\kappa p_i{\left(1+U\right)}^{\theta }+L\left[p_i{\left(1+U\right)}^{\theta }\right]\\ =\kappa p_i{\left(1+U\right)}^{\theta }+\theta p_i{\left(1+U\right)}^{\theta -1}{U}_1[-a_1\left(\left(r\left(t\right)\right)\right)+{\sigma }_1^2\left(\left(r\left(t\right)\right)\right)\begin{array}{c}\stackrel{}{}\\ \end{array}\\ +\frac{b_1\left(\left(r\left(t\right)\right)\right)}{U_1}+\frac{c_1\left(\left(r\left(t\right)\right)\right)U_1}{U_2+e\left(\left(r\left(t\right)\right)\right)U_1}]+\theta p_i{\left(1+U\right)}^{\theta -1}{U}_2[a_2\left(\left(r\left(t\right)\right)\right)\begin{array}{c}\stackrel{}{}\\ \end{array}\\ +{\sigma }_1^2\left(\left(r\left(t\right)\right)\right)+\frac{b_2\left(\left(r\left(t\right)\right)\right)}{U_2}-\frac{c_2\left(\left(r\left(t\right)\right)\right)U_2}{U_2+e\left(\left(r\left(t\right)\right)\right)U_1}]\end{array}$

$\begin{array}{l}+\theta p_i{\left(1+U\right)}^{\theta -1}U{\displaystyle {\int }_{\mathbb{Y}}}\frac{{\gamma }_1\left(u,r\left(t\right)\right)U_1+{\gamma }_2\left(u,r\left(t\right)\right)U_2}{U}\lambda \left(\text{d}u\right)\\ +\frac{\theta \left(\theta -1\right)}{2}{p}_i{\left(1+U\right)}^{\theta -2}{U}^2\frac{{\left[{\sigma }_1\left(\left(r\left(t\right)\right)\right)U_1+{\sigma }_2\left(\left(r\left(t\right)\right)\right)U_2\right]}^2}{U^2}\\ +{\displaystyle \underset{j=1}{\overset{S}{\sum }}}{\gamma }_{ij}{p}_j{\left(1+U\right)}^{\theta }+{\displaystyle {\int }_{\mathbb{Y}}}\left[p_i{\left(1+{\displaystyle \underset{j=1}{\overset{2}{\sum }}}\frac{U_j}{1+{\gamma }_1\left(u,r\left(t\right)\right)}\right)}^{\theta }-p_i{\left(1+U\right)}^{\theta }\right]\lambda \left(\text{d}u\right)\\ =\left[\mathcal{T}\left(U^{\theta }\right)\right]U^{\theta }+\mathcal{K}\left(U\right),\end{array}$ (4.8)

where

$\begin{array}{c}\mathcal{K}\left(U\right)=\left[{\left(1+U\right)}^{\theta }-U^{\theta }\right]\left\{\kappa p_i+{\displaystyle \underset{j=1}{\overset{S}{\sum }}}{\gamma }_{ij}{p}_j-{\displaystyle {\int }_{\mathbb{Y}}}{p}_i\lambda \left(\text{d}u\right)\right\}\\ +p_i\theta \left[{\left(1+U\right)}^{\theta -1}U-U^{\theta }\right]\{\frac{c_1\left(r\left(t\right)\right)U_1^2-c_2\left(r\left(t\right)\right)U_1^2}{U\left(U_2+e\left(r\left(t\right)\right)U_1\right)}\\ -\frac{a_1\left(r\left(t\right)\right)U_1-a_2\left(r\left(t\right)\right)U_2}{U}+\frac{{\sigma }_1^2\left(r\left(t\right)\right)U_1+{\sigma }_2^2\left(r\left(t\right)\right)U_2}{U}\end{array}$

$\begin{array}{c}+{\displaystyle {\int }_{\mathbb{Y}}}\frac{{\gamma }_1\left(u\mathrm{,}r\left(t\right)\right)U_1+{\gamma }_2\left(u\mathrm{,}r\left(t\right)\right)U_2}{U}\lambda \left(\text{d}u\right)\}\\ +\left[{\left(1+U\right)}^{\theta -2}{U}^2-U^{\theta }\right]\frac{\theta \left(\theta -1\right)}{2}{p}_i\left(\frac{{\sigma }_1\left(r\left(t\right)\right)U_1+{\sigma }_2\left(r\left(t\right)\right)U_2}{U}\right)\\ +p_i{\displaystyle {\int }_{\mathbb{Y}}}\left[{\left(1+{\displaystyle \underset{j=1}{\overset{2}{\sum }}}\frac{U_j}{\left(1+{\gamma }_j\left(u\mathrm{,}r\left(t\right)\right)\right)}\right)}^{\theta }-{\left({\displaystyle \underset{j=1}{\overset{2}{\sum }}}\frac{U_j}{\left(1+{\gamma }_j\left(u\mathrm{,}r\left(t\right)\right)\right)U}\right)}^{\theta }\right]\lambda \left(\text{d}u\right)\mathrm{.}\end{array}$

It is obvious that ${\lim }_{U\to +\infty }\frac{\mathcal{K}\left(U\right)}{U^{\theta }}=0$. Based on Jensen’s inequality, we derive

$\mathcal{T}\left[U^{\theta }\right]\le \kappa p_i+{\displaystyle \underset{j\mathrm{<mo>\; =\; </mo>1}}{\overset{S}{\sum }}}{\gamma }_{ij}{p}_j-p_i\theta B_1\left(r\left(t\right)\right)\frac{U_1}{U}-p_i\theta B_2\left(r\left(t\right)\right)\frac{U_2}{U}+p_i\frac{{\theta }^2{\sigma }^2}{2}$

$\begin{array}{c}+p_i\theta {\displaystyle {\int }_{\mathbb{Y}}}\left[\frac{U_1}{U}\ln \left(1+{\gamma }_1\left(u\mathrm{,}r\left(t\right)\right)\right)+\frac{U_2}{U}\ln \left(1+{\gamma }_2\left(u\mathrm{,}r\left(t\right)\right)\right)\right]\lambda \left(\text{d}u\right)\\ +p_i{\displaystyle {\int }_{\mathbb{Y}}}\left[{\left(\frac{U_1}{U\left(1+{\gamma }_1\left(u,r\left(t\right)\right)\right)}+\frac{U_2}{U\left(1+{\gamma }_2\left(u,r\left(t\right)\right)\right)}\right)}^{\theta }-1\right]\lambda \left(\text{d}u\right)\\ +p_i\theta \frac{c_1\left(i\right)}{e\left(r\left(t\right)\right)}+p_i\theta \frac{{\sigma }_1^2\left(r\left(t\right)\right)U_1+{\sigma }_2^2\left(r\left(t\right)\right)U_2}{2U}\\ -\frac{p_i\theta }{2}{\left(\frac{{\sigma }_1\left(r\left(t\right)\right)U_1+{\sigma }_2\left(r\left(t\right)\right)U_2}{U}\right)}^2\end{array}$

$\begin{array}{l}\le \kappa p_i+{\displaystyle \underset{j=1}{\overset{S}{\sum }}}{\gamma }_{ij}{p}_j-p_i\theta \underset{j=1,2}{\min }\left\{B_j\left(r\left(t\right)\right)\right\}+p_i{\theta }^2\frac{{\sigma }^2}{2}\\ +p_i{\displaystyle {\int }_{\mathbb{Y}}}\left[{\left(\frac{U_1}{\left(1+{\gamma }_1\right)U_3}+\frac{U_2}{\left(1+{\gamma }_2\right)U_3}\right)}^{\theta }-1+\theta \ln \left(1+\frac{{\gamma }_1{U}_1+{\gamma }_2{U}_2}{U_3}\right)\right]\lambda \left(\text{d}u\right)\\ +p_i\theta \frac{c_1\left(r\left(t\right)\right)}{e\left(r\left(t\right)\right)}+p_i\theta {\left[{\sigma }_1\left(r\left(t\right)\right)-{\sigma }_2\left(r\left(t\right)\right)\right]}^2\frac{U_1{U}_2}{2U^2}\end{array}$

$\begin{array}{l}\le \kappa p_i+{\displaystyle \underset{j=1}{\overset{S}{\sum }}}{\gamma }_{ij}{p}_j-p_i\theta B\left(i\right)+p_i{\theta }^2\frac{{\sigma }^2}{2}+p_i\theta \frac{c_1\left(r\left(t\right)\right)}{e\left(r\left(t\right)\right)}+p_i\theta \frac{{\left[{\sigma }_1\left(i\right)-{\sigma }_2\left(i\right)\right]}^2}{8}\\ +p_i{\displaystyle {\int }_{\mathbb{Y}}}\left[{\left(1+\underset{j=1,2}{\min }\left\{{\gamma }_{j*}\left(i\right)\right\}\right)}^{-\theta }-1+\theta \ln \left(1+\underset{j=1,2}{\min }\left\{{\gamma }_j^{*}\left(i\right)\right\}\right)\right]\lambda \left(\text{d}u\right)\\ :=\kappa p_i+{\displaystyle \underset{j=1}{\overset{S}{\sum }}}{\gamma }_{ij}{p}_j-p_i{\xi }_i\left(\theta \right),\end{array}$

that is

$\mathcal{T}\left[U^{\theta }\right]:=\kappa p_i+{\displaystyle \underset{j=1}{\overset{S}{\sum }}}{\gamma }_{ij}{p}_j-p_i{\xi }_i\left(\theta \right).$ (4.9)

In view of (4.7), (4.8) and (4.9), there exists $\mathcal{H}\left(\theta \right)>0$ such that

$\mathcal{L}\left[{\text{e}}^{\kappa t}{p}_i{\left(1+U\right)}^{\theta }\right]\le \mathcal{H}\left(\theta \right){\text{e}}^{\kappa t}.$ (4.10)

Integrating (4.10) from 0 to t and taking expectation of both sides, we have

$E\left[{\text{e}}^{\kappa t}{p}_i{\left(1+U\left(t\right)\right)}^{\theta }\right]\le \frac{\mathcal{H}\left(\theta \right)}{\kappa }\left({\text{e}}^{\kappa t}-1\right)+p_i{\left(1+U\left(0\right)\right)}^{\theta }\mathrm{.}$

This means that

$\begin{array}{l}\underset{t\to \infty }{\lim }E\left[{\left(1+U\left(t\right)\right)}^{\theta }\right]\le \frac{\mathcal{H}\left(\theta \right)}{\kappa \underset{i\in \mathbb{S}}{\min }{p}_i}\mathrm{.}\hfill \end{array}$ (4.11)

Base on (4.11), we obtain

$\underset{t\to \infty }{\lim \sup }E\left[{\left(U_1\left(t\right)+U_2\left(t\right)\right)}^{\theta }\right]\le \underset{t\to \infty }{\lim \sup }E\left[{\left(1+U\left(t\right)\right)}^{\theta }\right]\le \frac{\mathcal{H}\left(\theta \right)}{\kappa \underset{i\in \mathbb{S}}{\min }{p}_i}.$

From (4.6), we have

$\begin{array}{l}\underset{t\to \infty }{\lim \sup }E\left(\frac{1}{x^{\theta }\left(t\right)}\right)\le \frac{\mathcal{H}\left(\theta \right)}{\kappa \underset{i\in \mathbb{S}}{\min }{p}_i}=\mathrm{<mo>\; :\; </mo>}{H}_1\mathrm{,}\\ \underset{t\to \infty }{\lim \sup }E\left(\frac{1}{y^{\theta }\left(t\right)}\right)\le \frac{\mathcal{H}\left(\theta \right)}{\kappa \underset{i\in \mathbb{S}}{\min }{p}_i}=\mathrm{<mo>\; :\; </mo>}{H}_2\mathrm{,}\end{array}$ (4.12)

which is the required assertion.

Theorem 4.1. Under assumptions 2 and 4, if $\mathcal{B}>0$ , then system (1.6) is stochastically permanent.

Proof. By Lemma 4.4 and Lemma 4.5, we know

$\underset{t\to \infty }{\lim \sup }E\left[x{\left(t\right)}^{-\theta }\right]\le H_1\left(\theta \right)\mathrm{,}\underset{t\to \infty }{\lim \sup }E\left[y{\left(t\right)}^{-\theta }\right]\le H_2\left(\theta \right)\mathrm{.}$

Based on Chebyshev’s inequality, for any $ϵ\in \left(\mathrm{0,1}\right)$, there exists ${\delta }_{*}={\left(\frac{ϵ}{H\left(\theta \right)}\right)}^{\frac{1}{\theta }}>0$ such that

$p\left(\left|x\left(t\right)\right|<{\delta }_{*}\right)=p\left(\frac{1}{\left|x\left(t\right)\right|}>{\delta }_{*}\right)\le {\left({\delta }_{*}\right)}^{\theta }E\left(\frac{1}{x^{\theta }\left(t\right)}\right).$

So

$\underset{t\to \infty }{\lim \sup }P\left(\left|x\left(t\right)\right|<{\delta }_{*}\right)\le {\delta }_{*}^{\theta }H\left(\theta \right)<ϵ,$

Proof of $y\left(t\right)$ using the same method, namely,

$\underset{t\to \infty }{\lim \inf }P\left(\left|x\left(t\right)\right|\ge {\delta }_{\mathrm{<mo>\; *\; </mo>}}\right)\ge 1-ϵ\mathrm{,}\underset{t\to \infty }{\lim \inf }P\left(\left|y\left(t\right)\right|\ge {\delta }_{\mathrm{<mo>\; *\; </mo>}}\right)\ge 1-ϵ\mathrm{.}$

Finally, (4.2) is obtained by combining lemma 4.1 and Chebyshev inequality. Therefore, the system (1.6) is stochastically permanent.

In the same spirit as in the proof of Theorem 4.1, we yield the result of the subsystem (1.8) as follows.

Corollary 4.1. Under assumption 3, if $\beta \left(i\right)>0$ , then the subsystem (8) is stochastically permanent.

Proof. Note that for some $i\in \mathbb{S}$, $\beta \left(i\right)>0$. We can choose ${\theta }_0>0$ so small that for $0<\theta <{\theta }_0$,

$\begin{array}{l}{\xi }_i\left(\theta \right)=\theta \beta \left(i\right)-\frac{1}{2}{\theta }^2{\sigma }^2-\theta {\displaystyle {\int }_{\mathbb{Y}}}\ln \left(1+\underset{j=1,2}{\max }\left\{{\gamma }_j^{\mathrm{<mo>\; *\; </mo>}}\left(i\right)\right\}\right)\lambda \left(\text{d}u\right)\\ -{\displaystyle {\int }_{\mathbb{Y}}}\left[{\left(1+\underset{j=1,2}{\min }\left\{{\gamma }_{j\mathrm{<mo>\; *\; </mo>}}\left(i\right)\right\}\right)}^{-\theta }-1\right]\lambda \left(\text{d}u\right)>0.\end{array}$

So, $G\left(\theta \right)$ is a positive constant and is also a nonsingular M-matrix. According to Lemma 4.4 and Theorem 4.1, subsystem (1.8) is stochastically permanent.

5. Stationary Distribution

As far as we know, the existence of the ergodic stationary distribution of a stochastic competition model with high order stochastic perturbations has not been obtained theoretically. Therefore, this section mainly studies the existence and uniqueness of the stationary distribution of system (1.6) without Lévy jumps, namely system (5.1).

$\{\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),\\ \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).\end{array}$ (5.1)

Next, from the theorem in [30] we have the following lemma, we will use this lemma to prove ergodic stationary distribution.

Lemma 5.1. If the following conditions are satisfied:

1) ${\gamma }_{ij}>0$ for any $i\ne j$ ;

2) for each $k\in \mathbb{S}$, $D\left(x\mathrm{,}k\right)$ is symmetric and satisfies

$\lambda {\left|\zeta \right|}^2\le {\zeta }^{\text{T}}D\left(x\mathrm{,}k\right)\zeta \le {\lambda }^{-1}{\left|\zeta \right|}^2\text{for}\text{all}\zeta \in {ℝ}^n\mathrm{,}$

with some constant $\lambda \in \left(\mathrm{0,1}\right]$ for all $x\in {ℝ}^n$ ;

3) there exists a bounded open subset $\mathcal{D}$ of ${ℝ}^n$ with a regular boundary satisfying that, for each $k\in \mathbb{S}$ there exists a nonnegative function $V\left(\cdot \mathrm{,}k\right)\mathrm{<mo>\; :\; </mo>}{\mathcal{D}}^c\to ℝ$ such that $V\left(\cdot \mathrm{,}k\right)$ is twice continuously differential and that for some $\varsigma >0$,

$LV\left(x\mathrm{,}k\right)\le -\varsigma \mathrm{,}\left(x\mathrm{,}k\right)\in {\mathcal{D}}^c\times \mathbb{S}\mathrm{.}$

then $\left(x\left(t\right)\mathrm{,}y\left(t\right)\mathrm{,}r\left(t\right)\right)$ of system (5.1) is positive recurrent and ergodic. That is to say, there exsist a unique stationary distribution $\mu \left(\cdot \mathrm{,}\cdot \right)$ such that for any Borel measurable function $f\left(\cdot \mathrm{,}\cdot \right)\mathrm{<mo>\; :\; </mo>}{ℝ}^n\times \mathbb{S}\to ℝ$ satisfying

${\displaystyle \underset{k\in \mathbb{S}}{\sum }}{\displaystyle {\int }_{{ℝ}^n}}\left|f\left(x\mathrm{,}k\right)\right|\mu \left(x\mathrm{,}k\right)\text{d}x<\infty \mathrm{,}$

we have

$\mathbb{P}\left(\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t}f\left(x\left(s\right)\mathrm{,}r\left(s\right)\right)\text{d}s={\displaystyle \underset{k\in \mathbb{S}}{\sum }}{\displaystyle {\int }_{{ℝ}^n}}f\left(x\mathrm{,}k\right)\mu \left(x\mathrm{,}k\right)\text{d}x\right)=1.$

Theorem 5.1. Under $\pi {\varphi }_1>0$ , then for any given initial value $\left(\left(x\left(0\right)\mathrm{,}y\left(0\right)\right)\mathrm{,}r\left(0\right)\right)\in {ℝ}_{+}^2\times \mathbb{S}$ , system (5.1) has a unique ergodic stationary distribution $\mu (\cdot )$ on ${ℝ}_{+}^2\times \mathbb{S}$.

Proof. Define a $C^2$ -function $V\mathrm{<mo>\; :\; </mo>}{ℝ}_{+}^2\times \mathbb{S}\to {ℝ}_{+}$ by

$\begin{array}{l}V\left(x\mathrm{,}y\mathrm{,}i\right)=V_1\left(x\mathrm{,}y\mathrm{,}i\right)+V_2\left(x\mathrm{,}y\mathrm{,}i\right)\mathrm{,}\hfill \end{array}$

where

$V_1\left(x,y,i\right)=x^{\vartheta }+y^{\vartheta },0<\vartheta <1;V_2\left(x,y,i\right)=\frac{1}{x^{\theta }}+\frac{1}{y^{\theta }}.$

Then

$\begin{array}{l}L{V}_1\left(x,y,i\right)\\ =\vartheta x^{\vartheta }\left[a_1\left(i\right)-b_1\left(i\right)x-\frac{c_1\left(i\right)y}{x+e\left(i\right)y}+\frac{1}{2}\left(\vartheta -1\right){\sigma }_1^2\left(i\right)\right]\end{array}$

$\begin{array}{l}+\vartheta y^{\vartheta }\left[-a_2\left(i\right)-b_2\left(i\right)y+\frac{c_2\left(i\right)x}{x+e\left(i\right)y}+\frac{1}{2}\left(\vartheta -1\right){\sigma }_2^2\left(i\right)\right]\\ +{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}{x}^{\vartheta }+{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}{y}^{\vartheta }\end{array}$

$\begin{array}{l}=-b_1\left(i\right)\vartheta x^{\vartheta +1}+\vartheta x^{\vartheta }\left[a_1\left(i\right)-\frac{c_1\left(i\right)y}{x+e\left(i\right)y}+\frac{1}{2}\left(\vartheta -1\right){\sigma }_1^2\left(i\right)\right]\\ -b_2\left(i\right)\vartheta y^{\vartheta +1}+\vartheta y^{\vartheta }\left[-a_2\left(i\right)+\frac{c_2\left(i\right)x}{x+e\left(i\right)y}+\frac{1}{2}\left(\vartheta -1\right){\sigma }_2^2\left(i\right)\right]\\ +{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}{x}^{\vartheta }+{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}{y}^{\vartheta }\\ \le -b_1\left(i\right)\vartheta x^{\vartheta +1}+x^{\vartheta }\left[\vartheta a_1\left(i\right)+\frac{\vartheta }{2}\left(\vartheta -1\right){\sigma }_1^2\left(i\right)+{\displaystyle \underset{j\ne i}{\sum }}{\gamma }_{ij}\right]\\ -b_2\left(i\right)\vartheta y^{\vartheta +1}+y^{\vartheta }\left[-\vartheta a_2\left(i\right)+\vartheta c_2\left(i\right)-\frac{\vartheta }{2}\left(\vartheta -1\right){\sigma }_2^2\left(i\right)+{\displaystyle \underset{j\ne i}{\sum }}{\gamma }_{ij}\right],\end{array}$ (5.2)

and

$\begin{array}{l}L{V}_2\left(x,y,i\right)\\ =-\theta x^{-\theta }\left(a_1\left(i\right)-\frac{c_1\left(i\right)y}{x+e\left(i\right)y}-\frac{\left(\theta +1\right)}{2}{\sigma }_1^2\left(i\right)\right)+x^{1-\theta }\theta b_1\left(i\right)+{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}{x}^{-\vartheta }\\ -\theta y^{-\theta }\left(-a_2\left(i\right)+\frac{c_2\left(i\right)x}{x+e\left(i\right)y}-\frac{\left(\theta +1\right)}{2}{\sigma }_2^2\left(i\right)\right)+y^{1-\theta }\theta b_2\left(i\right)+{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}{y}^{-\vartheta }\\ =\theta x^{-\theta }\left[-\left(a_1\left(i\right)-\frac{c_1\left(i\right)y}{x+e\left(i\right)y}-b_1\left(i\right)x\right)+\frac{\left(\theta +1\right)}{2}{\sigma }_1^2\left(i\right)+\frac{1}{\theta }{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}\right]\\ +\theta y^{-\theta }\left[-\left(-a_2\left(i\right)+\frac{c_2\left(i\right)x}{x+e\left(i\right)y}-b_2\left(i\right)y\right)+\frac{\left(\theta +1\right)}{2}{\sigma }_2^2\left(i\right)+\frac{1}{\theta }{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}\right]\end{array}$

$\begin{array}{l}\le \theta x^{-\theta }\left\{-\left[{\varphi }_1\left(i\right)-\frac{\theta }{2}{\sigma }_1^2\left(i\right)\right]+\frac{c_1\left(i\right)}{e\left(i\right)}+b_1\left(i\right)x+\frac{1}{\theta }{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}\right\}\\ +\theta y^{-\theta }\left[-\left({\varphi }_2\left(i\right)-\frac{\theta }{2}{\sigma }_2^2\left(i\right)\right)+b_2\left(i\right)y+\frac{1}{\theta }{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}\right]\\ =\theta x^{-\theta }\left\{-\left[\pi {\varphi }_1-\frac{\theta }{2}{\sigma }_1^2\left(i\right)+\frac{1}{\theta }{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}\right]+\frac{c_1\left(i\right)}{e\left(i\right)}\right\}+\theta b_1\left(i\right)x^{1-\theta }\\ +\theta y^{-\theta }\left[-\left(\pi {\varphi }_2-\frac{\theta }{2}{\sigma }_2^2\left(i\right)+\frac{1}{\theta }{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}\right)\right]+\theta b_2\left(i\right)y^{1-\theta },\end{array}$ (5.3)

where

${\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).$

Choose a constant $0<\tilde{\theta }\le {\theta }_0$ such that for any $0<\theta \le \tilde{\theta }$,

${\varpi }_l:=\pi {\varphi }_l-\frac{\theta }{2}{\sigma }_l^2\left(i\right)+\frac{1}{\theta }{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}>0,i\in \mathbb{S},l=1,2.$

Then

$L{V}_2\left(x\mathrm{,}y\mathrm{,}i\right)\le -\theta x^{-\theta }\left({\varpi }_1+\frac{c_1\left(i\right)}{e\left(i\right)}\right)+\theta b_1\left(i\right)x^{1-\theta }-\theta y^{-\theta }{\varpi }_2+\theta b_2\left(i\right)y^{1-\theta }\mathrm{.}$

Combining (5.2) and (5.3), we obtain

$\begin{array}{c}LV\left(x,y,i\right)\le -b_1\left(i\right)\vartheta x^{\vartheta +1}+x^{\vartheta }\left[\vartheta a_1\left(i\right)+\frac{\vartheta }{2}\left(\vartheta -1\right){\sigma }_1^2\left(i\right)+{\displaystyle \underset{j\ne i}{\sum }}{\gamma }_{ij}\right]\\ -\theta x^{-\theta }\left({\varpi }_1+\frac{c_1\left(i\right)}{e\left(i\right)}\right)+\theta b_1{x}^{1-\theta }-b_2\left(i\right)\vartheta y^{\vartheta +1}\\ +y^{\vartheta }\left[-\vartheta a_2\left(i\right)+\vartheta c_2\left(i\right)-\frac{\vartheta }{2}\left(\vartheta -1\right){\sigma }_2^2\left(i\right)+{\displaystyle \underset{j\ne i}{\sum }}{\gamma }_{ij}\right]\\ -\theta y^{-\theta }{\varpi }_2+\theta b_2\left(i\right)y^{1-\theta }\\ :=h\left(x\right)+g\left(y\right),\end{array}$

where

$\begin{array}{c}h\left(x\right)=-b_1\left(i\right)\vartheta x^{\vartheta +1}+\vartheta x^{\vartheta }\left[a_1\left(i\right)+\frac{1}{2}\left(\vartheta -1\right){\sigma }_1^2\left(i\right)+{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}\right]\\ -\theta x^{-\theta }\left({\varpi }_1+\frac{c_1\left(i\right)}{e\left(i\right)}\right)+\theta b_1\left(i\right)x^{1-\theta },\end{array}$

$\begin{array}{c}g\left(y\right)=-b_2\left(i\right)\vartheta y^{\vartheta +1}+\vartheta y^{\vartheta }\left[-a_2\left(i\right)+c_2\left(i\right)-\frac{1}{2}\left(\vartheta -1\right){\sigma }_2^2\left(i\right)+{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\gamma }_{ij}\right]\\ -\theta y^{-\theta }{\varpi }_2+\theta b_2\left(i\right)y^{1-\theta }.\end{array}$

Then we have

$LV\left(x\mathrm{,}y\mathrm{,}i\right)\le \{\begin{array}{l}h\left(x\right)+g^{\mathcal{L}}\left(y\right)\to -\infty \mathrm{,}\text{as}x\to 0^{+}\mathrm{,}\hfill \\ h^{\mathcal{L}}\left(x\right)+g\left(y\right)\to -\infty \mathrm{,}\text{as}y\to 0^{+}\mathrm{,}\hfill \\ h\left(x\right)+g^{\mathcal{L}}\left(y\right)\to -\infty \mathrm{,}\text{as}x\to +\infty \mathrm{,}\hfill \\ h^{\mathcal{L}}\left(x\right)+g\left(y\right)\to -\infty \mathrm{,}\text{as}y\to +\infty \mathrm{,}\hfill \end{array}$

where

$h^{\mathcal{L}}={\sup }_{t\in {ℝ}_{+}}h\left(t\right),g^{\mathcal{L}}={\sup }_{t\in {ℝ}_{+}}g\left(t\right).$

To sum up, when $x\to 0^{+}$ or $y\to 0^{+}$, $x\to 0^{+}$ or $x\to +\infty$ or $y\to +\infty$, we can derive $LV\to -\infty$. Consequently, take $k>0$ sufficiently small and let $K:=\left[\varsigma ,\frac{1}{\varsigma }\right]\times \left[\varsigma ,\frac{1}{\varsigma }\right]$, one can see that

$LV\left(x\mathrm{,}y\mathrm{,}i\right)\le -\mathrm{1,}\left(x\mathrm{,}y\mathrm{,}i\right)\in K^c\times \mathbb{S}\mathrm{.}$

That is to say, the condition (3) in Lemma 5.1 is satisfied. On the other hand, Assumption 3 indicates the condition (1) that satisfies Lemma 5.1, and the diffusion matrix $D\left(x\mathrm{,}y\mathrm{,}i\right)=\text{diag}\left\{{\sigma }_1^2\left(i\right)x^2\mathrm{,}{\sigma }_2^2\left(i\right)y^2\right\}$ of system (5.1) is positive definite, which implies that condition (2) in Lemma 5.1 holds. According to lemma 5.1, system (5.1) is ergodic and has a unique stationary distribution. The proof is completed.

6. Conclusions and Example

The stochastic persistence and extinction of a stochastic ratio-dependent predator-prey system with Markovian switching and Lévy noise are studied. Our main results are as follows: Theorem 3.1 gives sufficient conditions on extinction and non-persistent in mean of each population. Theorem 3.3 gives sufficient conditions for each population to be weakly persistent. Further, Corollaries 3.1 tells us that for some $i\in S$, if $B_1\left(i\right)>c_1\left(i\right)$, system (1.8) is weakly persistent, and if $B_1\left(i\right)<0,B_2\left(i\right)<0$, system (1.8) is extinct. Theorem 4.1 gives sufficient condition for the stochastic permanence of (1.6). Theorem 4.1 and Corollaries 4.1 tell us that if some subsystems are stochastically persistent and others are extinctive, then as the result of the Markovian switch, the system (1.6) is still stochastically persistent. Finally, when there is no Lévy jump, a sufficient condition for the existence of ergodic stationary distribution of system (5.1) is established under certain conditions.

Next, two examples are introduced to support our theoretical analysis results.

Example 6.1. First of all, we will discuss the effect of Markov switches on population dynamics. Consider (1.8) with the Markov chain $r\left(t\right)$ taking value in the state space value $\mathbb{S}=\left\{1,2\right\}$, be deemed to have been the result of the switching between

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

and

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

where $\lambda \left(\mathbb{Y}\right)=1$, the initial value $x\left(0\right)=4,y\left(0\right)=2$ and the coefficients

$\{\begin{array}{l}{a}_1\left(1\right)=2,b_1\left(1\right)=1,c_1\left(1\right)=4,e\left(1\right)=2,{\sigma }_1\left(1\right)=2.5,{\gamma }_1\left(1,u\right)=1,\hfill \\ a_2\left(1\right)=3,b_2\left(1\right)=2,c_2\left(1\right)=4,e\left(1\right)=3,{\sigma }_2\left(1\right)=2,{\gamma }_2\left(1,u\right)=1,\hfill \\ a_1\left(2\right)=5,b_1\left(2\right)=1,c_1\left(2\right)=4,e\left(2\right)=2,{\sigma }_1\left(2\right)=1,{\gamma }_1\left(2,u\right)=1,\hfill \\ a_2\left(2\right)=0.5,b_2\left(2\right)=3,c_2\left(2\right)=5,e\left(2\right)=3,{\sigma }_2\left(2\right)=1,{\gamma }_2\left(2,u\right)=1.\hfill \end{array}$ (6.3)

Case 1. Let the generator of Markov chain $r\left(t\right)$ be

$\Gamma =\left(\begin{array}{cc}-4& 4\\ 1& -1\end{array}\right).$

Solving (2.3) yields the unique stationary distribution

$\phi =\left({\phi }_1,{\phi }_2\right)=\left(\frac{1}{5},\frac{4}{5}\right).$

Based on (6.3), we compute

$\{\begin{array}{l}{B}_1\left(1\right)=-1.4319<0,B_2\left(1\right)=-1.3069<0,B\left(1\right)=-1.4319<0,\hfill \\ B_1\left(2\right)=4.1931>0,B_2\left(2\right)=3.6931>0,B\left(2\right)=3.6931>0,\hfill \\ \mathcal{B}=2.6681>0.\hfill \end{array}$ (6.4)

From (6.4) and Corollary 4.2, it is easy to find that the subsystem (6.1) is extinct, but subsystem (6.2) is persistent. According to Theorem 4.1, system (1.6) is stochastically permanent. This example shows that although some subsystems are impermanent, the overall behavior of system (1.6) is stochastically permanent as a result of Markovian switching. Thus, Markovian switching may contribute to permanence to some extent. The numerical result is shown in Figure 1(a), Figure 1(b) and Figure 1(c).

Case2. Let $\Gamma =\left(\begin{array}{cc}-1& 1\\ 3& -3\end{array}\right)$ and $\phi =\left({\phi }_1,{\phi }_2\right)=\left(\frac{3}{4},\frac{1}{4}\right)$. We obtain

$\{\begin{array}{l}{\mathcal{B}}_1=-0.02565<0,\\ {\mathcal{B}}_2=-0.0569<0.\end{array}$

According to Theorem 4.1, system (1.6) is extinctive (see Figure 1(d)).

Example 6.2. Next, we explain the impacts of Lévy jumps on population dynamics. Consider (1.6) with Markov chain $r\left(t\right)$ taking value in state space $\mathbb{S}=\left\{1,2\right\}$. Let $\Gamma =\left(\begin{array}{cc}-1& 1\\ 2& -2\end{array}\right)$ and $\phi =\left({\phi }_1,{\phi }_2\right)=\left(\frac{2}{3},\frac{1}{3}\right)$, $\lambda \left(\mathbb{Y}\right)=1$, and the coefficients

$\{\begin{array}{l}{a}_1\left(1\right)=0.5,b_1\left(1\right)=1,c_1\left(1\right)=1,e\left(1\right)=2,{\sigma }_1\left(1\right)=0.7,\hfill \\ a_1\left(2\right)=0.3,b_1\left(2\right)=1,c_1\left(2\right)=1,e\left(2\right)=2,{\sigma }_1\left(2\right)=0.7,\hfill \\ a_2\left(1\right)=0.1,b_2\left(1\right)=2,c_2\left(1\right)=2,e\left(1\right)=3,{\sigma }_2\left(1\right)=0.5,\hfill \\ a_2\left(2\right)=0.25,b_2\left(2\right)=3,c_2\left(2\right)=2,e\left(2\right)=3,{\sigma }_2\left(2\right)=\mathrm{0.5.}\hfill \end{array}$

Case 1. Under the condition of without Lévy jumps effect, a simple computation yields

$\{\begin{array}{l}{B}_1\left(1\right)=0.2550,B_2\left(1\right)=1.7750,\hfill \\ B_1\left(2\right)=0.0550,B_2\left(2\right)=1.6250,\hfill \\ {\mathcal{B}}_1=0.188333>0,\hfill \\ {\mathcal{B}}_2=1.71167>0.\hfill \end{array}$

Figure 1. Time series of system (6) for ${\gamma }_1\left(1,u\right)={\gamma }_2\left(1,u\right)={\gamma }_1\left(2,u\right)={\gamma }_2\left(2,u\right)=1$, $\mathbb{Y}=\left(0,\infty \right)$, $\lambda \left(\mathbb{Y}\right)=1$, step size $\Delta t=0.02$, the initial data $x\left(0\right)=4$, $y\left(0\right)=2$. (a) is with $a_1\left(1\right)=2$, $b_1\left(1\right)=1$, $c_1\left(1\right)=c_2\left(1\right)=4$, $e\left(1\right)=2$, ${\sigma }_1\left(1\right)=2.5$, $a_2\left(1\right)=3$, $b_2\left(1\right)=2$, $e\left(1\right)=3$. (b) is with $a_1\left(2\right)=5$, $b_1\left(2\right)=1$, $c_1\left(2\right)=4$, $e\left(2\right)=2$, $a_2\left(2\right)=0.5$, $b_2\left(2\right)=3$, $c_2\left(2\right)=5$, $e\left(2\right)=3$, ${\sigma }_2\left(2\right)={\sigma }_1\left(2\right)=1$. ((c), (d)) The parameters are the same as in (a) and (b). These two sets of figures suggest that populations may switch between extinction and persistence as a result of the Markovian switch.

Therefore, by Theoerm 4.1, system (6) is stochastically permanent (see Figure 2(a)).

Case2. When ${\gamma }_1\left(1,u\right)=-0.4$, ${\gamma }_1\left(2,u\right)=-0.6$, ${\gamma }_2\left(1,u\right)=-0.5$, ${\gamma }_2\left(2,u\right)=-0.3$, that is, there is jump noise. By calculation

$\{\begin{array}{l}{B}_1\left(1\right)=-1.0613,B_2\left(1\right)=0.5819,\hfill \\ B_1\left(2\right)=-1.1381,B_2\left(2\right)=-0.8608,\hfill \\ {\mathcal{B}}_1=-0.9243<0,{\mathcal{B}}_2=-0.0300667<0,\hfill \\ \phi B_1=-0.243233<0,\phi B_2=-1.04567<0.\hfill \end{array}$

Figure 2. Time series of system (1.6) for $a_1\left(1\right)=0.5$, $a_1\left(2\right)=0.3$, $b_1\left(1\right)=b_1\left(2\right)=1$, $c_1\left(1\right)=c_1\left(2\right)=1$, $e\left(1\right)=e\left(2\right)=2$, ${\sigma }_1\left(1\right)={\sigma }_1\left(2\right)=0.7$, $a_2\left(1\right)=0.1$, $a_2\left(2\right)=0.25$, $b_2\left(2\right)=3$, $b_2\left(1\right)=2$, $c_2\left(1\right)=c_2\left(2\right)=2$, $e\left(1\right)=e\left(2\right)=3$, ${\sigma }_2\left(1\right)={\sigma }_2\left(2\right)=0.5$. step size $\Delta t=0.02$, the initial data $x\left(0\right)=4$, $y\left(0\right)=2$. (a) is with ${\gamma }_1\left(1,u\right)={\gamma }_1\left(2,u\right)={\gamma }_2\left(1,u\right)={\gamma }_2\left(2,u\right)=0$. (b) is with ${\gamma }_1\left(1,u\right)=-0.4$, ${\gamma }_1\left(2,u\right)=-0.6$, ${\gamma }_2\left(1,u\right)=-0.5$, ${\gamma }_2\left(2,u\right)=-0.3$. These two sets of figures show that Lévy jumps suppression population of persistent.

By Theoerm 3.1, system (1.6) is extinctive. This example suggests that Lévy jumps may suppress the permanence (see Figure 2(b)). The results suggest that Lévy jumps may suppress for the persistence of species.

Acknowledgements

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

Conflicts of Interest

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

References

[1]	Jan Maiti, A.M.M. and Samanta, G.P. (2013) Deterministic and Stochastic Analysis of a Ratio-Dependent Predator-Prey System with Delay. Nonlinear Analysis: Modelling and Control, 3, 383-398. https://doi.org/10.15388/NA.2007.12.3.14700
[2]	Arditi, R. and Ginzburg, L.R. (1989) Coupling in Predator-Prey Dynamics: Ratio-Dependence. Journal of Theoretical Biology, 139, 311-326. https://doi.org/10.1016/S0022-5193(89)80211-5
[3]	Arditi, R. and Saiah, H. (1992) Empirical Evidence of the Role of Heterogeneity in Ratio-Dependent Consumption. Ecology, 73, 1544-1551. https://doi.org/10.2307/1940007
[4]	Ji, C. and Jiang, D. (2009) Analysis of a Predator-Prey Model with Modified Leslie-Gower and Holling-Type II Schemes with Stochastic Perturbation. Journal of Mathematical Analysis and Applications, 359, 482-498. https://doi.org/10.1016/j.jmaa.2009.05.039
[5]	Holling, C.S. (1965) The Functional Response of Predator to Prey Density and Its Role in Mimicry and Population Regulation. Memoirs of the Entomological Society of Canada, 97, 5-60. https://doi.org/10.4039/entm9745fv
[6]	Liu, M. and Wang, K. (2013) Analysis of a Stochastic Autonomous Mutualism Model. Journal of Mathematical Analysis and Applications, 402, 392-403. https://doi.org/10.1016/j.jmaa.2012.11.043
[7]	Nguyen, T.H.L. and Ta, V.T. (2011) Dynamics of a Stochastic Ratio-Dependent Predator-Prey Model. Analysis and Applications, 9, 329-344. https://doi.org/10.1142/S0219530511001868
[8]	Tapaswi, P.K. and Mukhopadhyay, A. (1999) Effects of Environmental Fluctuation on Plankton Allelopathy. Journal of Mathematical Biology, 39, 39-58. https://doi.org/10.1007/s002850050162
[9]	Liu, M., Deng, M. and Du, B. (2015) Analysis of a Stochastic Logistic Model with Diffusion. Applied Mathematics and Computation, 266, 169-182. https://doi.org/10.1016/j.amc.2015.05.050
[10]	Zhang, X., Li, W., Liu, M. and Wang, K. (2015) Dynamics of a Stochastic Holling II One-Predator Two-Prey System with Jumps. Physica A, 421, 571-582. https://doi.org/10.1016/j.physa.2014.11.060
[11]	Zou, X. and Wang, K. (2014) Optimal Harvesting for a Stochastic Regime-Switching Logistic Diffusion System with Jumps. Nonlinear Analysis: Hybrid Systems, 13, 32-44. https://doi.org/10.1016/j.nahs.2014.01.001
[12]	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
[13]	Liu, M. and Wang, K. (2013) Dynamics of a Leslie-Gower Holling-Type II Predator-Prey System with Lévy Jumps. Nonlinear Analysis, 85, 204-213. https://doi.org/10.1016/j.na.2013.02.018
[14]	Zhu, C. and Yin, G. (2009) On Hybrid Competitive Lotka-Volterra Ecosystems. Nonlinear Analysis, 71, e1370-e1379. https://doi.org/10.1016/j.na.2009.01.166
[15]	Wang, S., Hu, G., Wei, T. and Wang, L. (2020) Permanence of Hybrid Competitive Lotka-Volterra System with Lévy Noise. Physica A: Statistical Mechanics and Its Applications, 495, Article ID: 123116. https://doi.org/10.1016/j.physa.2019.123116
[16]	Li, X., Alison, A., Jiang, D. and Mao, X. (2011) Sufficient and Necessary Conditions of Stochastic Permanence and Extinction for Stochastic Logistic Populations under Regime Switching. Journal of Mathematical Analysis and Applications, 376, 11-28. https://doi.org/10.1016/j.jmaa.2010.10.053
[17]	Li, X., Jiang, D. and Mao, X. (2009) Population Dynamical Behavior of Lotka-Volterra System under Regime Switching. Journal of Computational and Applied Mathematics, 232, 427-448. https://doi.org/10.1016/j.cam.2009.06.021
[18]	Mao, X. and Yuan, C. (2006) Stochastic Differential Equations with Markovian Switching. Imperial College Press, London. https://doi.org/10.1142/p473
[19]	Takeuchi, Y., Du, N.H., Hieu, N.T. and Sato, K. (2006) Evolution of Predator-Prey Systems Described by a Lotka-Volterra Equation under Random Environment. Journal of Mathematical Analysis and Applications, 323, 938-957. https://doi.org/10.1016/j.jmaa.2005.11.009
[20]	Zhao, X. and Zeng, Z. (2020) Stationary Distribution and Extinction of a Stochastic Ratio-Dependent Predator-Prey System with Stage Structure for the Predator. Physica A, 545, Article ID: 123310. https://doi.org/10.1016/j.physa.2019.123310
[21]	Ji, C., et al. (2011) Qualitative Analysis of a Stochastic Ratio-Dependent Predator-Prey System. Journal of Computational and Applied Mathematics, 235, 1326-1341. https://doi.org/10.1016/j.cam.2010.08.021
[22]	Ouyang, M. and Li, X. (2015) Permanence and Asymptotical Behavior of Stochastic Prey-Predator System with Markovian Switching. Applied Mathematics and Computation, 266, 539-559. https://doi.org/10.1016/j.amc.2015.05.083
[23]	Bai, L., et al. (2014) Analysis of a Stochastic Ratio-Dependent Predator-Prey Model Driven by Lévy Noise. Applied Mathematics and Computation, 233, 480-493. https://doi.org/10.1016/j.amc.2013.12.187
[24]	Mandal, P. and Banerjee, M. (2012) Stochastic Persistence and Stationary Distribution in a Holling-Tanner Type Prey-Predator Model. Physica A, 391, 1216-1233. https://doi.org/10.1016/j.physa.2011.10.019
[25]	Mao, X. (1997) Stochastic Differential Equations and Applications. Horwood Publishing, Chichester.
[26]	Samanda, G.P. (1996) Influence of Environmental Noises on the Gomatam Model of Interacting Species. Ecological Modelling, 91, 283-291. https://doi.org/10.1016/0304-3800(95)00195-6
[27]	Li, X., Mao, X. and Shen, Y. (2010) Approximation Solutions of Stochastic Differential Delay Equations with Regime Switching. Journal of Difference Equations and Applications, 19, 195-207. https://doi.org/10.1080/10236190802695456
[28]	Gao, H. and Wang, Y. (2019) Stochastic Mutualism Model under Regime Switching with Lévy Jumps. Physica A, 515, 355-375. https://doi.org/10.1016/j.physa.2018.09.189
[29]	Lipster, R. (1980) A Strong Law of Large Numbers for Local Martingales. Stochastics, 3, 217-228. https://doi.org/10.1080/17442508008833146
[30]	Yin, G. and Zhu, C. (2007) Asymptotic Properties of Hybrid Diffusion Systems. SIAM Journal on Control and Optimization, 46, 1155-1179. https://doi.org/10.1137/060649343