A Stochastic Holling-Type II Predator-Prey Model with Stage Structure and Refuge for Prey

In this paper, we study a stochastic Holling-type II predator-prey model with stage structure and refuge for prey. Firstly, the existence and uniqueness of the global positive solution of the system are proved. Secondly, the stochastically ultimate boundedness of the solution is discussed. Next, su ﬃ cient conditions for the existence and uniqueness of ergodic stationary distribution of the positive solution are established by constructing a suitable stochastic Lyapunov function. Then, su ﬃ cient conditions for the extinction of predator population in two cases and that of prey population in one case are obtained. Finally, some numerical simulations are presented to verify our results.


Introduction
Population ecology is one of the most important research fields in biomathematics. There are three basic relationships among different species living in the same natural environment: (i) competitive relationship [1,2], for example, crops and weeds compete for fertilizer, sunlight, and other resources for their own growth; (ii) reciprocal relationship [3,4], such as rhinoceros birds feed on the parasites from rhinoceroses. For rhinoceroses, rhinoceros birds can get rid of the parasites and warn them when they are in danger; (iii) predation relationship [5,6], as the saying goes, "big fish eat small fish, small fish eat shrimps." Among the relationships of population interaction, the predator-prey system has been extensively studied. Lotka and Volterra put forward the predator-prey model in 1925 and 1926, respectively. Since then, many scholars have devoted themselves to the study of various predator-prey models and obtained abundant research results [7][8][9][10].
In all the above predator-prey systems, it is generally assumed that predators in the same population have identical predation ability, and the prey in the same population has identical viability and fertility. However, in the real world, the life of many creatures should be divided into two stages: immature and mature. There are recognizable morphological and behavioral differences that may exist between these stages. For example, the immature creatures have no fertility and predation ability, and their survival ability and defense ability are relatively weak. On the contrary, the mature creatures not only have reproductive ability and predator ability but also have strong survival ability and defensive ability. These will have a great impact on the dynamic behavior of the population. Thus, it is more practical to study the predator-prey model with stage structure. In view of this, many scholars have studied assorted predator-prey systems with stage structure and given corresponding dynamic analysis [11][12][13][14][15][16][17]. In particular, the following predator-prey model with stage structure for prey has been considered in [16].
where α, r 1 , β, η, m, b, r 2 , r, k, and η 1 are positive constants and _ x i = dx i /dt, i = 1, 2, 3:x 1 , x 2 and x 3 represent the size of the immature prey population, mature prey population, and predator population, respectively. And r 1 , r 2 and r indicate the mortality rates of immature prey population, mature prey population, and predator population, respectively. The parameter α denotes the birth rate of the immature prey population; β is the transformation rate of the mature prey population. η and η 1 denote the intraspecific competition coefficient of the immature prey population and predator population, respectively. k is the digestion constant. Predators only prey on immature prey. Functional response plays an important role in describing the predator-prey model, which refers to the response of predators' predation rate to the density of prey, that is, the predation effect of predators on prey. Model (1) uses bilinear functional response to describe the interaction between predator and prey. However, many researchers [18][19][20] have pointed out that nonlinear functional response can more accurately describe the trend of population density. Specifically, the growth rate of predators described by Holling-type II functional response is increasing with the increase of the number of prey [21]. When the survival space is far from saturated and the resources are sufficient, the growth rate of predators is relatively fast, but when the survival space and resources are limited, the growth rate changes less significantly and finally approaches a fixed value. Holling-type II functional response is according to some practical phenomena and takes into account the influence of density constraint. Therefore, it has been discussed by many scholars and extensively applied [13,[22][23][24][25][26]. In addition, it is well known that in many cases, prey can avoid predators through refuges. For example, crabs on the beach hide under sand or stones to avoid seabirds, thus increasing their survival chance. In recent years, different predator-prey systems with refuges have attracted much attention [6,[27][28][29][30]. In [27], Qi and Meng studied the threshold behavior of a stochastic predator-prey system with prey refuge and fear effect. They concluded that the survival rate of prey can be improved by increasing the strength of refuge. In [29], Ghosh et al. have studied a prey-predator model incorporating prey refuge and additional food for predators and found that slightly higher refuge is beneficial to the coexistence of species, but strong refuge will lead to the extinction of the predator population.
Inspired by the above motivations, we consider the following Holling-type II predator-prey model with stage struc-ture and refuge for prey: where m denotes the influence of predators on prey, a represents the half-saturation constant, b ∈ ½0, 1 denotes the refuge rate of immature prey, ð1 − bÞx 1 is the number of immature prey that predators can capture, and mð1 − bÞ is the capture rate of predators.
On the other hand, in nature, the biological population is inevitably affected by environmental noise. May [31] pointed out that due to the continuous fluctuation of the environment, the birth rate, mortality, environmental capacity, and other parameters in the model will show random fluctuations in varying degrees. However, the deterministic model does not consider the impact of environmental disturbance, so to some extent, it cannot accurately predict the dynamic behavior of the population. In order to describe the environmental noise better, many scholars consider the environmental noise as white noise, for instance, [22-27, 32, 33]. Considering the sensitivity of mortality to environmental noise, in this paper, we assume that the mortality rates r 1 , r 2 and r are disturbed by environmental noise, i.e., where σ i ði = 1, 2, 3Þ are the intensities of the white noise and B i ðtÞði = 1, 2, 3Þ represent independent standard Brownian motions which are defined on the complete probability space ðΩ, F, fFg t≥0 , ℙÞ with a filtration fFg t≥0 satisfying the usual conditions (that means it is increasing and right continuous while F 0 contains all ℙ-null sets). Then, the stochastic model corresponding to the deterministic model (2) is derived in the following form: Based on the analysis above, the paper focuses on the analysis of some dynamic behaviors of the stochastic Holling-type II predator-prey model with stage structure and refuge for prey. The major contributions are presented as follows. First, for the first time, the Holling-type II functional response is applied to establish a predator-prey model with stage structure and refuge for prey. Second, the existence and uniqueness as well as the stochastically ultimate 2 Advances in Mathematical Physics boundedness of positive solution of the stochastic model are analyzed. Third, some conditions for the existence and uniqueness of ergodic stationary distribution of the stochastic model are established. Finally, some extinction conditions of the stochastic model are studied. Moreover, in [26], Xu et al. studied a stochastic twopredator one-prey model with modified Leslie-Gower and Holling-type II schemes. The authors proved the existence and uniqueness of global positive solution of their stochastic model by using comparison theorem. In this paper, we consider it by constructing a Lyapunov function. Besides, in [4], Liu et al. discussed a mutualism system in random environments. They obtained the existence and uniqueness of a stable stationary distribution by means of Markov semigroup theory and Fokker-Planck equation. In this paper, we consider the existence and uniqueness of an ergodic stationary distribution of the stochastic model (4) by the stochastic Lyapunov function method.
The rest of this paper is arranged as follows. Some useful definitions and lemmas are given in Section 2. The existence and uniqueness of the global positive solution are discussed in Section 3. In Section 4, we obtain sufficient conditions for stochastically ultimate bounded of the prey and predator. In Section 5, the sufficient conditions for the existence of a unique ergodic stationary distribution of the positive solution are established. In Section 6, we establish sufficient conditions for the extinction of the predator and prey in two cases. The first case is the extinction of the predator, and another case is that both the prey and the predator are extinct. To verify our results, some numerical simulations are presented in Section 7. Finally, the conclusion is given in Section 8.

