Hopf Bifurcation and Control of a Fractional-Order Delay Stage Structure Prey-Predator Model with Two Fear Effects and Prey Refuge

: A generalized delay stage structure prey-predator model with fear effect and prey refuge is considered in this paper via introducing fractional-order and fear effect induced by immature predators. Hopf bifurcation and control of this system are investigated though regarding the delay as the parameter. Firstly, by using the method of linearization and Laplace transform, the roots of the characteristic equation of the linearized system of the original system are discussed, and the sufﬁcient conditions for the system exhibits an unstable state of symmetrical periodic oscillation (Hopf bifurcation) are explored. Secondly, a linear delay feedback controller is added to the system to increase the stability domain successfully. Thirdly, numerical simulations are performed to validate the theoretical analysis, and the various impacts on the dynamical behavior of the system occurring by fear effects, prey refuge, and each fractional-order are illustrated, respectively. Furthermore, the inﬂuence of feedback gain on the bifurcation critical point is analyzed. Finally, an analysis based on the results and in-depth research about this system under the biological background is stated in the conclusion.


Introduction
The series of prey-predator models are paid widespread attention because they reflect the ecological phenomenon existing in the real world generally. The dynamic behavior of those models are investigated in-depth and a large number of valuable results have been obtained in the past few decades [1][2][3][4]. To keep the ecological balance, it is necessary to increase the survival rate of the prey in some ecosystems. The prey refuge is a suitable method to protect the prey population [5]. Prey-predator models with prey refuge are brought into focus and many worthy results are obtained [6][7][8].
Generally, the growth of many species are divided into two stages, immaturity and maturity, and the characteristics and behaviors of the different stage are quite distinguishing. For example, the predatory ability of mature predators is stronger than immature predators. In the past few years, researchers found that the prey-predator model with multi-stage structure is more reasonable than the one-stage model for describing the relationship between species [9][10][11].
Fear effect is another factor that impacts the dynamic behaviors of prey-predator models besides stage structure. For instance, the fear to predators could affect the birth rate of the prey, thereby affecting the population density of the prey [12]. Investigations have shown that the fear induced by the predator has an even greater effect on the prey than the direct killing [13][14][15]. Biologists discovered that many prey species have the inborn ability to identify predators in addition to acquired learning [16]. This means the fear effect may come from not only mature predators but also immature.
The organization of this paper is as follows. In Section 2, some definitions of fractional calculus and some basic knowledge are presented. In Section 3, the mathematical model is generalized, and the existence of the coexistence equilibrium point of the model is analyzed. In Section 4, Hopf bifurcation of the generalized system and the control of bifurcation of the controlled system are explored, respectively. In Section 5, numerical simulations are performed to further illustrate our theoretical results, and the influences of two fear effects, prey refuge, and fractional-order to the bifurcation of the system are given. Finally, a necessary conclusion explains the results and in-depth research about this system under the biological background.

Preliminary Knowledge
In this section, some basic definitions about fractional-order calculus and Hopf bifurcation used in the following sections are given.

