Extinction Analysis of Stochastic Predator–Prey System with Stage Structure and Crowley–Martin Functional Response

Abstract

In this paper, we researched some dynamical behaviors of a stochastic predator–prey system, which is considered under the combination of Crowley–Martin functional response and stage structure. First, we obtained the existence and uniqueness of the global positive solution of the system. Then, we studied the stochastically ultimate boundedness of the solution. Furthermore, we established two sufficient conditions, which are separately given to ensure the stochastic extinction of the prey and predator populations. In the end, we carried out the numerical simulations to explain some cases.

1. Introduction

Population dynamics is one of the main parts of biological mathematics. The predator–prey model is a classical problem in population research. Lotka and Volterra [1] researched the origin and theory of predator–prey model, which is given by

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\frac{dx}{dt}=a_{0}x(t)-c_{0}x(t)y(t),\hfill \\ \frac{dy}{dt}=e_0{c}_{0}x(t)y(t)-d_{0}y(t),\hfill \end{array}\right.\end{array}$$

(1)

where $x(t)$ and $y(t)$ represent the population density of the prey and predator, respectively; $a_0$ and $d_0$ denote the intrinsic growth rate and death rate, respectively; and $c_0$ and $e_0$ are the predation rate of a predator and nutrient-conversion rate, respectively.

An important feature of the predator–prey relationship is the functional response (i.e., the rate of prey consumption by an average predator). Mukherjee [2] discussed persistence and bifurcation on the predator–prey system of Holling Type II. Liu and Zhong [3] researched permanence and extinction for the delayed periodic predator–prey system with Holling Type II response function and diffusion. Zhang and Yang [4] studied Hopf bifurcation in the predator–prey system with Holling Type III functional response and time delays. The functional responses of Holling Types I–III are prey-dependent, which have been researched by many scholars. However, the functional response is inevitably influenced by the behavior of a predator, such as foraging and competing. Therefore, many scholars studied various types of predator-dependent functions. Gilliam and Skalski [5] claimed that the predator-dependent can provide better descriptions of predators feeding over a range of predator–prey abundances by comparing the statistical evidence from 19 predator–prey systems with the three predator-dependent functional responses (Hassell–Varley [6], Beddington–DeAngelis [7,8,9], and Crowley–Martin [10,11]), and, in some cases, the Crowley–Martin functional response is better. On the other hand, compared with the Hassell–Varley and Beddington–DeAngelis functional responses, the Crowley–Martin functional response is more suitable for the case that the predator feeding rate is decreased by higher predator density even when prey density is high. Thus, we consider the Crowley–Martin functional response in this paper.

In the classical predator–prey model, it is always assumed that each predator has the same predation capacity, and each prey has the same risk of predation. This assumption is unrealistic for many species. In nature, there are many species whose individuals have a life history that takes them through two stages, immature and mature. Individuals of different age groups exhibit different biological behaviors. In view of this, many scholars have studied the predator–prey system with stage structure [12,13,14,15,16]. Sun and Huo [17] considered bifurcation and stability in the predator–prey model with stage structure for the predator. Xu [18] discussed the global dynamics of the predator–prey model with time delay and stage structure for the prey. Lu [19] studied the stage-structured predator–prey model with predation over juvenile prey. However, few researchers have studied the predator–prey model with a stage structure for predator and prey. In nature, we know that immature predators have no predatory capacity. Meanwhile, many species hatch from eggs. For example, the Saltcedar leaf beetle is such a pest. In view of its eggshell, pathogens may not be effective against an immature pest. Based on this situation, it is reasonable to assume that immature prey does not run the risk of being preyed on. In terms of pest and disease control, the stage-structure model can better describe the dynamic behavior of some species. Therefore, in this paper, we mainly consider the predator–prey system with stage structure for both predator and prey.

We also consider the impact of environmental noise. Many scholars have studied various types of stochastic predator–prey systems with stage structure and functional-response functions [20,21,22]. Liu and Jiang [23] researched the dynamics of a stochastic predator–prey model with stage structure for predator and Holling Type II functional response. Chen and You [24] studied permanence, extinction, and periodic solution of the predator–prey system with a Beddington–DeAngelis functional response and stage structure for prey. Liu and Zhong [25] discussed the asymptotic properties of a stochastic predator–prey model with a Crowley–Martin functional response.

The main contributions of our work can be summarized as follows. The predator–prey model with random perturbation and Crowley–Martin functional response is established, which consider stage structure on both prey and predator. The existence and uniqueness of the global positive solution of the system is proved. Some sufficient conditions are given, which ensure the solutions of the system are stochastically ultimate boundedness. Then, sufficient conditions for the extinction of prey and predator are given, respectively. Finally, the conclusion is verified by numerical simulation results.

The paper is organized as follows. In Section 2, we give two prey–predator models with stage structure and a Crowley–Martin functional response. One is deterministic, and another is stochastic, which is discussed through the manuscript. In Section 3, we prove the existence and uniqueness of the global positive solution. In Section 4, we obtain sufficient conditions for stochastically ultimate boundedness of the prey and predator. In Section 5, we establish sufficient conditions for extinction of the predator and prey in two cases. The first case is the prey and predator extinction; another case is the predator extinction. In Section 6, numerical simulations illustrate the theoretical results. Section 7 gives the conclusions and future research directions.

