Dynamics of a stochastic modified Leslie–Gower predator–prey system with hunting cooperation

Abstract

In this paper, we consider a stochastic two-species predator–prey system with modified Leslie–Gower. Meanwhile, we assume that hunting cooperation occurs in the predators. By using Itô formula and constructing a proper Lyapunov function, we first show that there is a unique global positive solution for any given positive initial value. Furthermore, based on Chebyshev inequality, the stochastic ultimate boundedness and stochastic permanence are discussed. Then, under some conditions, we prove the persistence in mean and extinction of system. Finally, we verify our results by numerical simulations.

KEYWORDS:

• Stochastic predator–prey system
• hunting cooperation
• stochastic ultimate boundedness
• stochastic permanence
• extinction

2020 MATHEMATICS SUBJECT CLASSIFICATIONS:

• 60H10
• 92D25
• 34F05

1. Introduction

It is well known that the predator–prey relationship widely exists between two species. Therefore, predator–prey systems have been studied extensively, and the rich results have been achieved. Leslie  [1] and Leslie et al. [2] assumed that the carrying capacity of the predator's environment is proportional to the number of prey and introduced the following predator–prey system: (1) $\{\begin{array}{l}\dot{x}\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t\left)\right]-p\left(x\right(t\left)\right)y\left(t\right),\\ {\displaystyle \dot{y}\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{x\left(t\right)}).}\end{array}$(1) Here $x\left(t\right)$ and $y\left(t\right)$ denote the population densities of prey and predator at time t. $a,b,r$, and f are all positive parameters. a and r are the intrinsic growth rate of the prey and predator, respectively. b means the competitive strength among individuals of the prey. f is the maximum value of the per capita reduction rate of predator. The term $p\left(x\right(t\left)\right)$ is called the functional response function, which describes the number of prey successfully attacked by individual predators. The term $\frac{\mathrm{fy}\left(t\right)}{x\left(t\right)}$ is called the Leslie–Gower term, and it implies that the carrying capacity of the predator's environment is proportional to the number of prey, and not a constant (Logistic form).

Aziz-Alaoui [3] and Nindjin et al. [4] explained that the predator may convert to catch other populations when the prey x is severely scarce. However, its growth will be limited even if its most favourite food x is not available in abundance. In order to solve such a problem, they thought it was necessary to add a positive constant m to the denominator. Thus, system (1(1) $\{\begin{array}{l}\dot{x}\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t\left)\right]-p\left(x\right(t\left)\right)y\left(t\right),\\ {\displaystyle \dot{y}\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{x\left(t\right)}).}\end{array}$(1) ) becomes the following predator–prey system with modified Leslie–Gower: (2) $\{\begin{array}{l}\dot{x}\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t\left)\right]-p\left(x\right(t\left)\right)y\left(t\right),\\ {\displaystyle \dot{y}\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)}).}\end{array}$(2) In general, the functional response function $p\left(x\right(t\left)\right)$ represents the attack power of predator on prey. In the classical Lotka–Volterra system, $p\left(x\right(t\left)\right)$ takes the form $p\left(x\right(t\left)\right)=\mathrm{px}\left(t\right)$. This is reasonable since the number of prey hunted by the individual predator is naturally proportional to the number of the prey. Further, this function was revised by Holling, and three functional response functions were proposed to describe the predator behaviour on prey. One of Holling functional response functions is Holling II functional response, which takes the following form $p\left(x\right(t\left)\right)=\frac{\mathrm{\rho x}\left(t\right)}{1+\mathrm{h\rho x}\left(t\right)}.$ρ is the encounter rate between predators and their prey, and h is the handling time of predator on prey. In recent years, the Leslie–Gower predator–prey model with Holling II functional response has been widely studied by many scholars. Arancibia-Ibarra et al. [5] considered a modified Leslie–Gower system with Holling-type II schemes and Allee effect. It was revealed that system with a strong Allee effect had two positive equilibrium points and system with a weak Allee effect had three positive equilibrium points. Abid et al. [6] investigated a predator–prey model with modified Leslie–Gower and Holling-type II functional response. They discussed the existence and boundedness of solutions, then analysed the existence and local stability of the possible steady states. Moreover, the conditions of global stability of the interior equilibrium were established.

In nature, cooperative behaviour within populations is a very widespread phenomenon. Alves et al. [7] advised that cooperative predators may benefit from their behaviour since the success of predator's attacks on prey increases. They added a cooperation term to describe the attack of a predator on a prey and studied a predator–prey system with hunting cooperation. The attack rate (or encounter rate) of the predator was assumed to be $\rho =\rho \left(y\right)=\lambda +\mathrm{\alpha y}$. Therefore, the corresponding Holling II functional response function with hunting cooperation is as follows: $p\left(x\right(t),y(t\left)\right)=\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)},$where $\lambda (>0)$ is the attack rate of the per predator on the prey, $\alpha (>0)$ represents the intensity of the predator cooperation in hunting. Then, based on system (2(2) $\{\begin{array}{l}\dot{x}\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t\left)\right]-p\left(x\right(t\left)\right)y\left(t\right),\\ {\displaystyle \dot{y}\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)}).}\end{array}$(2) ), we have the following system (3) $\{\begin{array}{l}{\displaystyle \dot{x}\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t\left)\right]-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)},}\\ {\displaystyle \dot{y}\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)}).}\end{array}$(3) On the other hand, in the real world, environmental noise exists everywhere, which is bound to have an impact on the development of populations. Therefore, many scholars introduced white noise into predator–prey systems and studied the influence of white noise on predator–prey systems. Xu et al. [8] considered a stochastic predator–prey system with modified Leslie–Gower and Holling-type IV schemes. They proved the stochastic boundedness and stochastic permanence, then they investigated the persistence in mean and extinction. Meanwhile, the existence of a ergodic stationary distribution was discussed. Zhou et al. [9] studied a stochastic predator–prey model with modified Leslie–Gower and Holling-type II schemes. Especially, sufficient conditions of persistence in mean and extinction were given. The rich results for the stochastic population systems can also see [10–12].

In this paper, we assume that in system (3(3) $\{\begin{array}{l}{\displaystyle \dot{x}\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t\left)\right]-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)},}\\ {\displaystyle \dot{y}\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)}).}\end{array}$(3) ), the intrinsic growth rates a and r of the prey and predator are affected by white noise, that is, $a\to a+{\sigma }_1{\dot{B}}_1\left(t\right),r\to r+{\sigma }_2{\dot{B}}_2\left(t\right).$where $B_1\left(t\right),B_2\left(t\right)$ are mutually independent Brownian motions, and ${\sigma }_1^2,{\sigma }_2^2$ represent the intensities of the white noises. Then, corresponding to system (3(3) $\{\begin{array}{l}{\displaystyle \dot{x}\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t\left)\right]-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)},}\\ {\displaystyle \dot{y}\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)}).}\end{array}$(3) ), a stochastic two-species predator–prey system with modified Leslie–Gower and hunting cooperation among predators can be described as follows: (4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) The dynamic behaviours of the stochastic predator–prey system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) may provide effective ideas for population protection. From the biological sense, the first concern is that the solution of the system is non-negative, then the stochastic permanence and extinction are also very significant issues. Therefore, our paper is organized as follows. In Section 2, by Itô formula and proper Lyapunov function, it is shown that there is a unique global positive solution $\left(x\right(t),y(t\left)\right)$ to system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) for any positive initial value $\left(x\right(0),y(0\left)\right)\in R_{+}^2$. In Section 3, the stochastic ultimate boundedness of solution is studied. In Section 4, it is proved that system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) is stochastically permanent under Assumption A. In Section 5, some sufficient conditions for permanence in mean and extinction of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) are obtained. Finally, we present numerical examples to verify our theoretical results.

In this paper, unless otherwise specified, we let $(\mathrm{\Omega },\mathcal{F},\{{\mathcal{F}}_t{\}}_{t⩾0},P)$ be a complete probability space with a filtration $\{{\mathcal{F}}_t{\}}_{t⩾0}$ satisfying the usual conditions (i.e. it is increasing and right continuous all P-null sets). Also let $R_{+}^2=\left\{\right(x\left(t\right),y\left(t\right)):x(t)>0,y(t)>0\}$, $\Vert \left(x\right(t),y(t\left)\right)\Vert =\sqrt{x^2\left(t\right)+y^2\left(t\right)}$.