Model Description
Xiao et al. [41] studied the following Beddington-DeAngelis prey-predator model with stage structure and prey refuge where x(t), y 1 (t), and y 2 (t) represent the population densities of prey, immature predator, and mature predator at time t, respectively. r is the birth rate of prey. d 1 and d 2 represent the natural mortality of immature predators and mature predators, respectively. c is the intraspecific competition rate of the prey. m ∈ [0, 1) is the prey refuge rate, and n represents the proportion of immature predators that grow into mature predators. a and b refer to the processing time of the mature predator and the strength of the interaction. α and β refer to the capture rate of the prey and the conversion rate of nutrients into the production of predator species, respectively. τ represents the delay due to the gestation of the mature predator. The authors investigated the local stability of the equilibrium point of the system and the influence of prey refuge on the densities of predator species and prey species.
In 2021, Wang and Hu [35] improved this model and discussed a Crowley-Martin prey-predator model with fear effect and prey refuge as follows where k is the prey's fear factor induced by mature predators and d 0 represents the natural mortality of prey. The existence and stability of the equilibrium point of the system Equation (8) have been established in [35].
Considering that the evolution of prey-predators system related to both the current moment and some past state of the species, and each species has different degree of dependence on the past, incommensurate fractional-orders are added to the system (8). Otherwise, inspired by the reference [16], the prey's fear effect is thought about not only mature predator but also immature predator. Thus, the system Equation (8) is generalized as the follows: where q i ∈ (0, 1] (i = 1, 2, 3) is fractional-order. k 1 and k 2 are the prey's fear factors induced by immature predators and mature predators, respectively. Obviously, the system has a zero equilibrium point E 0 = (0, 0, 0). When r > d 0 , the system has a predator-extinction equilibrium point E 1 = ( r−d 0 c , 0, 0).In fact, we are interested in the stability and stability switch at the coexistence equilibrium point of the system (9). Thus, it is necessary to find the conditions in that system (9) has a positive value equilibrium point. Lemma 1. When the following four conditions are satisfied, the system (9) has a unique coexistence equilibrium point E * = (x * , y * 1 , Proof. In fact, system (9) exists an coexistence equilibrium point E * = (x * , y * 1 , y * 2 ) means the following equation set has positive solution It is easy to obtain y 1 = d 2 n y 2 , and substituting y 1 into the first and second equations of Equation (10), one has if curve F(x, y 2 ) = 0 intersects curve G(x, y 2 ) = 0 in the first quadrant, then system (10) has positive solution. According to the first equation of Equation (12), one has In the first quadrant, if From F(x, y 2 ) = 0, one can get Substituting Equation (15) into inequation Equation (14), one has where For all x > 0, then the inequation Equation (16) is established, and it means dy 2 dx < 0. If y 2 = 0, then x (1) = r−d 0 c , and when x (1) > 0, there is r > d 0 . where Equation (17) has only one positive root y 2 (1) . According to second equation of Equation (12), in the first quadrant, one has If it is easy to obtain 0 < x (2) < x (1) . Therefore, if (H1)−(H4) are satisfied, F(x, y 2 ) and G(x, y 2 ) have a unique intersection (x * , y 2 * ) in the first quadrant, and then y 1 * > 0 can be obtained by the third equation of Equation (10).

Hopf Bifurcation Analysis and Control of System (9)
We are interested in the dynamical properties at the coexistence equilibrium point (x * , y 1 * , y 2 * ) of system (9). In this section, the Hopf bifurcation and control are analyzed in details.

Hopf Bifurcation Analysis of System (9)
Using the transformation * , w(t) = y 2 (t) − y 2 * , and linearizing the converted system, we can have where According to Definition 4, we analyze the conditions of emergence of Hopf bifurcation for system (9) at the coexistence equilibrium point.
In order to find the critical value of delay that the stability of system (19) switches, one can assume (23), substituting it into Equation (23) and separating the real and imaginary parts, we have where Solving Equation (24), it is easy to obtain the following result, By Equation (23), we can get It is easy to know that According to assumption (A1), the equation |U 1 (iω)| = |U 2 (iω)| has at least one positive root.
Combining with the formula sin 2 ωτ + cos 2 ωτ = 1, the value of ω can be solved. Without loss of generality, we assume that all positive roots are ω k (k = 1, 2, ..., K). By substituting each ω k into Equation (25) and the corresponding critical value of τ k can be obtained (for the exact mathematical expressions, please refer to [40]). In relation to the actual meaning of delay, we only pay attention to the value of τ when Hopf bifurcation occurs firstly, so the bifurcation critical value of delay is the critical value of frequency corresponding to τ 0 is denoted as ω 0 . According to Definition 4, we need to verify the transversality condition at the critical point (τ 0 , ω 0 ). Thus, it is necessary to give the following hypothesis Lemma 3. If the hypothesis (H8) holds, let s(τ) = γ(τ) + iω(τ) be the root of Equation (23) near τ = τ j satisfying γ(τ j ) = 0, ω(τ j ) = ω 0 , then the following transversality conditions established Proof. According to the implicit function derivation rule, deriving τ on both sides of Equation (23) respectively, one gets where It can be deduced from Equation (28) that where A 1 , A 2 are the real and imaginary parts of A(s), B 1 , B 2 are the real and imaginary parts of B(s). In terms of (H8), one has The proof of Lemma 3 is finished.

Hopf Bifurcation Control of System (9)
In this section, we focus on the control of Hopf bifurcation of system (9). From an ecological point of view, it is more effective to control the stability of the system by regulating the population density of mature predators than by regulating the population density of immature predators, as the mature predators play a dominant role in the ecosystem. A linear delay feedback controller L[y 2 (t) − y 2 (t − τ)] is added to the third equation of system (9) to control the emergence of Hopf bifurcation, i.e., the stability domain is regulated by controlling the population density of mature predators. The controlled system is and doing the linearization at zero equilibrium point to Equation (30), the linearization system of controlled system (30) can be achieved where a ij (i, j = 1, 2, 3) is same as Equation (19).
It is necessary to get the transversality condition, thus we make the following assumption (H9) the expressions of C i , D i (i = 1, 2) are in Appendix D.
Obviously, if hypothesis (H9) is true, then the transversality condition is true. The proof of Lemma 4 is finished.

Remark 2.
The influence of the linear delay feedback controller on the stability domain of system (9) is direct because the formula for calculating the delay critical value includes L (see Appendixes C and D).

Numerical Simulations
In this section, we use the Adama-Bashforth-Moulton predictive correction method [43] to validate the feasibility of theoretical analysis.

Case 1. The influences of fear factors on the stability region
In this paper, the fear factor is considered as two cases caused by mature predators and immature predators, respectively, i.e., k 1 and k 2 . We are interested in which one makes an important role in the stability of system (38). When all parameters and fractional-orders remain unchanged except k 1 , let k 1 increases continuously, we can get different bifurcation critical points (τ 0 , ω 0 ) presented in Table 1:  Next, in a similar way, remain all parameters and fractional-orders unchanged except k 2 , and let k 2 increases continuously, we also get different bifurcation critical points (τ 0 , ω 0 ) presented in Table 2: It can be viewed in Figure 3 that the occurrence of the Hopf bifurcation is put off slightly as k 1 increases. However, the relationship between the critical value of delay τ 0 and fear factor k 2 shows a U-shaped curve. What calls for special attention is when k 2 increases from 0.1 to 0.6, τ 0 descends by 9%. From the perspective of ecology, if the fear of predators is greater, the instability of the system will be more obvious, and the critical value of delay will decrease. In other words, the occurrence of Hopf bifurcation is advanced and the stability state is broken. However, when the intensity of the fear effect reaches a certain level, the limit effect will be produced, and the critical value of delay will decrease more and more weakly. On the other hand, Figure 3 shows that the fear effect on the stability region of the system mainly comes from mature predators.

Case 2. The influence of prey refuge rate on the stability region
Prey refuge is an effective measure of ecosystem regulation. We are interested in how the prey refuge rate m influences the stability of system (38). Same as Case 1, remain all the parameters and fractional-orders unchanged except m, let m increases continuously, we can get different bifurcation critical points (τ 0 , ω 0 ) presented in Table 3.
It can be noticed easily from Table 3 and Figure 4 that when m increases from 0.05 to 0.12, the critical value of delay τ 0 increases from 35.533 to 262.802, i.e., the stability region of the system becomes 7.4 times the original. That is to say, to system (38), m has extremely influence on stability at the coexistence

Case 3. The influence of fractional-orders on the stability region
Let q 1 = 1, q 2 = 1, q 3 = 1 and k 1 = 0, the system (38) becomes an integer-order system corresponding to system (8). The coexistence equilibrium point is (0.2131, 0.0247, 0.1974), and the bifurcation critical point is (τ 0 = 22.77, ω 0 = 0.0501), which is consistent with the results in Wang [35]. If q 1 = 0.98, q 2 = 0.92, q 3 = 0.95 and k 1 = 0, the bifurcation critical point is (τ 0 = 35.46, ω 0 = 0.0389). These results validate that when other parameters of the model are consistent, the delay critical value of the emergence of Hopf bifurcation in the fractional-order system is obviously larger than that in the integer-order system, and the stability domain of the system expands from [0, 22.77) to [0, 35.46). Otherwise, Tables 4-6 further illustrate that fractional-order can effectively expand the stability domain of the system.
Next, we want to know which fractional-order has the obvious effect on the stability of the system (38). The main idea is to keep two fractional-orders unchanged and vary the third one. Tables 4-6 show the different bifurcation critical points (τ 0 , ω 0 ) along with q i (i = 1, 2, 3) varying, respectively. The varying curves are drawn in Figure 5 to compare the distinguishing influences to the stability region.  Table 6. The relationship between q 3 and τ 0 (q 1 = 0.98, q 2 = 0.92). It can be seen from Tables 4-6 and Figure 5 that q 1 has an important influence over q 2 and q 3 on the stability of system (38). That is to say, in an ecosystem such as model Equation (38), the prey is the main fact that affects the stability of the ecosystem. It can be expressed that the change of prey affects the population density not only of prey but also of predator, which intensifies the turbulence of the ecosystem. Moreover, immature predators have more influence than mature predators when fractional-order is less than 0.8 in this ecosystem. It can be seen that when the three fractional-orders tend to 1, respectively, the critical value of delay changes gradually converge. This further shows that it is more practical to use fractional-order to explain the evolution of the ecosystem.

Example 2
Now, we add a linear delay feedback controller to the system to control the emergence of Hopf bifurcation. All the parameters of the controller system are the same as system (38), and the feedback gain is selected as L = −0.01, then the controlled system can be described as The bifurcation critical point of controlled system (39) is (ω * 0 = 0.028071, τ * 0 = 61.544). This means the emergence of Hopf bifurcation is put off obviously, and the stability region is enlarged successfully. The influence of feedback gain L on the bifurcation critical point is illustrated by numeric simulations in Table 7 and Figure 6. It can be seen from Figure 6 that as the feedback gain decreases, the system converges to a steady state faster. In other words, the smaller the feedback gain, the better the control effect of the controller to the Hopf bifurcation.

Conclusions
In this paper, the dynamic behaviors of a fractional-order delay stage structure preypredator model with two fear effects and prey refuge are explored by the linearized method and Laplace transform of fractional-order delay differential equation. Firstly, the conditions for the existence of the coexistence equilibrium point of the system (9) is deduced through the implicit function derivation rule and the function monotonicity theory. Secondly, the stability of the coexistence equilibrium point of the system (9) is investigated with the delay as parameter, and sufficient conditions for the emergence of Hopf bifurcation of the system (9) are obtained. Thirdly, a linear delay feedback controller is added to the system (9) to control the emergence of Hopf bifurcation, and the result states that the system can be controlled successfully by selecting an appropriate feedback gain. Finally, two examples are introduced to validate the theoretical results with the help of numerical simulation.
Moreover, some numerical simulations are performed to explore the influence facts of stability of system (9). The results show that fear factors k 1 , k 2 , prey refuge rate m and fractional-orders q i (i = 1, 2, 3) have distinguish effects to the bifurcation critical value of delay τ, and then affect the stability region of the system. These results have some implications for the regulation and management of ecosystems described in system (9).
However, because of lacking the complete theory of fractional-order differential equation, all the theoretical analyses in this paper are performed on the linearization system of original system (9) and the rationality of theoretical analysis is verified by numerical simulation. We are trying to study the stability of fractional differential equations theoretically in the next work as it is a challenging problem. Acknowledgments: The authors are greatly thankful to the editor and the referees for their valuable suggestions, which help to improve this paper considerably.

Conflicts of Interest:
The authors declare no conflict of interest.