2. Preliminaries

Before giving the main results, we first introduce some mathematical symbols and formulas in this paper. Throughout this paper, we define

$${\mathbb{R}}_{+}^q=\{x=(x_1,x_2,\cdots ,x_q)\in {\mathbb{R}}^q:x_i>0,1\le i\le q\}.$$

Consider the process of q-dimensional $It\widehat{o}$

$$\begin{array}{c}\hfill dX(t)=LV(X(t),t)dt+g(X(t),t)dB(t),\end{array}$$

(2)

where $B(t)=(B_1(t),B_2(t),\cdots ,B_q(t)$) denotes independent standard Brownian motions defined on a complete probability space $(\mathsf{\Omega },\mathcal{F},{\left\{{\mathcal{F}}_t\right\}}_{t\ge 0},P)$ with a filtration ${\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}$, and assume that the constant initial value $X_0\in {\mathbb{R}}^q$. Differential operator L of Formula (4) is given by

$$L=\frac{\partial }{\partial t}+\sum _{i=1}^q{f}_i(X(t),t)\frac{\partial }{\partial X_i}+\frac{1}{2}\sum _{i,j=1}^q{[g(X(t),t),g{(X(t),t)}^T]}_{ij}\frac{{\partial }^2}{\partial x_i\partial x_j}.$$

We denote a function $V(x(t),t)$ defined on $C^{2,1}({\mathbb{R}}^q,\mathbb{R})$. Applying L on $V(X(t),t)$, one has

$$LV(X(t),t)=V_t(X(t),t)+V_X(X(t),t)f(X(t),t)+\frac{1}{2}trace[g^T(X(t),t)V_{XX}(X(t),t)g(X(t),t)],$$

where $V_t=\frac{\partial V}{\partial t}$, $V_X=(\frac{\partial V}{\partial X_1},\frac{\partial V}{\partial X_2},\cdots ,\frac{\partial V}{\partial X_q})$, $V_{XX}={(\frac{{\partial }^{2}V}{\partial X_i\partial X_j})}_{d\times d}$. By $It\widehat{o}$ formula, we can obtain

$$dV(X(t),t)=LV(X(t),t)dt+V_X(X(t),t)g(X(t),t)dB(t).$$

If the system is autonomous, the definition of operator L and $It\widehat{o}$ formula discussed above can be found in Reference [26].

Based on the statement in Section 1, consider the following model:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\frac{d{x}_1(t)}{dt}=a{x}_2(t)-d_1{x}_1(t)-p{x}_1(t),\hfill \\ \frac{d{x}_2(t)}{dt}=p{x}_1(t)-d_2{x}_2(t)-b_1{x}_2^2(t)-\frac{c{x}_2(t)y_2(t)}{1+\alpha x_2(t)+\beta y_2(t)+\alpha \beta x_2(t)y_2(t)},\hfill \\ \frac{d{y}_1(t)}{dt}=\frac{ec{x}_2(t)y_2(t)}{1+\alpha x_2(t)+\beta y_2(t)+\alpha \beta x_2(t)y_2(t)}-d_3{y}_1(t)-h{y}_1(t),\hfill \\ \frac{d{y}_2(t)}{dt}=h{y}_1(t)-d_4{y}_2(t)-b_2{y}_2^2(t),\hfill \end{array}\right.\end{array}$$

(3)

where $x_1(t)$ and $x_2(t)$ denote the densities of immature and mature prey at time t, respectively; $y_1(t)$ and $y_2(t)$ represent the densities of immature and mature predators at time t, respectively; the parameters a, $d_1$, $d_2$, $d_3$, $d_4$, $b_1$, $b_2$, p, h, e and c are positive constants, a is the birth rate of immature prey, p and h indicate maturity rate of immature prey and immature predator, respectively; $b_1$ and $b_2$ express the competition rate between a mature prey population and mature predator population, respectively; $d_1$ and $d_2$ are the death rates of immature and mature prey, respectively; and $d_3$ and $d_4$ represent the death rates of immature and mature predators, respectively.

May [27] pointed out that due to continuous fluctuation in the environment, the birth rate, death rates, carrying capacity, competition coefficients, and all other parameters involved with the model exhibit random fluctuation. Thus, we consider environmental random disturbance as follows:

$$-d_1\to -d_1+{\sigma }_1\dot{B_1}(t),-d_2\to -d_2+{\sigma }_2\dot{B_2}(t),-d_3\to -d_3+{\sigma }_3\dot{B_3}(t),-d_4\to -d_4+{\sigma }_4\dot{B_4}(t),$$

where $B_i(t)(i=1,2,3,4)$ represent independent standard Brownian motions, and ${\sigma }_i(i=1,2,3,4)$ are the intensities of the environmental random disturbance. Then, we can obtain the following system:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}d{x}_1(t)=[a{x}_2(t)-d_1{x}_1(t)-p{x}_1(t)]dt+{\sigma }_1{x}_1(t)d{B}_1(t),\hfill \\ d{x}_2(t)=[p{x}_1(t)-d_2{x}_2(t)-b_1{x}_2^2(t)-\frac{c{x}_2(t)y_2(t)}{1+\alpha x_2(t)+\beta y_2(t)+\alpha \beta x_2(t)y_2(t)}]dt+{\sigma }_2{x}_2(t)d{B}_2(t),\hfill \\ d{y}_1(t)=[\frac{ec{x}_2(t)y_2(t)}{1+\alpha x_2(t)+\beta y_2(t)+\alpha \beta x_2(t)y_2(t)}-d_3{y}_1(t)-h{y}_1(t)]dt+{\sigma }_3{y}_1(t)d{B}_3(t),\hfill \\ d{y}_2(t)=[h{y}_1(t)-d_4{y}_2(t)-b_2{y}_2^2(t)]dt+{\sigma }_4{y}_2(t)d{B}_4(t).\hfill \end{array}\right.\end{array}$$