Preliminaries
In this section, some definitions and lemmas are given to prepare for further work.

Existence and Uniqueness of the Global Positive Solution
The existence and uniqueness of global positive solution are the premises of studying the properties of population dynamics. From [34], for any given initial data ðx 1 ð0Þ, x 2 ð0Þ, x 3 ð0ÞÞ ∈ ℝ 3 + , if the coefficients of SDE model (4) satisfy the linear growth condition and the local Lipschitz condition, then there exists a unique global positive solution ðx 1 ðtÞ, x 2 ðtÞ, x 3 ðtÞÞ ∈ ℝ 3 + . But the coefficients of model (4) only satisfy the local Lipschitz condition; then, the solution of the system will explode in a finite time. In this section, we will prove that SDE model (4) has a unique global positive solution.

Advances in Mathematical Physics
Proof. Since the coefficients of system (4) satisfy the local Lipschitz condition, then for any given initial value ðx 1 ð0Þ, x 2 ð0Þ, x 3 ð0ÞÞ ∈ ℝ 3 + , there is a unique local solution ðx 1 ðtÞ, x 2 ðtÞ, x 3 ðtÞÞ ∈ ℝ 3 + on t ∈ ð0, τ e Þ (see [34]), where τ e denotes the explosion time. To show this solution is global, we only need to prove τ e = ∞ a:s: Let n 0 > 0 be sufficiently large such that x i ð0Þ ∈ ½1/n 0 , n 0 ði = 1, 2, 3Þ. For any integer n ≥ n 0 , we define the following stopping time: We let inf ∅ = ∞ (∅ denotes the empty set). From the definition of stopping time, it is easy to know that τ n is increasing as n ⟶ ∞. Set τ ∞ = lim n⟶∞ τ n ; hence, τ ∞ ≤ τ e a:s: If we can show that τ ∞ = ∞ a:s:, then τ e = ∞ a:s: and ðx 1 ð0Þ, x 2 ð0Þ, x 3 ð0ÞÞ ∈ ℝ 3 + a:s: for all t ≥ 0. Thus, to complete the proof, we only need to prove τ ∞ = ∞ a:s: If this statement is false, for any ε ∈ ð0, 1Þ, there exists a constant T > 0 such that Therefore, there is an integer n 1 ≥ n 0 , for all n ≥ n 1 , we have Now, we define a C 2 -function V: where c is a positive constant that will be determined later. Since thus, V is a nonnegative function. Applying Itô formula to V, we have where LV : Choose c = ar/kmð1 − bÞ , then where K is a positive constant. The rest of the proof is similar to [37] and hence is omitted here. The proof is completed.