2. Positive and global solutions

It is necessary that proving the population density of prey and predator is non-negative. In general, in order to prove that there exists a unique global solution for any given initial value for a stochastic differential equation, the coefficients of the equation need to satisfy the linear growth condition and the local Lipschitz condition. However, the coefficients of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) neither satisfy the linear growth condition, nor local Lipschitz condition. In this section, by making the change of variables and applying Itô formula, we will firstly show there is a unique positive local solution for system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ).

Lemma 2.1

For any given initial value $\left(x\right(0),y(0\left)\right)\in R_{+}^2$, there is a unique local positive solution $\left(x\right(t),y(t\left)\right)$ to system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) for $t\in [0,{\tau }_e)$ a.s., where ${\tau }_e$ is the explosion time.

Proof.

Letting $u\left(t\right)=\ln x\left(t\right),v\left(t\right)=\ln y\left(t\right)$, and applying Itô formula, we get the following system (5) $\{\begin{array}{l}{\displaystyle \mathrm{d}u\left(t\right)=[a-\frac{{\sigma }_1^2}{2}-b{e}^{u\left(t\right)}-\frac{(\lambda +\alpha {\mathrm{e}}^{v\left(t\right)}){\mathrm{e}}^{v\left(t\right)}}{1+h(\lambda +\alpha {\mathrm{e}}^{v\left(t\right)}){\mathrm{e}}^{u\left(t\right)}}]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}v\left(t\right)=(r-\frac{{\sigma }_2^2}{2}-\frac{f{e}^{v\left(t\right)}}{m+{\mathrm{e}}^{u\left(t\right)}})\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2\left(t\right),}\end{array}$(5) and the initial value $u\left(0\right)=\ln x\left(0\right),v\left(0\right)=\ln y\left(0\right)$. Obviously, the coefficients of system (5(5) $\{\begin{array}{l}{\displaystyle \mathrm{d}u\left(t\right)=[a-\frac{{\sigma }_1^2}{2}-b{e}^{u\left(t\right)}-\frac{(\lambda +\alpha {\mathrm{e}}^{v\left(t\right)}){\mathrm{e}}^{v\left(t\right)}}{1+h(\lambda +\alpha {\mathrm{e}}^{v\left(t\right)}){\mathrm{e}}^{u\left(t\right)}}]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}v\left(t\right)=(r-\frac{{\sigma }_2^2}{2}-\frac{f{e}^{v\left(t\right)}}{m+{\mathrm{e}}^{u\left(t\right)}})\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2\left(t\right),}\end{array}$(5) ) are locally Lipschitz continuous, so there exists a unique local solution $\left(u\right(t),v(t\left)\right)$ of system (5(5) $\{\begin{array}{l}{\displaystyle \mathrm{d}u\left(t\right)=[a-\frac{{\sigma }_1^2}{2}-b{e}^{u\left(t\right)}-\frac{(\lambda +\alpha {\mathrm{e}}^{v\left(t\right)}){\mathrm{e}}^{v\left(t\right)}}{1+h(\lambda +\alpha {\mathrm{e}}^{v\left(t\right)}){\mathrm{e}}^{u\left(t\right)}}]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}v\left(t\right)=(r-\frac{{\sigma }_2^2}{2}-\frac{f{e}^{v\left(t\right)}}{m+{\mathrm{e}}^{u\left(t\right)}})\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2\left(t\right),}\end{array}$(5) ) on $t\in [0,{\tau }_e)$, where ${\tau }_e$ is the explosion time. Therefore, there is also a unique local positive solution $\left(x\right(t),y(t\left)\right)$ of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) for any given initial value $\left(x\right(0),y(0\left)\right)\in R_{+}^2$ on $t\in [0,{\tau }_e)$, here $x\left(t\right)={\mathrm{e}}^{u\left(t\right)},y\left(t\right)={\mathrm{e}}^{v\left(t\right)}$.

Next, we will prove system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) has a global positive solution for any given initial value $\left(x\right(0),y(0\left)\right)\in R_{+}^2$.

Theorem 2.2

For any given initial value $\left(x\right(0),y(0\left)\right)\in R_{+}^2$, there is a unique solution $\left(x\right(t),y(t\left)\right)$ for system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) on $t\ge 0$, and the solution will remain in $R_{+}^2$ with probability one, namely $\left(x\right(t),y(t\left)\right)\in R_{+}^2$ for all $t\ge 0$ almost surely.

Proof.

By Lemma 2.1, we know that there is a unique local positive solution for system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ). Further, we show that this solution is global, i.e. ${\tau }_e=\mathrm{\infty }$. Let $n_0\ge 0,$ and $n_0$ be sufficiently large such that $x\left(0\right)$ and $y\left(0\right)$ all lie within the interval $[\frac{1}{n_0},n_0]$. For each integer $n\ge n_0$, define the stopping time: ${\tau }_n=inf\{t\in [0,{\tau }_e):x(t)\notin (\frac{1}{n},n)\mathrm{or}y(t)\notin (\frac{1}{n},n)\},$where $inf\varnothing =\mathrm{\infty }$ (∅ is the empty set). Obviously, ${\tau }_n$ is increasing as $n\to \mathrm{\infty }$. Let ${\tau }_{\mathrm{\infty }}=\underset{n\to \mathrm{\infty }}{lim}{\tau }_n$. Easy to get ${\tau }_{\mathrm{\infty }}\le {\tau }_{e}a.s.$ If we can show ${\tau }_{\mathrm{\infty }}=\mathrm{\infty }$ $a.s.$, then ${\tau }_e=\mathrm{\infty }$ $a.s.$ Now, we only need to prove ${\tau }_{\mathrm{\infty }}=\mathrm{\infty }$ $a.s.$ If not, we assume ${\tau }_{\mathrm{\infty }}\ne \mathrm{\infty }$ a.s., there are two constants T>0 and $ϵ\in (0,1)$, such that $P\{{\tau }_{\mathrm{\infty }}\le T\}\ge ϵ.$ Therefore, there exists an integer $n_1(\ge n_0),$ such that when $n\ge n_1$, we have (6) $P\{{\tau }_n\le T\}\ge ϵ.$(6) In order to finish our proof, we will construct $C^2$-function $V_1:R_{+}^2\to R_{+}$ under two different cases.