(4)

In this paper, we mainly research some population characteristics of System (3).

3. Existence and Uniqueness of Global Positive Solution

As we know, the density of population $x_1(t)$, $x_2(t)$, $y_1(t)$ and $y_2(t)$ should be positive. Therefore, we give the following theorem to ensure that the system has a unique positive solution.

Theorem 1.

For any given initial data $x_1(0)>0$, $x_2(0)>0$, $y_1(0)>0$ and $y_2(0)>0$, there is a unique solution $(x_1(t),x_2(t),y_1(t),y_2(t))$ to System (3), and the solution remains in ${\mathbb{R}}_{+}^4$ with probability 1.

Proof of Theorem 1.

Since System (3) satisfies the local Lipschitz continuous condition, there is a local unique solution $\{x_1(t),x_2(t),y_1(t),y_2(t)\}\in {\mathbb{R}}_{+}^4$ for any initial data $\{x_1(0),x_2(0),y_1(0),y_2(0)\}\in {\mathbb{R}}_{+}^4$ on $t\in [0,{\tau }_e)$ (with probability 1), where ${\tau }_e$ is the explosion time. To show that the solution is global, we only need to prove ${\tau }_e=\infty$ a.s. Give the following conditions for the initial value:

$$\frac{1}{l_0}<min\{x_1(0),x_2(0),y_1(0),y_2(0)\}\le max\{x_1(0),x_2(0),y_1(0),y_2(0)\}<l_0,$$

where $l_0$ is a sufficiently large number. For each integer $l\ge l_0$, define the stopping time

$${\tau }_l=inf\{t\in (0,{\tau }_e):x_1(t)\notin (\frac{1}{l},l)or{x}_2(t)\notin (\frac{1}{l},l)or{y}_1(t)\notin (\frac{1}{l},l)or{y}_2(t)\notin (\frac{1}{l},l)\},$$

where, in this paper, we set $inf\varnothing =\infty$. According to the definition of ${\tau }_l$, it is clear that ${\tau }_l$ increases as $l\to \infty$. Set ${\tau }_{\infty }=li{m}_{l\to \infty }{\tau }_l$, whence ${\tau }_{\infty }\le {\tau }_e$. That is to say, in order to prove the solution is global, it is sufficient to show that ${\tau }_{\infty }=\infty$ a.s. Then, we define a $C^2$-function V: ${\mathbb{R}}_{+}^4\to {\mathbb{R}}_{+}$ by

$$\begin{array}{c}\hfill V(x_1(t),x_2(t),y_1(t),y_2(t))=\sum _{i=1}^2[x_i(t)-1-ln{x}_i(t)]+\sum _{i=1}^2[y_i(t)-1-ln{y}_i(t)].\end{array}$$

(5)

The non-negative of this function can be seen from

$$u-1-lnu\ge 0,\forall u>0.$$

Let $T>0$, for $0\le t\le {\tau }_m\wedge T$. Applying It$\widehat{o}$’s formula to $V(x_1(t),x_2(t),y_1(t),y_2(t))$, we have

$$\begin{array}{c}\hfill d(V(x_1,x_2,y_1,y_2))=L(V(x_1,x_2,y_1,y_2))dt+\sum _{i=1}^4[(x_i-1){\sigma }_{i}d{B}_i(t)].\end{array}$$

(6)

According to the definition of operator L, we have

$$L(V(x_1,x_2,y_1,y_2))=[F_1(x_1,x_2,y_1,y_2)+F_2(x_1,x_2,y_1,y_2)+H_1(x_1,x_2,y_1,y_2)+H_2(x_1,x_2,y_1,y_2)]dt,$$

where $F_1(x_1,x_2,y_1,y_2)=(1-\frac{1}{x_1})(a{x}_2-d_1{x}_1-p{x}_1)+\frac{{\sigma }_1^2}{2}$, $F_2(x_1,x_2,y_1,y_2)=(1-\frac{1}{x_2})(p{x}_1-d_2{x}_2-b_1{x}_2^2-\frac{c{x}_2{y}_2}{1+\alpha x_2+\beta y_2+\alpha \beta x_2{y}_2})+\frac{{\sigma }_2^2}{2}$, $H_1(x_1,x_2,y_1,y_2)=(1-\frac{1}{y_1})(\frac{ec{x}_2{y}_2}{1+\alpha x_2+\beta y_2+\alpha \beta x_2{y}_2}-d_3{y}_1-h{y}_1)+\frac{{\sigma }_3^2}{2}$ and $H_2(x_1,x_2,y_1,y_2)=(1-\frac{1}{y_2})(h{y}_1-d_4{y}_2-b_2{y}_2^2)+\frac{{\sigma }_4^2}{2}$.