Stochastically Ultimate Boundedness
Theorem 5 shows that the solution of system (4) remains in the positive cone ℝ 3 + . However, this nonexplosion property in a population dynamical system is often not good enough. Hence, we need to consider another important asymptotic 4 Advances in Mathematical Physics property: stochastically ultimate boundedness, which means that the solution will be ultimately bounded with large probability.
Theorem 6. If the conditions σ 2 1 < 2r 1 + 2β − 1, σ 2 2 < 2r 2 − 1, and σ 2 3 < 2r − 2km − 1 hold, then, the solution of system (4) is stochastically ultimately bounded for any initial data By Itô formula, we have where Let According to the conditions, we can find that function f ðx 1 , x 2 , x 3 Þ has an upper bound. We assume that its upper bound is as follows: Let M = M 1 + 1 and noticing f ð0, 0, 0Þ = 0, thus, M > 0: According to (19), we can obtain By Itô formula, we have Integrating both sides of (24) from 0 to t and then taking expectations, we get Hence, we have namely, For any ε > 0, let H = ffiffiffiffi ffi M p / ffiffi ε p . By Chebyshev's inequality, we can obtain Then, The proof of Theorem 6 is completed.□

Existence of Ergodic Stationary Distribution
In this section, we will study the sufficient conditions for the existence and uniqueness of an ergodic stationary distribution of the positive solutions of system (4). The existence of an ergodic stationary distribution means that all species can coexist for a long time and are stochastic weakly persistent [17].