Case 1: If $\mathrm{h\lambda m}\le 1$, we define a function $V_1\left(x\right(t),y(t\left)\right)$ as follows: $V_1\left(x\right(t),y(t\left)\right)=\mathrm{fh\lambda }\left[x\right(t)-\ln x(t)-1]+\alpha \left[y\right(t)-\ln y(t)-1].$Since $u-\ln u-1\ge 0$ for all u>0, it is known that the function $V_1:R_{+}^2\to R_{+}$. Letting $n\ge n_1$, for $0\le t\le {\tau }_n$, and applying Itô formula to $V_1\left(x\right(t),y(t\left)\right)$, we obtain $\mathrm{d}{V}_1\left(x\right(t),y(t\left)\right)=L{V}_1\left(x\right(t),y(t\left)\right)\mathrm{d}t+\mathrm{fh\lambda }{\sigma }_1\left(x\right(t)-1)\mathrm{d}{B}_1\left(t\right)+\alpha {\sigma }_2\left(y\right(t)-1)\mathrm{d}{B}_2\left(t\right),$here $\begin{array}{rl}L{V}_1\left(x\right(t),y(t\left)\right)& =\mathrm{fh\lambda x}\left(t\right)(a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)})\\ & -\mathrm{fh\lambda }(a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)})+\frac{\mathrm{fh\lambda }{\sigma }_1^2}{2}\\ & +\mathrm{\alpha y}\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})-\alpha (r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})+\frac{\alpha {\sigma }_2^2}{2}\\ & \le \mathrm{fh\lambda }\left[\right(a+b\left)x\right(t)-b{x}^2(t\left)\right]+\mathrm{fh}{\lambda }^{2}y\left(t\right)+\frac{\mathrm{fh\lambda \alpha }{y}^2\left(t\right)}{1+\mathrm{h\lambda x}\left(t\right)}-\frac{\mathrm{\alpha f}{y}^2\left(t\right)}{m+x\left(t\right)}\\ & +(-\mathrm{afh\lambda }-\mathrm{\alpha r}+\frac{\mathrm{fh\lambda }{\sigma }_1^2}{2}+\frac{\alpha {\sigma }_2^2}{2})+\alpha (r+\frac{f}{m})y\left(t\right)\\ & \le K_1+(\mathrm{fh}{\lambda }^2+\mathrm{\alpha r}+\frac{\mathrm{\alpha f}}{m})y\left(t\right)-\frac{\mathrm{\alpha f}(1-\mathrm{h\lambda m})y^2\left(t\right)}{[1+\mathrm{h\lambda x}(t\left)\right][m+x(t\left)\right]}\\ & \le K_1+K_{2}y\left(t\right),\end{array}$where $K_1$ and $K_2$ are positive constants. Since $y\left(t\right)\le 2\left[y\right(t)-1-\ln y(t\left)\right]+\ln 4\le \frac{2}{\alpha }{V}_1\left(x\right(t),y(t\left)\right)+\ln 4$ for $y\left(t\right)>0$, we have $L{V}_1\left(x\right(t),y(t\left)\right)\le K_1+K_2\ln 4+\frac{2}{\alpha }{K}_2{V}_1\left(x\right(t),y(t\left)\right)\le K_3(1+V_1\left(x\right(t),y(t\left)\right)),$where $K_3=max\{K_1+K_2\ln 4,\frac{2}{\alpha }{K}_2\}$, consequently $\begin{array}{rl}\mathrm{d}{V}_1\left(x\right(t),y(t\left)\right)& \le K_3(1+V_1\left(x\right(t),y(t\left)\right))\mathrm{d}t+\mathrm{fh\lambda }{\sigma }_1\left(x\right(t)-1)\mathrm{d}{B}_1\left(t\right)\\ & +\alpha {\sigma }_2\left(y\right(t)-1)\mathrm{d}{B}_2\left(t\right).\end{array}$Therefore, integrating sides of the above inequality from 0 to ${\tau }_n\wedge T$ and then taking expectation, we get $\begin{array}{rl}& E{V}_1\left(x\right({\tau }_n\wedge T),y({\tau }_n\wedge T\left)\right)\\ & \le V_1\left(x\right(0),y(0\left)\right)+K_{3}E{\int }_0^{{\tau }_n\wedge T}(1+V_1\left(x\right(t),y(t\left)\right))\mathrm{d}t\\ & \le V_1\left(x\right(0),y(0\left)\right)+K_{3}T+K_{3}E{\int }_0^{{\tau }_n\wedge T}{V}_1\left(x\right(t),y(t\left)\right)\mathrm{d}t\\ & \le V_1\left(x\right(0),y(0\left)\right)+K_{3}T+K_3{\int }_0^{T}E{V}_1\left(x\right({\tau }_n\wedge t),y({\tau }_n\wedge t\left)\right)\mathrm{d}t.\end{array}$Further, by Gronwall inequality, we get $E{V}_1\left(x\right({\tau }_n\wedge T),y({\tau }_n\wedge T\left)\right)\le K_4,$where $K_4=\left(V_1\right(x\left(0\right),y\left(0\right))+K_{3}T)){\mathrm{e}}^{K_{3}T}.$ The remainder of the proof is similar to Theorem 1 in [8], we omit the details.

Case 2: If $\mathrm{h\lambda m}\ge 1$, we define a $C^2$-function $V_1\left(x\right(t),y(t\left)\right)$ as follows: $V_1\left(x\right(t),y(t\left)\right)=f\left[x\right(t)-\ln x(t)-1]+\mathrm{\alpha m}\left[y\right(t)-\ln y(t)-1].$Applying Itô formula to $V_1\left(x\right(t),y(t\left)\right)$, we obtain that $\mathrm{d}{V}_1\left(x\right(t),y(t\left)\right)=L{V}_1\left(x\right(t),y(t\left)\right)\mathrm{d}t+f{\sigma }_1\left(x\right(t)-1)\mathrm{d}{B}_1\left(t\right)+\mathrm{\alpha m}{\sigma }_2\left(y\right(t)-1)\mathrm{d}{B}_2\left(t\right),$where $\begin{array}{rl}& L{V}_1\left(x\right(t),y(t\left)\right)\\ & \le f\left[\right(a+b\left)x\right(t)-b{x}^2(t\left)\right]+\frac{\mathrm{f\lambda y}\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}+\frac{\mathrm{f\alpha }{y}^2\left(t\right)}{1+\mathrm{h\lambda x}\left(t\right)}\\ & -\frac{\mathrm{\alpha mf}{y}^2\left(t\right)}{m+x\left(t\right)}+\mathrm{\alpha mry}\left(t\right)+(-\mathrm{af}-\mathrm{\alpha mr}+\frac{f{\sigma }_1^2}{2}+\frac{\mathrm{\alpha m}{\sigma }_2^2}{2})+\frac{\mathrm{\alpha mfy}\left(t\right)}{m+x\left(t\right)}.\\ & \le K_1+(\mathrm{f\lambda }+\mathrm{\alpha mr}+\mathrm{\alpha f})y\left(t\right)-\frac{\mathrm{f\alpha }(\mathrm{h\lambda m}-1)x\left(t\right)y^2\left(t\right)}{[1+\mathrm{h\lambda x}(t\left)\right][m+x(t\left)\right]}\\ & \le K_1+K_{2}y\left(t\right).\end{array}$Here $K_1$ and $K_2$ are positive constants. Similar to the proof under Case 1, we obtain the existence of a global and positive solution.

3. Stochastically ultimate boundedness

Definition 3.1

[13]

The solution $\left(x\right(t),y(t\left)\right)$ of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) is said to be stochastically ultimately bounded, if for any $ϵ\in (0,1)$, there is a positive constant $H=H\left(ϵ\right)$, such that for any given initial value $\left(x\right(0),y(0\left)\right)\in R_{+}^2$, the solution $\left(x\right(t),y(t\left)\right)$ of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) has theproperty that (7) $\underset{t\to \mathrm{\infty }}{lim sup}P\{\Vert \left(x\right(t),y(t\left)\right)\Vert >H\}<ϵ.$(7)

Before proving boundedness, we first prove the following lemma.

Lemma 3.2

There exists a constant $C_1>0$ independent of the initial value $\left(x\right(0),y(0\left)\right)\in R_{+}^2$, such that the solution $\left(x\right(t),y(t\left)\right)$ of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) has the following property (8) $\underset{t\to \mathrm{\infty }}{lim sup}E\Vert \left(x\right(t),y(t\left)\right){\Vert }^2\le C_1.$(8)

Proof.