Then, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill LV(x_1,x_2,y_1,y_2)& \le -b_1{x}_2^2+(a+b_1)x_2-b_2{y}_2^2+b_2{y}_2+p+h+\frac{c}{\beta }+\frac{ec}{\alpha \beta }+\sum _{i=1}^4{d}_i+\sum _{i=1}^4\frac{{\sigma }_i^2}{2}\le K,\hfill \end{array}\end{array}$$

(7)

where $K=\frac{{(a+b_1)}^2}{4b_1}+\frac{b_2}{4}+p+h+\frac{c}{\beta }+\frac{ec}{\alpha \beta }+{\sum }_{i=1}^4{d}_i+{\sum }_{i=1}^4\frac{{\sigma }_i^2}{2}>0$. It can be obtained from Formulas (6) and (7) that

$$\begin{array}{c}\hfill d(V(x_1,x_2,y_1,y_2))\le Kdt+\sum _{i=1}^4[(x_i-1){\sigma }_{i}d{B}_i(t)].\end{array}$$

(8)

Integrating both sides of Formula (8) from 0 to ${\tau }_l\wedge T$, we have

$$\begin{array}{cc}\hfill V(x_1({\tau }_l\wedge T),x_2({\tau }_l\wedge T),& y_1({\tau }_l\wedge T),y_2({\tau }_l\wedge T))\le V(x_1(0),x_2(0),y_1(0),y_2(0))+KT\hfill \\ \hfill & +{\int }_0^{{\tau }_l\wedge T}{\sigma }_1(x_1-1)d{B}_1(t)+{\int }_0^{{\tau }_l\wedge T}{\sigma }_2(x_2-1)d{B}_2(t)\hfill \\ \hfill & +{\int }_0^{{\tau }_l\wedge T}{\sigma }_3(y_1-1)d{B}_3(t)+{\int }_0^{{\tau }_l\wedge T}{\sigma }_4(y_2-1)d{B}_4(t).\hfill \end{array}$$

(9)

Taking expectations at both sides of Formula (9), it is easy to obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill EV(x_1({\tau }_l\wedge T),x_2({\tau }_l\wedge T),y_1({\tau }_l\wedge T),y_2({\tau }_l\wedge T))& \le M,\hfill \end{array}\end{array}$$

(10)

where $M=V(x_1(0),x_2(0),y_1(0),y_2(0))+KT$. According to the definition of ${\tau }_l$, there is some i ($i=1,2$), such that $x_i({\tau }_l,\omega )$ and $y_i({\tau }_l,\omega )$ equal either $\frac{1}{l}$ or l. Then, $V(x_1({\tau }_l,\omega ),x_2({\tau }_l,\omega ),y_1({\tau }_l,\omega ),y_2({\tau }_l,\omega ))$ is no less than either

$$l-1-lnlor\frac{1}{l}-1-ln\frac{1}{l}.$$

Then, one has

$$V(x_1({\tau }_l,\omega ),x_2({\tau }_l,\omega ),y_1({\tau }_l,\omega ),y_2({\tau }_l,\omega ))\ge [(l-1-lnl)\wedge (\frac{1}{l}-1-ln\frac{1}{l})].$$

According to Formula (10), we can obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill M& \ge EV(x_1({\tau }_{l\wedge T}),x_2({\tau }_{l\wedge T}),y_1({\tau }_{l\wedge T}),y_2({\tau }_{l\wedge T}))\hfill \\ \hfill & \ge E[1_{{\tau }_l\le T}(\omega )V(x_1({\tau }_l),x_2({\tau }_l),y_1({\tau }_l),y_2({\tau }_l))]\hfill \\ \hfill & \ge P\{{\tau }_l\le T\}[(l-1-lnl)\wedge (\frac{1}{l}-1-ln\frac{1}{l})].\hfill \end{array}\end{array}$$

(11)

Letting $l\to \infty$, we have

$$\underset{l\to \infty }{lim}P\{{\tau }_l\le T\}=0.$$

Since $T>0$ is arbitrary, we have

$$P\{{\tau }_{\infty }<\infty \}=0.$$

Then,

$$P\{{\tau }_{\infty }=\infty \}=1.$$

The proof of Theorem 1 is completed. □

4. Stochastically Ultimate Boundedness

Theorem 1 shows that the solution of System (3) remains in the positive cone $R_{+}^4$. However, this nonexplosion property in a population dynamical system is often not good enough. Therefore, the property of ultimate boundedness is more desired. First, we give the definition of stochastically ultimate boundedness.

Definition 1

([28]). With respect to System (3), the solution is said to be stochastically ultimate bounded, if for $ϵ\in (0,1)$, there is a positive constant $H=H(ϵ)$ such that for any initial data $\{x_1(0),x_2(0),y_1(0),y_2(0)\}\in {\mathbb{R}}_{+}^4$, the solution $\{x_1(t),x_2(t),y_1(t),y_2(t)\}$ has the property that

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim\; sup}P\left\{\right|X(t)|\ge H\}\le ϵ,\end{array}$$

(12)

where $\left|X(t)\right|={(x_1^2+x_2^2+y_1^2+y_2^2)}^{\frac{1}{2}}$.

Assumption 1.