Advances in Mathematical Physics
Thus, Then, we obtain lim t⟶∞ x 3 ðtÞ = 0: Remark 9. Theorem 8 shows that the larger noise intensity or high mortality of predators as well as high refuge rate of immature prey can lead to the extinction of the predator population. However, under these conditions, we cannot judge whether the immature and mature prey populations are extinct or not.
Theorem 10. Assume that ðx 1 ðtÞ, x 2 ðtÞ, x 3 ðtÞÞ is the solution of system (4) with initial data ðx 1 ð0Þ, then, all the populations will die out exponentially with probability one, i.e., Proof. Let V = ln ðx 1 + x 2 Þ. By Itô formula, we have Rewriting as the following form Let the matrix it is clear that Thus, A is negative definite. Define λ max as the maximum eigenvalue of matrix A, then we have Thus, we obtain Integrating (72) from 0 to t and then dividing by t on both sides, we get By the strong law of large number of the martingale, we have According to the nonnegativity of x 1 and x 2 , we can obtain that ∃T 0 = T 0 ðωÞ and a set Ω ε ⊂ Ω such that ℙðΩ ε Þ > 1 − ε and By Itô formula, we have

Advances in Mathematical Physics
Integrating (77) from 0 to t and then dividing by t on both sides, we get Then, we obtain Let ε ⟶ 0, then lim t⟶∞ sup ððln x 3 ðtÞ − ln x 3 ð0ÞÞ/tÞ ≤ −r − ð1/2Þσ 2 3 < 0, which implies that lim t⟶∞ x 3 = 0.□ Remark 11. According to Theorem 10, we can find that the extinction of immature prey population leads to the extinction of predator population. This is true, since model (4) shows that the predator population has no additional food source.
Finally, we use the numerical simulation method to analyze the effect of prey refuge on system (4). The refuge rate b of immature prey varies from 0 to 1. b = 0 means that the immature prey population has no refuge, that is to say, the immature prey population is completely exposed to the predatory environment of the predator population. On the contrary, b = 1 indicates that the immature prey is totally protected by refuge. Now, we take b = 0 and b = 1, respectively, and other parameters are as (i). The numerical simulation results are shown in Figure 4. Figure 4 shows that when b = 0, the immature prey population, mature prey population, and predator population are persistent, that is, the predator population and the prey population coexist, which can be verified by Theorem 7. But when b = 1, the predator population becomes extinct, because the predator population has no additional food source. Theorem 8 also proves this result.

Conclusion
In this paper, a stochastic Holling-type II predator-prey model with stage structure and refuge for prey is analyzed. Firstly, the existence and uniqueness of the global positive solution for system (4) have been proved. Due to its nonexplosive property, we obtain the conditions for the stochastically ultimate boundedness of positive solution. Next, the sufficient conditions for the existence of a unique ergodic stationary distribution of the positive solution are established. Then, we obtain sufficient conditions for the extinction of predator population in two cases and that of prey population in one case. The theoretical results are proved by numerical simulations.
The following conclusions are verified by numerical simulations. Figure 1 shows that the existence of ergodic stationary distribution can allow all populations are coexistent and stochastic persistent in a long time. In Figure 2, it is shown that the large noise intensity will lead to the extinction of predator population, as shown in Figures 2(c) and 2(d), but the prey population can survive, as shown in Figures 2(a) and 2(b). Figure 3 shows that the extinction of immature prey population will lead to the extinction of predator population. This is due to the fact that predators do not have additional food sources. Finally, it can be seen from Figure 4(c) that if the immature prey population is completely protected by refuge, it will lead predator population towards extinction.
In addition, based on our results, from the perspective of ecological protection, we can provide additional food to the predators to ensure that the prey population and predator population maintain a high-density level. On the one hand, it can reduce the predation rate of the predator population to the prey [30]; on the other hand, it can ensure that the predator population will not be completely extinct when the refuge rate of the prey population is high.

Data Availability
Please contact the author for data requests.

Conflicts of Interest
The authors declare that they have no competing interests.