Define (9) $V_2\left(x\right(t),y(t\left)\right)=x^2\left(t\right)+y^2\left(t\right)+(H_1+x(t\left)\right)y^2\left(t\right),$(9) where $H_1=max\{0,m-1\}$. Hence, for $\left(x\right(t),y(t\left)\right)\in R_{+}^2$, by the Itô formula, we have (10) $\begin{array}{rl}\mathrm{d}{V}_2\left(x\right(t),y(t\left)\right)& =L{V}_2\left(x\right(t),y(t\left)\right)\mathrm{d}t+\left(2x\right(t)+y^2(t\left)\right){\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right)\\ & +2(1+H_1+x(t\left)\right){\sigma }_2{y}^2\left(t\right)\mathrm{d}{B}_2\left(t\right).\end{array}$(10) Here $\begin{array}{rl}& L{V}_2\left(x\right(t),y(t\left)\right)\\ & =2x^2\left(t\right)(a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)})+{\sigma }_1^2{x}^2\left(t\right)\\ & +2y^2\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})+{\sigma }_2^2{y}^2\left(t\right)+(H_1+x(t\left)\right)2y^2\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\\ & +x\left(t\right)y^2\left(t\right)(a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)})+(H_1+x(t\left)\right){\sigma }_2^2{y}^2\left(t\right)\\ & \le 2a{x}^2\left(t\right)-2b{x}^3\left(t\right)+{\sigma }_1^2{x}^2\left(t\right)+2r{y}^2\left(t\right)+{\sigma }_2^2{y}^2\left(t\right)-\frac{2f{y}^3\left(t\right)(H_1+1+x(t\left)\right)}{m+x\left(t\right)}\\ & +2r(H_1+x(t\left)\right)y^2\left(t\right)+(H_1+x(t\left)\right){\sigma }_2^2{y}^2\left(t\right)+\mathrm{ax}\left(t\right)y^2\left(t\right)-b{x}^2\left(t\right)y^2\left(t\right).\end{array}$Noting that $\frac{2f{y}^3\left(t\right)(H_1+1+x(t\left)\right)}{m+x\left(t\right)}=\frac{2f{y}^3\left(t\right)(1+x(t\left)\right)}{m+x\left(t\right)}\ge 2f{y}^3\left(t\right)\mathrm{for}m\le 1,$and $\frac{2f{y}^3\left(t\right)(H_1+1+x(t\left)\right)}{m+x\left(t\right)}=2f{y}^3\left(t\right)\mathrm{for}m>1,$it is followed by $\begin{array}{rl}& L{V}_2\left(x\right(t),y(t\left)\right)\\ & \le (2a+{\sigma }_1^2)x^2\left(t\right)-2b{x}^3\left(t\right)+(2r+{\sigma }_2^2)y^2-2f{y}^3\left(t\right)\\ & +2r(H_1+x(t\left)\right)y^2\left(t\right)+(H_1+x(t\left)\right){\sigma }_2^2{y}^2\left(t\right)+\mathrm{ax}\left(t\right)y^2\left(t\right)-b{x}^2\left(t\right)y^2\left(t\right)\\ & =-V_2\left(x\right(t),y(t\left)\right)+(2a+{\sigma }_1^2+1)x^2\left(t\right)-2b{x}^3\left(t\right)+(2r+{\sigma }_2^2+1)y^2\left(t\right)-2f{y}^3\left(t\right)\\ & +(2r+1)(H_1+x(t\left)\right)y^2\left(t\right)+(H_1+x(t\left)\right){\sigma }_2^2{y}^2\left(t\right)+\mathrm{ax}\left(t\right)y^2\left(t\right)-b{x}^2\left(t\right)y^2\left(t\right)\\ & \le -V_2\left(x\right(t),y(t\left)\right)+(2a+{\sigma }_1^2+1)x^2\left(t\right)-2b{x}^3\left(t\right)+(2r+{\sigma }_2^2+1)y^2\left(t\right)\\ & -2f{y}^3\left(t\right)+K{y}^2\left(t\right),\end{array}$where $K=\underset{\{x\ge 0\}}{max}\left\{\right(2r+1+{\sigma }_2^2\left)\right(H_1+x\left(t\right))+\mathrm{ax}(t)-b{x}^2(t\left)\right\}$. Further, there is a constant $H_2>0$, such that $L{V}_2\left(x\right(t),y(t\left)\right)\le H_2-V_2\left(x\right(t),y(t\left)\right).$Substituting the above inequality into (10(10) $\begin{array}{rl}\mathrm{d}{V}_2\left(x\right(t),y(t\left)\right)& =L{V}_2\left(x\right(t),y(t\left)\right)\mathrm{d}t+\left(2x\right(t)+y^2(t\left)\right){\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right)\\ & +2(1+H_1+x(t\left)\right){\sigma }_2{y}^2\left(t\right)\mathrm{d}{B}_2\left(t\right).\end{array}$(10) ) gives $\begin{array}{rl}\mathrm{d}{V}_2\left(x\right(t),y(t\left)\right)& \le [H_2-V_2\left(x\right(t),y(t\left)\right)]\mathrm{d}t+\left(2x\right(t)+y^2(t\left)\right){\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right)\\ & +2(1+H_1+x(t\left)\right){\sigma }_2{y}^2\left(t\right)\mathrm{d}{B}_2\left(t\right).\end{array}$Once again by the Itô formula, we get $\begin{array}{rl}\mathrm{d}\left[{\mathrm{e}}^t{V}_2\left(x\right(t),y(t\left)\right)\right]& ={\mathrm{e}}^t[V_2\left(x\right(t),y(t\left)\right)\mathrm{d}t+\mathrm{d}{V}_2\left(x\right(t),y(t\left)\right)]\\ & \le H_2{\mathrm{e}}^t\mathrm{d}t+{\mathrm{e}}^t\left(2x\right(t)+y^2(t\left)\right){\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right)\\ & +2{\mathrm{e}}^t(1+H_1+x(t\left)\right){\sigma }_2{y}^2\left(t\right)\mathrm{d}{B}_2\left(t\right).\end{array}$Integrating both sides of the above inequality from 0 to t, and then taking expectations, we get ${\mathrm{e}}^{t}E{V}_2\left(x\right(t),y(t\left)\right)\le V_2\left(x\right(0),y(0\left)\right)+H_2({\mathrm{e}}^t-1),$This implies that $\underset{t\to \mathrm{\infty }}{lim sup}E{V}_2\left(x\right(t),y(t\left)\right)\le H_2.$On the other hand, we have $\Vert \left(x\right(t),y(t\left)\right){\Vert }^2\le 2max\left\{x^2\right(t),y^2(t\left)\right\}\le 2V_2\left(x\right(t),y(t\left)\right),$Therefore, we can obtain $\underset{t\to \mathrm{\infty }}{lim sup}E\Vert \left(x\right(t),y(t\left)\right){\Vert }^2\le 2\underset{t\to \mathrm{\infty }}{lim sup}E{V}_2\left(x\right(t),y(t\left)\right)\le 2H_2:=C_1.$

Based on the above conclusion, we can obtain the stochastically ultimate boundedness of the solution $\left(x\right(t),y(t\left)\right)$ of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ).

Theorem 3.3

The solution of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) is stochastically ultimately bounded.

Proof.

By Lemma 3.2, there is $C_1>0$ such that $\underset{t\to \mathrm{\infty }}{lim sup}E\Vert \left(x\right(t),y(t\left)\right){\Vert }^2\le C_1.$Now, for any $ϵ>0$, let $H=\frac{\sqrt{2C_1}}{\sqrt{ϵ}}$, then by Chebyshev inequality, we have $P\{\Vert \left(x\right(t),y(t\left)\right)\Vert >H\}\le \frac{E\Vert \left(x\right(t),y(t\left)\right){\Vert }^2}{H^2},$Hence $\underset{t\to \mathrm{\infty }}{lim sup}P\{\Vert \left(x\right(t),y(t\left)\right)\Vert >H\}\le C_1\cdot \frac{ϵ}{2C_1}=\frac{ϵ}{2}<ϵ.$By Definition 3.1, the solution of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) is stochastically ultimately bounded.

4. Stochastic permanence

Definition 4.1

[14]

The solution $X\left(t\right)=\left(x\right(t),y(t\left)\right)$ of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) is said to be stochastically permanent, if for any $ϵ\in (0,1)$, there exists a pair of positive constants $M_1=M_1\left(ϵ\right)$ and $M_2=M_2\left(ϵ\right)$, such that for any initial value $X\left(0\right)=\left(x\right(0),y(0\left)\right)\in R_{+}^2$, the solution $X\left(t\right)=\left(x\right(t),y(t\left)\right)$ to system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) has the properties that $\begin{array}{rl}\underset{t\to \mathrm{\infty }}{lim inf}P\{\Vert X(t)\Vert \ge M_1\}& \ge 1-ϵ,\\ \underset{t\to \mathrm{\infty }}{lim inf}P\{\Vert X(t)\Vert \le M_2\}& \ge 1-ϵ.\end{array}$

For convenient of discussion, we give the following assumption.

Assumption A

$\frac{1}{2}max\{{\sigma }_1^2,{\sigma }_2^2\}<min\{a,r-\frac{1}{h}\}$.

Theorem 4.2