${\sigma }_1^2-2d_1-2p+1<0$, ${\sigma }_2^2-2d_2+1<0$, ${\sigma }_3^2-2d_3-2h+1<0$ and ${\sigma }_4^2-2d_4+1<0$.

Theorem 2.

Under Assumption 1, the solution of System (3) is stochastically ultimately bounded for any initial data $\{x_1(0),x_2(0),y_1(0),y_2(0)\}\in {\mathbb{R}}_{+}^4$.

Proof of Theorem 2.

For $\{x_1(t),x_2(t),y_1(t),y_2(t)\}\in {\mathbb{R}}_{+}^4$, define $V(x_1(t),x_2(t),y_1(t),y_2(t))$ as the following

$$V(x_1(t),x_2(t),y_1(t),y_2(t))=\sum _{i=1}^2{x}_i^2(t)+\sum _{i=1}^2{y}_i^2(t).$$

By $It\widehat{o}$’s formula, we have

$$\begin{array}{c}\hfill dV(x_1,x_2,y_1,y_2)=LV(x_1,x_2,y_1,y_2)dt+2{\sigma }_1{x}_1^{2}d{B}_1(t)+2{\sigma }_2{x}_2^{2}d{B}_2(t)+2{\sigma }_3{y}_1^{2}d{B}_3(t)+2{\sigma }_4{y}_2^{2}d{B}_4(t).\end{array}$$

(13)

Therefore, it is easy to derive

$$\begin{array}{cc}\hfill LV(x_1,x_2,y_1,y_2)& =-2b_1{x}_2^3-2b_2{y}_2^3+({\sigma }_1^2-2d_1-2p)x_1^2+({\sigma }_2^2-2d_2)x_2^2+({\sigma }_3^2-2d_3-2h)y_1^2\hfill \\ \hfill & +({\sigma }_4^2-2d_4)y_2^2+2(a+p)x_1{x}_2+2h{y}_1{y}_2+\frac{2ec{x}_2{y}_1{y}_2}{1+\alpha x_2+\beta y_2+\alpha \beta x_2{y}_2}\hfill \\ \hfill & -\frac{2c{x}_2^2{y}_2}{1+\alpha x_2+\beta y_2+\alpha \beta x_2{y}_2}\hfill \\ \hfill & \le ({\sigma }_1^2-2d_1-2p+1)x_1^2+({\sigma }_2^2-2d_2+1)x_2^2+({\sigma }_3^2-2d_3-2h+1)y_1^2\hfill \\ \hfill & +({\sigma }_4^2-2d_4+1)y_2^2+2(a+p)x_1{x}_2+2h{y}_1{y}_2+\frac{2ec{y}_1}{\alpha \beta }-x_1^2-x_2^2-y_1^2-y_2^2.\hfill \end{array}$$

(14)

Let

$$\begin{array}{cc}\hfill f(x_1,x_2,y_1,y_2)=& ({\sigma }_1^2-2d_1-2p+1)x_1^2+({\sigma }_2^2-2d_2+1)x_2^2+({\sigma }_3^2-2d_3-2h+1)y_1^2\hfill \\ \hfill & +({\sigma }_4^2-2d_4+1)y_2^2+2(a+p)x_1{x}_2+2h{y}_1{y}_2+\frac{2ec{y}_1}{\alpha \beta }.\hfill \end{array}$$

(15)

Under Assumption 1, it is easy to find that function $f(x_1,x_2,y_1,y_2)$ has an upper bound. We assume that its upper bound is as follows

$$\begin{array}{cc}\hfill M=\underset{(x_1,x_2,y_1,y_2)\in {\mathbb{R}}_{+}^4}{sup}& \left\{f(x_1,x_2,y_1,y_2)\right\}.\hfill \end{array}$$

(16)

Letting $N=M+1$ and noticing $f(0,0,0,0)=0$, we have $N>0$. According to Formula (14), we can obtain

$$\begin{array}{cc}\hfill dV(x_1,x_2,y_1,y_2)& \le [N-(x_1^2+x_2^2+y_1^2+y_2^2)]dt+2{\sigma }_1{x}_1^{2}d{B}_1(t)\hfill \\ \hfill & +2{\sigma }_2{x}_2^{2}d{B}_2(t)+2{\sigma }_3{y}_1^{2}d{B}_3(t)+2{\sigma }_4{y}_2^{2}d{B}_4(t).\hfill \end{array}$$

(17)

By $It\widehat{o}$’s formula, we have

$$\begin{array}{cc}\hfill d[e^{t}V(x_1,x_2,y_1,y_2)]& =e^{t}V(x_1,x_2,y_1,y_2)dt+e^{t}dV(x_1,x_2,y_1,y_2)\hfill \\ \hfill & \le N{e}^{t}dt+2{\sigma }_1{x}_1^{2}d{B}_1(t)+2{\sigma }_2{x}_2^{2}d{B}_2(t)+2{\sigma }_3{y}_1^{2}d{B}_3(t)+2{\sigma }_4{y}_2^{2}d{B}_4(t).\hfill \end{array}$$

(18)

Integrating both sides of Formula (18) from 0 to t and then taking expectations, we have

$$e^{t}E[V(x_1,x_2,y_1,y_2)]\le V(x_1(0),x_2(0),y_1(0),y_2(0))+N{e}^t-N.$$

