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

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.


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 dx dt = a 0 x(t) − c 0 x(t)y(t), dy dt = e 0 c 0 x(t)y(t) − d 0 y(t), (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 Consider the process of q-dimensional Itô dX(t) = LV(X(t), t)dt + g(X(t), t)dB(t), (2) where B(t) = (B 1 (t), B 2 (t), · · · , B q (t)) denotes independent standard Brownian motions defined on a complete probability space (Ω, F , {F t } t≥0 , P) with a filtration {F t } t≥0 , and assume that the constant initial value X 0 ∈ R q . Differential operator L of Formula (4) is given by We denote a function V(x(t), t) defined on C 2,1 (R q , R). Applying L on V(X(t), t), one has If the system is autonomous, the definition of operator L and Itô formula discussed above can be found in Reference [26]. Based on the statement in Section 1, consider the following model: 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: where B i (t)(i = 1, 2, 3, 4) represent independent standard Brownian motions, and σ i (i = 1, 2, 3, 4) are the intensities of the environmental random disturbance. Then, we can obtain the following system: In this paper, we mainly research some population characteristics of System (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.
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)} ∈ R 4 + for any initial data {x 1 (0), x 2 (0), y 1 (0), y 2 (0)} ∈ R 4 + on t ∈ [0, τ e ) (with probability 1), where τ e is the explosion time. To show that the solution is global, we only need to prove τ e = ∞ a.s. Give the following conditions for the initial value: where l 0 is a sufficiently large number. For each integer l ≥ l 0 , define the stopping time where, in this paper, we set in f ∅ = ∞. According to the definition of τ l , it is clear that τ l increases as That is to say, in order to prove the solution is global, it is sufficient to show that τ ∞ = ∞ a.s. Then, we define a C 2 -function V: The non-negative of this function can be seen from According to the definition of operator L, we have where K = (a+b 1 ) 2 4b 1 It can be obtained from Formulas (6) and (7) that Integrating both sides of Formula (8) from 0 to τ l ∧ T, we have Taking expectations at both sides of Formula (9), it is easy to obtain where Then, one has According to Formula (10), we can obtain Letting l → ∞, we have lim l→∞ P{τ l ≤ T} = 0.
Since T > 0 is arbitrary, we have Then, The proof of Theorem 1 is completed.

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. (3), the solution is said to be stochastically ultimate bounded, if for ∈ (0, 1), there is a positive constant H = H( ) such that for any initial data {x 1 (0),

Definition 1 ([28]). With respect to System
where (3) is stochastically ultimately bounded for any initial

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

Letting the matrix
it is clear that matrix A 1 is negative definite under the condition in (i). Define λ max as the maximum eigenvalue of matrix A 1 . According to the condition in (i), we obtain Therefore, we have Integrating both sides of Formula (19), it can be obtained that Let 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: Therefore, we can get lim sup Then, we have lim 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 L(ln(y 1 + y 2 )) ≤ 1 2(y 1 + y 2 ) 2 [(y 1 (t), y 2 (t))A 2 (y 1 (t), y 2 (t)) T ], (23) where Letλ max be the maximum eigenvalue of matrix A 2 . Under the condition in (ii), we have Then, it can be obtained that lim sup The proof of Theorem 3 is completed.

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): where time increment ∆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 = 7 8 , p = 0.3, h = 0.5, α = 0.8, β = 0.5, d 1 = 0.9, d 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 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 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 Figures 2 and 3, the predator tends to become extinct and the prey survives for a long time. In nature, this situation is reasonable.

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]).
Author Contributions: In this paper, C.X. was in charge of theoretical deduction, numerical simulation, and writing the paper. Y.Y. and G.R. were responsible for verifying the conclusion of the paper.