If Assumption A holds, there exists a positive constant $F(\theta ,p)$, which is independent of the initial value $X\left(0\right)=\left(x\right(0),y(0\left)\right)\in R_{+}^2$, such that the solution $X\left(t\right)$ of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) has the property that (11) $\underset{t\to \mathrm{\infty }}{lim sup}E\left(\frac{1}{\Vert X\left(t\right){\Vert }^{\theta }}\right)\le F,$(11) where $0<\theta <2$ is an arbitrary positive constant satisfying (12) $\frac{\theta +1}{2}max\{{\sigma }_1^2,{\sigma }_2^2\}<min\{a,r-\frac{1}{h}\},$(12) and p>0 satisfying (13) $p+\frac{\theta (\theta +1)}{2}max\{{\sigma }_1^2,{\sigma }_2^2\}<\theta min\{a,r-\frac{1}{h}\}.$(13)

Proof.

Define $U\left(X\right(t\left)\right)=U\left(x\right(t),y(t\left)\right)=\frac{1}{x\left(t\right)+y\left(t\right)}$, by the Itô formula, we obtain $\begin{array}{rl}\mathrm{d}U\left(X\right(t\left)\right)& =-U^2\left(X\right(t\left)\right)\left[x\right(t)(a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)})\\ & +y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})]\mathrm{d}t+U^3\left(X\right(t\left)\right)\left[{\sigma }_1^2{x}^2\right(t)+{\sigma }_2^2{y}^2(t\left)\right]\mathrm{d}t\\ & -U^2\left(X\right(t\left)\right)\left[{\sigma }_{1}x\right(t\left)\mathrm{d}{B}_1\right(t)+{\sigma }_{2}y(t\left)\mathrm{d}{B}_2\right(t\left)\right]\\ & =\mathrm{LU}\left(X\right(t\left)\right)\mathrm{d}t-U^2\left(X\right(t\left)\right)\left({\sigma }_{1}x\right(t\left)\mathrm{d}{B}_1\right(t)+{\sigma }_{2}y(t\left)\mathrm{d}{B}_2\right(t\left)\right),\end{array}$where $\begin{array}{rl}\mathrm{LU}\left(X\right(t\left)\right)& =-U^2\left(X\right(t\left)\right)\left[x\right(t)(a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)})\\ & +y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})]+U^3\left(X\right(t\left)\right)\left[{\sigma }_1^2{x}^2\right(t)+{\sigma }_2^2{y}^2(t\left)\right].\end{array}$Under Assumption A, there certainly exists a constant θ such that it satisfies (12(12) $\frac{\theta +1}{2}max\{{\sigma }_1^2,{\sigma }_2^2\}<min\{a,r-\frac{1}{h}\},$(12) ). And denote $W\left(X\right(t\left)\right)=(1+U(X\left(t\right)){)}^{\theta }$, then by the Itô formula, we get $\begin{array}{rl}\mathrm{d}W\left(X\right(t\left)\right)& =\mathrm{LW}\left(X\right(t\left)\right)\mathrm{d}t-\theta {(1+U\left(X\right(t\left)\right))}^{\theta -1}{U}^2\left(X\right(t\left)\right)\\ & \cdot \left[{\sigma }_{1}x\right(t\left)\mathrm{d}{B}_1\right(t)+{\sigma }_{2}y(t\left)\mathrm{d}{B}_2\right(t\left)\right],\end{array}$where $\begin{array}{rl}\mathrm{LW}\left(X\right(t\left)\right)& =-\theta {(1+U\left(X\right(t\left)\right))}^{\theta -1}{U}^2\left(X\right(t\left)\right)x\left(t\right)(a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)})\\ & -\theta {(1+U\left(X\right(t\left)\right))}^{\theta -1}{U}^2\left(X\right(t\left)\right)y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\\ & +\theta {(1+U\left(X\right(t\left)\right))}^{\theta -1}{U}^3\left(X\right(t\left)\right)\left[{\sigma }_1^2{x}^2\right(t)+{\sigma }_2^2{y}^2(t\left)\right]\\ & +\frac{\theta (\theta -1)}{2}{(1+U\left(X\right(t\left)\right))}^{\theta -2}{U}^4\left(X\right(t\left)\right)\left[{\sigma }_1^2{x}^2\right(t)+{\sigma }_2^2{y}^2(t\left)\right].\end{array}$By (13(13) $p+\frac{\theta (\theta +1)}{2}max\{{\sigma }_1^2,{\sigma }_2^2\}<\theta min\{a,r-\frac{1}{h}\}.$(13) ), we can choose p>0 sufficiently small such that it satisfies $p+\frac{\theta (\theta +1)}{2}max\{{\sigma }_1^2,{\sigma }_2^2\}<\theta min\{a,r-\frac{1}{h}\}$. Thus by the Itô formula, we compute $\begin{array}{rl}& L\left({\mathrm{e}}^{\mathrm{pt}}W\right(X\left(t\right)\left)\right)\\ & ={\mathrm{e}}^{\mathrm{pt}}\mathrm{LW}\left(X\right(t\left)\right)+p{e}^{\mathrm{pt}}W\left(X\right(t\left)\right)\\ & =-{\mathrm{e}}^{\mathrm{pt}}\theta {(1+U\left(X\right(t\left)\right))}^{\theta -1}{U}^2\left(X\right(t\left)\right)x\left(t\right)(a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)})\\ & -{\mathrm{e}}^{\mathrm{pt}}\theta {(1+U\left(X\right(t\left)\right))}^{\theta -1}{U}^2\left(X\right(t\left)\right)y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\\ & +{\mathrm{e}}^{\mathrm{pt}}\theta {(1+U\left(X\right(t\left)\right))}^{\theta -1}{U}^3\left(X\right(t\left)\right)\left[{\sigma }_1^2{x}^2\right(t)+{\sigma }_2^2{y}^2(t\left)\right]\\ & +\frac{\theta (\theta -1)}{2}{\mathrm{e}}^{\mathrm{pt}}{(1+U\left(X\right(t\left)\right))}^{\theta -2}{U}^4\left(X\right(t\left)\right)\left[{\sigma }_1^2{x}^2\right(t)+{\sigma }_2^2{y}^2(t\left)\right]\\ & +p{e}^{\mathrm{pt}}{(1+U\left(X\right(t\left)\right))}^{\theta }\\ & ={\mathrm{e}}^{\mathrm{pt}}{(1+U\left(X\right(t\left)\right))}^{\theta -2}[p{(1+U\left(X\right(t\left)\right))}^2+G\left(X\right(t\left)\right)],\end{array}$where $\begin{array}{rl}G\left(X\right(t\left)\right)& =-\theta (1+U\left(X\right(t\left)\right))U^2\left(X\right(t\left)\right)\left[x\right(t)(a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)})\\ & +y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})]+\frac{\theta (\theta +1)}{2}{U}^4\left(X\right(t\left)\right)\left[{\sigma }_1^2{x}^2\right(t)+{\sigma }_2^2{y}^2(t\left)\right]\\ & +\theta U^3\left(X\right(t\left)\right)\left[{\sigma }_1^2{x}^2\right(t)+{\sigma }_2^2{y}^2(t\left)\right].\end{array}$It is easy to get that ${\sigma }_1^2{x}^2\left(t\right)+{\sigma }_2^2{y}^2\left(t\right)\le max\{{\sigma }_1^2,{\sigma }_2^2\}\left(x^2\right(t)+y^2(t\left)\right)\le max\{{\sigma }_1^2,{\sigma }_2^2\}\frac{1}{U^2\left(X\right(t\left)\right)},$and $\begin{array}{rl}& x\left(t\right)(a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)})+y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\\ & \ge \mathrm{ax}\left(t\right)-b{x}^2\left(t\right)-\frac{x\left(t\right)y\left(t\right)}{\mathrm{hx}\left(t\right)}+\mathrm{ry}\left(t\right)-\frac{f}{m}{y}^2\left(t\right)\\ & =\mathrm{ax}\left(t\right)+(r-\frac{1}{h})y\left(t\right)-b{x}^2\left(t\right)-\frac{f}{m}{y}^2\left(t\right)\\ & \ge min\{a,r-\frac{1}{h}\}\frac{1}{U\left(X\right(t\left)\right)}-max\{b,\frac{f}{m}\}\frac{1}{U^2\left(X\right(t\left)\right)}.\end{array}$Therefore $\begin{array}{rl}G\left(X\right(t\left)\right)& \le -\theta (1+U\left(X\right(t\left)\right))U^2\left(X\right(t\left)\right)[min\{a,r-\frac{1}{h}\}\frac{1}{U\left(X\right(t\left)\right)}\\ & -max\{b,\frac{f}{m}\}\frac{1}{U^2\left(X\right(t\left)\right)}]+\frac{\theta (\theta +1)}{2}{U}^2\left(X\right(t\left)\right)max\{{\sigma }_1^2,{\sigma }_2^2\}\\ & +\mathrm{\theta U}\left(X\right(t\left)\right)max\{{\sigma }_1^2,{\sigma }_2^2\}.\end{array}$By $\theta min\{a,r-\frac{1}{h}\}-\frac{\theta (\theta +1)}{2}max\{{\sigma }_1^2,{\sigma }_2^2\}-p>0$, it is no difficult to see that there is a constant $F_1>0$ such that $\begin{array}{rl}& L\left({\mathrm{e}}^{\mathrm{pt}}W\right(X\left(t\right)\left)\right)\\ & ={\mathrm{e}}^{\mathrm{pt}}{(1+U\left(X\right(t\left)\right))}^{\theta -2}(p{(1+U\left(X\right(t\left)\right))}^2+G\left(X\right(t\left)\right))\\ & \le {\mathrm{e}}^{\mathrm{pt}}{(1+U\left(X\right(t\left)\right))}^{\theta -2}[p+\theta max\{b,\frac{f}{m}\}\\ & +U\left(X\right(t\left)\right)(2p-\theta min\{a,r-\frac{1}{h}\}+\theta max\{b,\frac{f}{m}\}+\theta max\{{\sigma }_1^2,{\sigma }_2^2\})\\ & -U^2\left(X\right(t\left)\right)(\theta min\{a,r-\frac{1}{h}\}-\frac{\theta (\theta +1)}{2}max\{{\sigma }_1^2,{\sigma }_2^2\}-p)]\\ & \le {\mathrm{e}}^{\mathrm{pt}}{(1+U\left(X\right(t\left)\right))}^{\theta -2}[p+\theta max\{b,\frac{f}{m}\}\\ & +U\left(X\right(t\left)\right)(2p+\theta max\{b,\frac{f}{m}\}+\theta max\{{\sigma }_1^2,{\sigma }_2^2\})\\ & -U^2\left(X\right(t\left)\right)(\theta min\{a,r-\frac{1}{h}\}-\frac{\theta (\theta +1)}{2}max\{{\sigma }_1^2,{\sigma }_2^2\}-p)]\\ & \le F_1{\mathrm{e}}^{\mathrm{pt}}.\end{array}$Then $\begin{array}{rl}& \mathrm{d}\left({\mathrm{e}}^{\mathrm{pt}}W\right(X\left(t\right)\left)\right)=L\left({\mathrm{e}}^{\mathrm{pt}}W\right(X\left(t\right)\left)\right)\\ & -{\mathrm{e}}^{\mathrm{pt}}\theta {(1+U\left(X\right(t\left)\right))}^{\theta -1}{U}^2\left(X\right(t\left)\right)\left({\sigma }_{1}x\right(t\left)\mathrm{d}{B}_1\right(t)+{\sigma }_{2}y(t\left)\mathrm{d}{B}_2\right(t\left)\right)\\ & \le F_1{\mathrm{e}}^{\mathrm{pt}}\mathrm{d}t-{\mathrm{e}}^{\mathrm{pt}}\theta {(1+U\left(X\right(t\left)\right))}^{\theta -1}{U}^2\left(X\right(t\left)\right)\left({\sigma }_{1}x\right(t\left)\mathrm{d}{B}_1\right(t)+{\sigma }_{2}y(t\left)\mathrm{d}{B}_2\right(t\left)\right).\end{array}$Integrating both sides of the above inequality from 0 to t, and then taking the expectation, we obtain that $\begin{array}{rl}{\mathrm{e}}^{\mathrm{pt}}E{(1+U(X\left(t\right)\left)\right)}^{\theta }& \le {(1+U(X\left(0\right)\left)\right)}^{\theta }+\frac{F_1({\mathrm{e}}^{\mathrm{pt}}-1)}{p}\\ & \le {(1+U(X\left(0\right)\left)\right)}^{\theta }+F_2{\mathrm{e}}^{\mathrm{pt}},\end{array}$where $F_2=\frac{F_1}{p}$, so $\underset{t\to \mathrm{\infty }}{lim sup}E{U}^{\theta }\left(X\right(t\left)\right)\le \underset{t\to \mathrm{\infty }}{lim sup}E{(1+U\left(X\right(t\left)\right))}^{\theta }\le F_2.$Note that ${\left(x\right(t)+y(t\left)\right)}^{\theta }\le 2^{\theta }{\left(x^2\right(t)+y^2(t\left)\right)}^{\frac{\theta }{2}}=2^{\theta }\Vert X\left(t\right){\Vert }^{\theta },$Therefore $\underset{t\to \mathrm{\infty }}{lim sup}E\left(\frac{1}{\Vert X\left(t\right){\Vert }^{\theta }}\right)\le 2^{\theta }\underset{t\to \mathrm{\infty }}{lim sup}E{U}^{\theta }\left(X\right(t\left)\right)\le 2^{\theta }{F}_2:=F.$