Hence, we have

$$\underset{t\to \infty }{lim\; sup}E[V(X(t))]\le N,$$

where $X(t)=(x_1,x_2,y_1,y_2)$. Then, we have

$$\underset{t\to \infty }{lim\; sup}{E[|X(t)|}^2]\le N.$$

For any $ϵ>0$, let $H=\frac{\sqrt{N}}{\sqrt{ϵ}}$. By Chebyshev’s inequality, we can obtain

$$P\left\{\right|X(t)|>H\}\le \frac{E(|X(t){|}^2)}{H^2}.$$

Then,

$$\underset{t\to \infty }{lim\; sup}P\left\{\right|X(t)|>H\}\le \frac{N}{H^2}=ϵ.$$

The proof of Theorem 2 is completed. □

5. Stochastic Extinction

In this section, we show that the population becomes extinct with probability one.

Theorem 3.

Assume that $\{x_1(t),x_2(t),y_1(t),y_2(t)\}$ is the solution of System (3) with initial data $\{x_1(0),$ $x_2(0),y_1(0),y_2(0)\}$. Then,

(i) all prey and predators die out exponentially with probability one, if $(2d_1+{\sigma }_1^2)(2a-2d_2-{\sigma }_2^2)<{(a-d_1-d_2)}^2$;

(ii) predators $y_1(t)$ and $y_2(t)$ die out exponentially with probability one, if $(2d_3+{\sigma }_3^2)(\frac{2ec}{\alpha }-2d_4-{\sigma }_4^2)<{(\frac{ec}{\alpha }-d_3-d_4)}^2$.

Proof of Theorem 3.

According to System (3), we get

$$d(x_1+x_2)=d{x}_1+d{x}_2=L(x_1+x_2)dt+{\sigma }_1{x}_{1}d{B}_1(t)+{\sigma }_2{x}_{2}d{B}_2(t),$$

where $L(x_1+x_2)=a{x}_2-d_1{x}_1-b_1{x}_2^2-d_2{x}_2-\frac{c{x}_2{y}_2}{1+\alpha x_2+\beta y_2+\alpha \beta x_2{y}_2}$. Let $V(x_1,x_2)=ln(x_1+x_2)$. By $It\widehat{o}$’s formula, we can obtain

$$\begin{array}{cc}\hfill dV(x_1,x_2)& =\frac{1}{x_1+x_2}(d{x}_1+d{x}_2)-\frac{1}{2{(x_1+x_2)}^2}[{(d{x}_1)}^2+{(d{x}_2)}^2],\hfill \\ \hfill & =[\frac{1}{x_1+x_2}(a{x}_2-d_1{x}_1-d_2{x}_2-b_1{x}_2^2-\frac{c{x}_2{y}_2}{1+\alpha x_2+\beta y_2+\alpha \beta x_2{y}_2})\hfill \\ \hfill & -\frac{{\sigma }_1^2{x}_1^2}{2{(x_1+x_2)}^2}-\frac{{\sigma }_2^2{x}_2^2}{2{(x_1+x_2)}^2}]dt+\frac{{\sigma }_1{x}_1}{x_1+x_2}d{B}_1(t)+\frac{{\sigma }_2{x}_2}{x_1+x_2}d{B}_2(t).\hfill \end{array}$$

(19)

Then,

$$\begin{array}{cc}\hfill LV(x_1,x_2)& =\frac{1}{2{(x_1+x_2)}^2}[2(x_1+x_2)(a{x}_2-d_1{x}_1-d_2{x}_2-b_1{x}_2^2-\frac{c{x}_2{y}_2}{1+\alpha x_2+\beta y_2+\alpha \beta x_2{y}_2})\hfill \\ \hfill & -{\sigma }_1^2{x}_1^2-{\sigma }_2^2{x}_2^2],\hfill \\ \hfill & \le \frac{1}{2{(x_1+x_2)}^2}[2(x_1+x_2)(a{x}_2-d_1{x}_1-d_2{x}_2)-{\sigma }_1^2{x}_1^2-{\sigma }_2^2{x}_2^2].\hfill \end{array}$$

(20)

We can write term

$$2(x_1+x_2)(a{x}_2-d_1{x}_1-d_2{x}_2)-{\sigma }_1^2{x}_1^2-{\sigma }_2^2{x}_2^2,$$

in the following way:

$$(x_1(t),x_2(t))\left[\begin{array}{cc}-2d_1-{\sigma }_1^2& a-d_1-d_2\\ a-d_1-d_2& 2a-2d_2-{\sigma }_2^2\end{array}\right]{(x_1(t),x_2(t))}^T.$$

Letting the matrix

$$A_1=\left[\begin{array}{cc}-2d_1-{\sigma }_1^2& a-d_1-d_2\\ a-d_1-d_2& 2a-2d_2-{\sigma }_2^2\end{array}\right],$$

it is clear that matrix $A_1$ is negative definite under the condition in (i). Define ${\lambda }_{max}$ as the maximum eigenvalue of matrix $A_1$. According to the condition in (i), we obtain

$$(x_1(t),x_2(t))\left[\begin{array}{cc}-2d_1-{\sigma }_1^2& a-d_1-d_2\\ a-d_1-d_2& 2a-2d_2-{\sigma }_2^2\end{array}\right]{(x_1(t),x_2(t))}^T\le -\left|{\lambda }_{max}\right|(x_1^2(t)+x_2^2(t)).$$

Therefore, we have

$$\begin{array}{cc}\hfill dV(x_1,x_2)=d(ln(x_1+x_2))& \le [\frac{-|{\lambda }_{max}|(x_1^2+x_2^2)}{2{(x_1+x_2)}^2}]dt+\frac{{\sigma }_1{x}_1}{x_1+x_2}d{B}_1(t)+\frac{{\sigma }_2{x}_2}{x_1+x_2}d{B}_2(t),\hfill \\ \hfill & \le -\frac{1}{4}\left|{\lambda }_{max}\right|dt+\frac{{\sigma }_1{x}_1}{x_1+x_2}d{B}_1(t)+\frac{{\sigma }_2{x}_2}{x_1+x_2}d{B}_2(t).\hfill \end{array}$$

(21)

Integrating both sides of Formula (19), it can be obtained that

$$\begin{array}{c}\hfill ln(x_1+x_2)\le ln(x_1(0)+x_2(0))-\frac{1}{4}\left|{\lambda }_{max}\right|t+{\int }_0^t\frac{{\sigma }_1{x}_1}{x_1+x_2}d{B}_1(t)+{\int }_0^t\frac{{\sigma }_2{x}_2}{x_1+x_2}d{B}_2(t).\end{array}$$

(22)

Let $Z_1(t)={\int }_0^t\frac{{\sigma }_1{x}_1}{x_1+x_2}d{B}_1(t)$ and $Z_2(t)={\int }_0^t\frac{{\sigma }_2{x}_2}{x_1+x_2}d{B}_2(t)$, where $Z_1(t)$ and $Z_2(t)$ are local martingales. By the strong law of large numbers for local martingales (see, e.g., Reference [26]), we can obtain the following properties:

$$\underset{t\to \infty }{lim}\frac{Z_i(t)}{t}=0a.s.i=1,2.$$

Therefore, we can get

$$\underset{t\to \infty }{lim\; sup}\frac{ln(x_1+x_2)}{t}\le -\frac{1}{4}\left|{\lambda }_{max}\right|<0a.s.$$

Then, we have

$$\underset{t\to \infty }{lim}{x}_1(t)\to 0,a.s.\underset{t\to \infty }{lim}{x}_2(t)\to 0.a.s.$$

Then, with the extinction of the prey, we find the predator dies out according to System (3). The discussion of the predator population is similar. We have

$$\begin{array}{cc}\hfill L(ln(y_1+y_2))& \le \frac{1}{2{(y_1+y_2)}^2}[(y_1(t),y_2(t))A_2{(y_1(t),y_2(t))}^T],\hfill \end{array}$$

(23)

where

$$A_2=\left[\begin{array}{cc}-2d_3-{\sigma }_3^2& \frac{ec}{\alpha }-d_3-d_4\\ \frac{ec}{\alpha }-d_3-d_4& \frac{2ec}{\alpha }-2d_4-{\sigma }_4^2\end{array}\right].$$

Let ${\overline{\lambda }}_{max}$ be the maximum eigenvalue of matrix $A_2$. Under the condition in (ii), we have

$$\begin{array}{c}\hfill ln(y_1+y_2)\le ln(y_1(0)+y_2(0))-\frac{1}{4}\left|{\overline{\lambda }}_{max}\right|t+{\int }_0^t\frac{{\sigma }_3{y}_1}{y_1+y_2}d{B}_3(t)+{\int }_0^t\frac{{\sigma }_4{y}_2}{y_1+y_2}d{B}_4(t).\end{array}$$

(24)

Then, it can be obtained that

$$\underset{t\to \infty }{lim\; sup}\frac{ln(y_1+y_2)}{t}\le -\frac{1}{4}\left|{\overline{\lambda }}_{max}\right|<0a.s.$$

which implies

$$\underset{t\to \infty }{lim}{y}_1(t)\to 0,a.s.\underset{t\to \infty }{lim}{y}_2(t)\to 0.a.s.$$

The proof of Theorem 3 is completed. □

6. Numerical Simulations

In this section, we illustrate our theoretical results using the numerical simulations of System (3). We randomly selected the initial condition in $(0,1)$. The initial state of the system is $(0.6324,0.8147,0.127,0.2785)$. We used the Milstein method mentioned in Reference [29] to substantiate the analytical findings. Consider the discretization transformation of System (3):

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{x}_1^{j+1}=x_1^j+(a{x}_2^j-d_1{x}_1^j-f{x}_1^j)\Delta t+{\sigma }_1{x}_1^j\sqrt{\Delta t}{ϵ}_{1,j}+\frac{1}{2}{\sigma }_1^2{x}_1^j({ϵ}_{1,j}^2\Delta t-\Delta t),\hfill \\ x_2^{j+1}=x_2^j+(f{x}_1^j-d_2{x}_2^j-b_1{x}_2^{2j}-\frac{c{x}_2^j{y}_2^j}{1+\alpha x_2^j+\beta y_2^j+\alpha \beta x_2^j{y}_2^j})\Delta t+{\sigma }_2{x}_2^j\sqrt{\Delta t}{ϵ}_{2,j}+\frac{1}{2}{\sigma }_2^2{x}_2^j({ϵ}_{2,j}^2\Delta t-\Delta t),\hfill \\ y_1^{j+1}=y_1^j+(\frac{ec{x}_2^j{y}_2^j}{1+\alpha x_2^j+\beta y_2^j+\alpha \beta x_2^j{y}_2^j}-d_3{y}_1^j-h{y}_1^j)\Delta t+{\sigma }_3{y}_1^j\sqrt{\Delta t}{ϵ}_{3,j}+\frac{1}{2}{\sigma }_3^2{y}_1^j({ϵ}_{3,j}^2\Delta t-\Delta t),\hfill \\ y_2^{j+1}=y_2^j+(h{y}_1^j-d_4{y}_2^j-b_2{y}_2^{2j})\Delta t+{\sigma }_4{y}_2^j\sqrt{\Delta t}{ϵ}_{4,j}+\frac{1}{2}{\sigma }_4^2{y}_2^j({ϵ}_{4,j}^2\Delta t-\Delta t),\hfill \end{array}\right.\end{array}$$