Theorem 4.3

Under Assumption A, system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) is stochastically permanent.

Proof.

By Theorem 4.2, there exists F>0 such that $\underset{t\to \mathrm{\infty }}{lim sup}E\left(\frac{1}{\Vert X\left(t\right){\Vert }^{\theta }}\right)\le F.$Thus, for any $ϵ>0$, let $M_1={\left(\frac{ϵ}{2F}\right)}^{\frac{1}{\theta }}$, then by Chebyshev inequality, we obtain $P\{\Vert X(t)\Vert <M_1\}=P\{\frac{1}{\Vert X\left(t\right){\Vert }^{\theta }}>\frac{1}{M_1^{\theta }}\}\le M_1^{\theta }E\left(\frac{1}{\Vert X\left(t\right){\Vert }^{\theta }}\right).$It is straightforward to see that $\underset{t\to \mathrm{\infty }}{lim sup}P\{\Vert X(t)\Vert <M_1\}\le \frac{ϵ}{2F}\cdot F=\frac{ϵ}{2}<ϵ.$In other words $\underset{t\to \mathrm{\infty }}{lim inf}P\{\Vert X(t)\Vert \ge M_1\}\ge 1-ϵ.$On the other hand, by Theorem 3.3, we can see that there is a constant $M_2=H=\frac{\sqrt{2C_1}}{\sqrt{ϵ}}$, such that $\underset{t\to \mathrm{\infty }}{lim sup}P\{\Vert X(t)\Vert >M_2\}<ϵ,$that is $\underset{t\to \mathrm{\infty }}{lim inf}P\{\Vert X(t)\Vert \le M_2\}\ge 1-ϵ.$By Definition 4.1, the solution of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) is stochastically permanent.

5. Persistent in mean and extinction

In this section, we discuss the persistence in mean and extinction for system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ).

Definition 5.1

[15]

(1)	The population $X\left(t\right)$ is said to be persistent in mean if $\underset{t\to \mathrm{\infty }}{lim inf}\frac{1}{t}{\int }_0^{t}X\left(s\right)\mathrm{d}s>0,a.s.$
(2)	The population $X\left(t\right)$ is said to be extinct a.s. if $\underset{t\to \mathrm{\infty }}{lim}X\left(t\right)=0,a.s.$

Theorem 5.2

Suppose $\left(x\right(t),y(t\left)\right)$ be a solution of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) with the initial value $\left(x\right(0),y(0\left)\right)\in R_{+}^2$, then:

(i)	If $a>\frac{{\sigma }_1^2}{2}$ and $r<\frac{{\sigma }_2^2}{2}$, then x is persistent in mean and y is extinct a.s.
(ii)	If $a<\frac{{\sigma }_1^2}{2}$ and $r>\frac{{\sigma }_2^2}{2}$, then x is extinct a.s. and y is persistent in mean.
(iii)	If $a<\frac{{\sigma }_1^2}{2}$ and $r<\frac{{\sigma }_2^2}{2}$, then both the populations x and y are extinct a.s.

Proof.

(i) Applying Itô formula to $\ln y\left(t\right)$, we have $\mathrm{d}\ln y\left(t\right)=(r-\frac{{\sigma }_2^2}{2}-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2\left(t\right).$Integrating both sides from 0 to t the above equation, it yields $\ln y\left(t\right)=\ln y\left(0\right)+(r-\frac{{\sigma }_2^2}{2})t-f{\int }_0^t\frac{y\left(s\right)}{m+x\left(s\right)}\mathrm{d}s+{\sigma }_2{B}_2\left(t\right).$Further, we have (14) $\ln y\left(t\right)\le \ln y\left(0\right)+(r-\frac{{\sigma }_2^2}{2})t+{\sigma }_2{B}_2\left(t\right).$(14) Applying the strong law of large numbers of local martingales, we obtain $\underset{t\to \mathrm{\infty }}{lim}\frac{{\sigma }_2{B}_2\left(t\right)}{t}=0.$Then dividing by t on both sides of (14(14) $\ln y\left(t\right)\le \ln y\left(0\right)+(r-\frac{{\sigma }_2^2}{2})t+{\sigma }_2{B}_2\left(t\right).$(14) ), and letting $t\to \mathrm{\infty }$, we can easily derive $\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln y\left(t\right)}{t}\le r-\frac{{\sigma }_2^2}{2}<0,$therefore $\underset{t\to \mathrm{\infty }}{lim}y\left(t\right)=0a.s.$Hence for arbitrary ${ϵ}_1>0$, there exist $t_1$, such that $\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}\le (\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)\le {ϵ}_1$ for $t\ge t_1$. Thus, $x\left(t\right)(a-\mathrm{bx}(t)-{ϵ}_1)\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right)\le \mathrm{d}x\left(t\right)\le x\left(t\right)(a-\mathrm{bx}(t\left)\right)\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right).$Following Lemma A.1 of [16], and by Comparison Theorem for the stochastic equation, we have $\underset{t\to \mathrm{\infty }}{lim inf}\frac{1}{t}{\int }_0^{t}x\left(s\right)\mathrm{d}s\ge \frac{a-{ϵ}_1-\frac{{\sigma }_1^2}{2}}{b},\underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}{\int }_0^{t}x\left(s\right)\mathrm{d}s\le \frac{a-\frac{{\sigma }_1^2}{2}}{b}.$For the arbitrary of ${ϵ}_1$, we have $\underset{t\to \mathrm{\infty }}{lim}\frac{1}{t}{\int }_0^{t}x\left(s\right)\mathrm{d}s=\frac{a-\frac{{\sigma }_1^2}{2}}{b}>0.$Similar to the proof of the above (i) and Theorem 3.4 of [16], we may obtain conclusion (ii) and (iii). The details are omitted.

6. Numerical simulations

By using the Milstein method in [17], we can get the discretization system corresponding to system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) as follows: (15) $\begin{array}{rl}{x}_{k+1}& =x_k+x_k(a-b{x}_k-\frac{(\lambda +\alpha y_k)y_k}{1+h(\lambda +\alpha y_k)x_k})\mathrm{\Delta t}+{\sigma }_1{x}_k\sqrt{\mathrm{\Delta t}}{\xi }_k\\ & +\frac{{\sigma }_1^2}{2}{x}_k({\xi }_k^2-1)\mathrm{\Delta t},\\ y_{k+1}& =y_k+y_k(r-\frac{f{y}_k}{m+x_k})\mathrm{\Delta t}+{\sigma }_2{y}_k\sqrt{\mathrm{\Delta t}}{\eta }_k+\frac{{\sigma }_2^2}{2}{y}_k({\eta }_k^2-1)\mathrm{\Delta t},\end{array}$(15) where $\mathrm{\Delta t}$ is time increment, ${\xi }_k,{\eta }_k$ are two independent Gaussian random variables, and ${\xi }_k,{\eta }_k\sim N(0,1)$.

Example 6.1

For system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ), we choose $a=0.8,b=0.3,\lambda =0.5,\alpha =0.25,h=1,r=0.5,f=0.8,m=1$

By Theorem 3.3, we know that the solution of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) is stochastically ultimately bounded. In order to verify this conclusion, we draw the solution curve ‘$t-x\left(t\right),y\left(t\right)$’ with parameters (1) ${\sigma }_1=0,{\sigma }_2=0$, (2) ${\sigma }_1=0.15,{\sigma }_2=0.15$, and (3) ${\sigma }_1=0.3,{\sigma }_2=0.3$ respectively. In Figure 1, the figures of solution ‘$t-x\left(t\right),y\left(t\right)$’ with the initial value $(0.5,0.5)$ are drawn, which shows that the numerical result is identical to our theoretical result. Note that the solution of the deterministic system (4) (system (4) with ${\sigma }_1=0,{\sigma }_2=0$) is ultimately bounded since $\mathrm{d}x\left(t\right)\le x\left(t\right)[a-\mathrm{bx}(t\left)\right]\mathrm{d}t$ and $\mathrm{d}y\left(t\right)=y\left(t\right)[r-\mathrm{fy}(t\left)/\right(m+x\left(t\right)\left)\right]\mathrm{d}t$, which leads to $\underset{t\to +\mathrm{\infty }}{lim sup}x\left(t\right)\le a/b$ and $\underset{t\to +\mathrm{\infty }}{lim sup}y\left(t\right)\le r[m+a/b]/f$. Therefore, we conclude that noises do not influence the ultimate boundedness of the solution.

Figure 1. Solutions of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) with $a=0.8,b=0.3,\lambda =0.5,\alpha =0.25,h=1,r=0.5,f=0.8,m=1.$ Green lines: ${\sigma }_1={\sigma }_2=0$. Red lines: ${\sigma }_1={\sigma }_2=0.15$. Blue lines: ${\sigma }_1={\sigma }_2=0.3$.

Example 6.2

For system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ), we choose $a=1,b=0.6,\lambda =0.8,\alpha =0.25,h=1,r=1.5,f=1,m=\mathrm{0.2.}$

Considering system (4) with ${\sigma }_1={\sigma }_2=0.$ It is easily know that the solution $\left(x\right(t),y(t\left)\right)$ of the deterministic system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) with the initial value $\left(x\right(0),y(0\left)\right)\in R_{+}^2$ is positive. Therefore, it is true, for $r-\frac{1}{h}>0$, that $\begin{array}{rl}\mathrm{d}x\left(t\right)+\mathrm{d}y\left(t\right)& \ge \left[\mathrm{ax}\right(t)-b{x}^2(t)-\frac{1}{h}y(t)+\mathrm{ry}(t)-\frac{f}{m}{y}^2(t\left)\right]\mathrm{d}t\\ & \ge min\{a,r-\frac{1}{h}\}\left[x\right(t)+y(t\left)\right]-max\{b,\frac{f}{m}\}\left[x^2\right(t)+y^2(t\left)\right]\\ & \ge min\{a,r-\frac{1}{h}\}\left[x\right(t)+y(t\left)\right]-max\{b,\frac{f}{m}\}\left[x\right(t)+y(t){]}^2.\end{array}$Further, we have $\underset{t\to +\mathrm{\infty }}{lim inf}\left[x\right(t)+y(t\left)\right]\ge min\{a,r-\frac{1}{h}\}/max\{b,\frac{f}{m}\}$. Let $m_{x+y}=min\{a,r-\frac{1}{h}\}/max\{b,\frac{f}{m}\}$. Then we have $2\underset{t\to +\mathrm{\infty }}{lim inf}\Vert \left(x\right(t),y(t\left)\right){\Vert }^2\ge \underset{t\to +\mathrm{\infty }}{lim inf}\left[x\right(t)+y(t){]}^2=m_{x+y}^2$. Meanwhile, similar to Example 6.1, the ultimate boundedness of the solution for the deterministic system (4) is easily obtained. Hence, by Definition 4.1, we know that system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) with ${\sigma }_1={\sigma }_2=0$ is permanent.