(25)

where time increment $\Delta t$ is positive and ${ϵ}_{i,j}$($i=1,2,3,4$) are the Gaussian random variables that follow distribution $N(0,1)$. For System (3), the parameters are selected as follows:

(i). $a=0.7$, $b_1=0.9$, $b_2=0.9$, $c=0.8$, $e=\frac{7}{8}$, $p=0.3$, $h=0.5$, $\alpha =0.8$, $\beta =0.5$, $d_1=0.9$, $d_2=0.5$, $d_3=0.5$, $d_4=0.5$, ${\sigma }_1^2=0.5$, ${\sigma }_2^2=1$, ${\sigma }_3^2=0.1$, ${\sigma }_4^2=0.8$.

(ii). $a=0.3$, $b_1=0.2$, $b_2=0.8$, $c=0.7$, $e=0.7$, $p=0.6$, $h=0.5$, $\alpha =0.8$, $\beta =0.6$, $d_1=0.2$, $d_2=0.1$, $d_3=0.7$, $d_4=0.7$, ${\sigma }_1^2=0.01$, ${\sigma }_2^2=0.01$, ${\sigma }_3^2=0.2$, ${\sigma }_4^2=0.5$.

It is easy to verify that Parameters (i) satisfy the condition of the extinction of the prey population in Theorem 3. The corresponding numerical results are shown in Figure 1.

As can be clearly seen from Figure 1, $x_1(t)$, $x_2(t)$, $y_1(t)$ and $y_2(t)$ tend to zero in both the deterministic and stochastic models. Under Parameter (i), we have $(2d_1+{\sigma }_1^2)(2a-2d_2-{\sigma }_2^2)<{(a-d_1-d_2)}^2$. By Theorem 3, $x_1(t)$, $x_2(t)$, $y_1(t)$ and $y_2(t)$ tend to become extinct. Numerical simulations clearly support this result (see Figure 1). Therefore, Figure 1 provides evidence for the accuracy of Conclusion (i) in Theorem 3. Then, Under Parameter (ii), the corresponding numerical simulation results are as follows.

As can be clearly seen from Figure 2, $y_1(t)$ and $y_2(t)$ tend to zero in both the deterministic and stochastic models. By calculation, we can find that Parameter (ii) satisfies condition $(2d_3+{\sigma }_3^2)(\frac{2ec}{\alpha }-2d_4-{\sigma }_4^2)<{(\frac{ec}{\alpha }-d_3-d_4)}^2$. According to Theorem 3, $y_1(t)$ and $y_2(t)$ tend to become extinct. Numerical simulations clearly support this result (see Figure 2). Meanwhile, under Parameter (ii), we give the trajectories of $x_1(t)$ and $x_2(t)$ over a long period of time (see Figure 3). Figure 3 shows that the immature prey and mature prey are permanence for a long time.

Under the condition of Parameter (ii), according to Figure 2 and Figure 3, the predator tends to become extinct and the prey survives for a long time. In nature, this situation is reasonable.

7. Conclusions

In this paper, we researched the predator–prey system with a Crowley–Martin functional response function and environmental noise. In Reference [5], we found that the predator-dependent functional response is more reasonable than the prey-dependent functional response. In particular, the Crowley–Martin functional response is more suitable for the case that the predator feeding rate is decreased by higher predator density. Compared with Holling Types I–III functional responses, the Crowley–Martin functional response has more complex forms. From an analysis point of view, the theoretical analysis of predator–prey system with a Crowley–Martin functional response is more difficult, and the results are more complex. Meanwhile, we know that the system is inevitably affected by environmental noise. Therefore, we researched the predator–prey model with a Crowley–Martin functional response function and environmental noise. On this basis, we first attempted to consider the stage structure on both prey and predator. First, we proved the existence and uniqueness of the global positive solution of System (3). Next, we pointed out that the positive solution is stochastically bounded. Then, we gave sufficient conditions for the extinction of the predator and prey populations in two cases. Some interesting questions deserve further investigation; we will research the stability and stationary distribution of System (3) (see Reference [30]), and consider the impact of sudden changes and time delays on population characteristics (see Reference [31]) in the future. In addition, we will research the chaotic behavior of a predator–prey system and the Allee effect (see References [32,33]).