Next, the white noises are considered. Let ${\sigma }_1=0.2,{\sigma }_2=0.2$. By computing, we have $\frac{1}{2}max\{{\sigma }_1^2,{\sigma }_2^2\}-min\{a,r-\frac{1}{h}\}=-0.48<0,$which implies that all the parameters satisfy ‘Assumption A’. Similarly, we choose ${\sigma }_1=0.5,{\sigma }_2=\mathrm{0.5.}$ It is obvious that ‘Assumption A’ is also true. By Theorem 4.3, system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) is stochastically permanent. In Figure 2, we draw figures of solution ‘$t-x\left(t\right),y\left(t\right)$’ with the initial value $(0.2,0.2)$, and the stochastical permanence of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) is shown.

Figure 2. Solutions of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) with $a=1,b=0.6,\lambda =0.8,\alpha =0.25,h=1,r=1.5,f=1,m=\mathrm{0.2.}$ Green lines: ${\sigma }_1={\sigma }_2=0$. Red lines: ${\sigma }_1={\sigma }_2=0.2$. Blue lines: ${\sigma }_1={\sigma }_2=0.5$.

It is revealed by Theorem 4.3 that white noises do not influence the permanence when the intensity of white noises is not too strong and satisfies ‘Assumption A’.

Example 6.3

For system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ), we choose $a=1,b=0.6,\lambda =0.8,\alpha =0.25,h=1,r=1.5,f=1,m=\mathrm{0.2.}$

These parameters are the same as those in Example 6.2. In order to see the influence of white noises, in the following, we change the intensity of white noises.

(i)	Let ${\sigma }_1=0.2$, ${\sigma }_2=\mathrm{1.8.}$ Obviously, ${\sigma }_1$ is also the same as that in Example 6.2 but ${\sigma }_2$ becomes bigger than that in Example 6.2. Then, we get $\frac{{\sigma }_1^2}{2}=0.02<a,\frac{{\sigma }_2^2}{2}=1.62>r.$It implies that all the parameters satisfy the condition (i) of Theorem 5.2. In Figure 3, we draw the curve of solution ‘$t-x\left(t\right),y\left(t\right)$’ with initial value $(0.2,0.2)$. Meanwhile, the curve of solution ‘$t-x\left(t\right),y\left(t\right)$’ with initial value $(0.2,0.2)$ for the corresponding deterministic system are drawn. It is shown by this figure that conclusion (i) of Theorem 5.2 is true. Since the predator is influenced by the white noise with the strong intensity, its dynamic behaviour changes. The strong white noise influences the dynamic behaviour of the population while the weak white noise influences are weak. while influences of the weak white noise are tiny.
(ii)	Let ${\sigma }_1=1.5$, ${\sigma }_2=\mathrm{0.2.}$ Obviously, ${\sigma }_2$ is also the same as that in Example 6.2 but ${\sigma }_1$ becomes bigger than that in Example 6.2. Then, we get $\frac{{\sigma }_1^2}{2}=1.125>a,\frac{{\sigma }_2^2}{2}=0.02<r.$It implies that all the parameters satisfies the condition (ii) of Theorem 5.2. In Figure 4, we draw the curve of solution ‘$t-x\left(t\right),y\left(t\right)$’ with the initial value $(0.2,0.2)$. Meanwhile, the curve of solution ‘$t-x\left(t\right),y\left(t\right)$’ with the initial value $(0.2,0.2)$ for the corresponding deterministic system is drawn. It is shown by this figure that conclusion (ii) of Theorem 5.2 is true. Since the prey is influenced by the white noise with the strong intensity, its dynamic behaviour changes. The strong white noise influences the dynamic behaviour of the population while the weak white noise influences is weak.
(iii)	let ${\sigma }_1=1.5$, ${\sigma }_2=\mathrm{1.8.}$ Obviously, ${\sigma }_1$ and ${\sigma }_1$ become bigger than that in Example 6.2. Then, we get $\frac{{\sigma }_1^2}{2}=1.125>a,\frac{{\sigma }_2^2}{2}=1.62>r.$It implies that all the parameters satisfies the condition (iii) of Theorem 5.2. In Figure 5, we draw the curve of solution ‘$t-x\left(t\right),y\left(t\right)$’ with initial value $(0.2,0.2)$. Meanwhile the curve of solution ‘$t-x\left(t\right),y\left(t\right)$’ with initial value $(0.2,0.2)$ for the corresponding deterministic system are drawn. It shows that conclusion (iii) of Theorem 5.2 is true. Moreover, white noise can change the properties of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) when the intensity of white noises is sufficiently large. The stochastic system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) will be extinct even if the corresponding deterministic system is persistent.

Figure 3. Solutions of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) with $a=1,b=0.6,\lambda =0.8,\alpha =0.25,h=1,r=1.5,f=1,m=\mathrm{0.2.}$ Red lines: ${\sigma }_1=0.2,{\sigma }_2=1.8$. Blue lines: ${\sigma }_1={\sigma }_2=0$.

Figure 4. Solutions of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) with $a=1,b=0.6,\lambda =0.8,\alpha =0.25,h=1,r=1.5,f=1,m=\mathrm{0.2.}$ Red lines: ${\sigma }_1=1.5,{\sigma }_2=0.2$. Blue lines: ${\sigma }_1={\sigma }_2=0$.

Figure 5. Solutions of system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) with $a=1,b=0.6,\lambda =0.8,\alpha =0.25,h=1,r=1.5,f=1,m=\mathrm{0.2.}$ Red lines: ${\sigma }_1=1.5,{\sigma }_2=1.8$. Blue lines: ${\sigma }_1={\sigma }_2=0$.

7. Conclusion

In this paper, a stochastic predator–prey system with modified Leslie–Gower and hunting cooperation among predators is studied. We first prove the existence and uniqueness of a global positive solution, and we investigate the stochastically ultimate boundedness, and the stochastic permanence of the solution. Then several sufficient conditions are established for persistence in mean and extinction. In the end, using Matlab software, some numerical simulations are provided to verify our results.

Some interesting conclusions are obtained. First, the ultimate boundedness of solution is not influenced by noises no matter how strong the intensity of noises is. However, white noises may influence the persistence. Noting that all parameters are the same except the values of ${\sigma }_1$ and ${\sigma }_2$ in Example 6.2 and 6.3; therefore, the influence of white noise is revealed directly by comparing Figures 2 and 3–5. Figure 2 reveals a fact that the permanence do not influenced by white noises when the intensity of white noises is not too strong and satisfies ‘Assumption A’. But, when the intensity of white noises is strong enough, what will occur for the dynamic behaviours of system (4)? The changes of the solution for system (4) are clear in Figures 3–5. If the predator is influenced by strong noise, its persistence will be destroyed, and it may be go to extinct(see Figure 3). It is the same for the prey (see Figure 4). Further, If the predator and the prey are influenced by strong noises, their persistence will be destroyed, and they may be go to extinct(see Figure 5). Therefore, in order to the persistence of system (4), system should avoid being affected by white noise.

In the future, we will study other properties further. For instance, we will investigate the existence of the stationary distribution for system (4(4) $\{\begin{array}{l}{\displaystyle \mathrm{d}x\left(t\right)=x\left(t\right)[a-\mathrm{bx}(t)-\frac{(\lambda +\mathrm{\alpha y}(t\left)\right)y\left(t\right)}{1+h(\lambda +\mathrm{\alpha y}(t\left)\right)x\left(t\right)}]\mathrm{d}t+{\sigma }_{1}x\left(t\right)\mathrm{d}{B}_1\left(t\right),}\\ {\displaystyle \mathrm{d}y\left(t\right)=y\left(t\right)(r-\frac{\mathrm{fy}\left(t\right)}{m+x\left(t\right)})\mathrm{d}t+{\sigma }_{2}y\left(t\right)\mathrm{d}{B}_2\left(t\right).}\end{array}$(4) ) and others.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The author(s) reported there is no funding associated with the work featured in this article.