Bogdanov–Takens bifurcation of a Holling IV prey–predator model with constant-effort harvesting

A prey–predator model with constant-effort harvesting on the prey and predators is investigated in this paper. First, we discuss the number and type of the equilibria by analyzing the equations of equilibria and the distribution of eigenvalues. Second, with the rescaled harvesting efforts as bifurcation parameters, a subcritical Hopf bifurcation is exhibited near the multiple focus and a Bogdanov–Takens bifurcation is also displayed near the BT singularity by analyzing the versal unfolding of the model. With the variation of bifurcation parameters, the system shows multi-stable structure, and the attractive domains for different attractors are constituted by the stable and unstable manifolds of saddles and the limit cycles bifurcated from Hopf and Bogdanov–Takens bifurcations. Finally, a cusp point and two generalized Hopf points are found on the saddle-node bifurcation curve and the Hopf bifurcation curves, respectively. Several phase diagrams for parameters near one of the generalized Hopf points are exhibited through the generalized Hopf bifurcation.


Introduction
The prey-predator model based on Lotka-Volterra model is one of the most popular models in mathematical ecology and has been widely applied in understanding population dynamics of the species, which is characterized by the complicated interaction among the species and the interaction between the species and their surroundings. Functional response is the rate of prey consumption by the predators, which can characterize the interaction between the prey and the predators. The initial functional responses are presented by Holling in [11,12] and have been classified into three types called Holling I, II and III. Generally speaking, Holling type I response is seen mostly at filter-feeding predators, such as mollusks, algae, and cells. Holling types II and III responses are applicable to the invertebrates and complicated vertebrates, respectively. The common feature of these functions is all monotonic increasing and bounded, which means that the more the prey is in the environment, the better for the predators [6]. However, many experiments and observations indicate that monotonicity does not always hold and a non-monotonic response may exist when the nutrient concentration reaches a high level [3,9,31]. It is also explained that when the number of the prey is large enough, their group defense ability and camouflage ability are strengthened to decrease the predation ability of the predators. For example, large swarms of insects make individual identification difficult for their predators. A lone musk ox can be successfully attacked by wolves, while large groups of musk oxen are less likely to be attacked successfully [28]. Filamentous algae are often qualified as inedible by herbivorous zooplankton [8]. To model such an inhibitory effect, a Holling IV response mx a+bx+x 2 and a simplified Holling IV response mx a+x 2 have been proposed by Andrews [1] and Sokoland [25], respectively. Since then the prey-predator system with Holling IV functional response has been widely studied by many researchers [4,15,20,30,32].
For commercial and economical purposes, the exploitation of biological resources and harvesting of populations are commonly practiced in fishery, forestry and wildlife management [22]. It is known that harvesting, as a scientific management of renewable resources, has a strong effect on the dynamic evolution of the biological models, which can determine the survival and sustainable development of the species. In order to reasonably govern renewable resources, the study of the exploitation of biological resources and harvesting of populations is becoming a hot topic in the field of biological economics, which is related to the optimal management of renewable resources [5].
Recently, it has become very popular among many researchers in using bifurcation theory to analyze the harvested prey-predator model to reach an optimal management in the field of mathematical biology [10, 13, 15-18, 21, 29, 33, 34], which can predict the species' evolutionary direction, including the constant state, periodic state and chaotic state, as the critical biological parameters vary.
From the present work, we find that most of the harvested prey-predator models which have been considered by many authors are confined to three aspects: (i) constant-yield harvesting, (ii) harvesting subjected to only one species or (iii) the numerical response is proportional to the functional response, where the numerical response is the change in predator density as a function of change in prey density [26]. It is the reproduction rate of a consumer as a function of food density. That is, when the predators consume their prey, the numerical response is dynamic relation between the increasing number of the predators and the consumed prey. From the convention, we know that harvesting frequently varies with season, market demand, the species quantity, harvesting cost and so on; the prey and the predators can be both valuable from the business point; the numerical response can be uncorrelated with the functional response. Based on above analysis and the virtue of the Holling IV response, in this paper we extend previous work from three aspects and establish a prey-predator model with Leslie-Gower and Holling IV schemes with constant-effort harvesting, which is given bẏ (1.1) The logistic model r 1 x(1 -x K ) is used to describe the growth of the prey in the absence of predators in the limited environmental carrying capacity. The second term mx b+x 2 y describes the change in the density of the prey attacked per unit time as the prey density changes.
The last term c 1 x is the constant-effort harvesting of the prey. r 2 y(1 -y sx ) is the predator growth term, which is described by the Leslie-Gower function. y sx measures the loss in the predator population due to rarity (per capita of y x ) of its favorite term [19]. The term c 2 y is the constant-effort harvesting of the predators. The meanings of the parameters are given as follows: (i) x and y denote the densities of the prey and the predators.
(ii) K is the environmental carrying capacity for the prey.
(iii) m denotes the maximal predation rate.
(iv) b is the so-called half-saturation constant.
(v) s is a measure of food quality that the prey provides for conversion into the predator birth.
(vi) r 1 and r 2 are intrinsic growth rates of the prey and predators, respectively.
(vii) c 1 and c 2 measure the effort being spent by a harvesting agency (harvesting efforts).
The model (1.1) without harvesting (i.e., c 1 = c 2 = 0) has been discussed in [7,14,20]. Li and Xiao [20] proved that the model could simultaneously undergo a Bogdanov-Takens bifurcation and a subcritical Hopf bifurcation in the small neighborhoods of two different equilibria, respectively. It was shown that for different parameters the model could have a stable limit cycle enclosing two equilibria, or an unstable limit cycle enclosing a hyperbolic equilibria, or two limit cycles enclosing a hyperbolic equilibrium. In Ref. [14], for the same model, Huang et al. investigated the degenerate Bogdanov-Takens bifurcation of codimension 3 around the degenerate positive equilibrium. Some other abundant dynamical phenomena could emerge, such as the coexistence of three hyperbolic positive equilibria, a homoclinic loop enclosing an unstable limit cycle, or a stable limit cycle enclosing three unstable hyperbolic positive equilibria. From the above-mentioned work, Hopf cyclicity and the global dynamics of the same model were investigated by Dai and Zhao [7] for the cases of one non-degenerate positive equilibrium and three distinct positive equilibria, respectively. Some explicit conditions for the globally stable equilibrium were established by applying Dulac's criterion and constructing the Lyapunov function. It was shown that the Hopf bifurcation could occur at each weak focus and more complicated and new dynamics were observed. When the system had a unique positive equilibrium, there existed parameter values such that the system had two limit cycles around it. When the system had three positive equilibria, one limit cycle could bifurcate from each of the two positive anti-saddles simultaneously.
In order to discuss if the added harvesting terms can affect the bifurcation structure of the prey-predator model, we also do a series of bifurcation analysis for the model with liner harvesting. We find that, besides some original dynamical behaviors, some new phenomena and bifurcations appear in the harvested model, such as a globally stable equilibrium, a cusp bifurcation and a generalized Hopf bifurcation. From the cusp point, two saddle-node bifurcation curves are bifurcated, through which three equilibria collide to form an equilibrium or one equilibrium splits into three equilibria. From the generalized Hopf bifurcation point, a fold bifurcation curve of limit cycles is bifurcated, on which two limit cycles with different stability collide to form a semi-stable cycle. These dynamical behaviors have not been found or discussed in the original model in Ref. [7,14,20].
Before entering the topic, we carry out the following scaling transformations: Dropping the bars system (1.1) becomeṡ The rest of this paper is organized as follows. In Sect. 2, we analyze the number and property of the equilibria by using the root formula of the cubic equation. In Sect. 3, we discuss the subcritical Hopf bifurcation near the multiple focus and the Bogdanov-Takens bifurcation near a BT singularity. The characteristics of the cusp point and the generalized Hopf points are analyzed simply. Numerical simulations support our results of theoretical analysis. Finally, we end this paper with a conclusion in Sect. 4.

Types of equilibria
In this section, for system (1.2) we will study the number and type of equilibria in the region D. It is clear that if h 1 < 1 system (1.2) has a boundary equilibrium E 0 = (1h 1 , 0), which is a sink for δ < h 2 and a saddle for δ > h 2 . The corresponding phase portraits are displayed in Fig. 1.
Obviously, a positive equilibrium E * (x * , y * ) of system (1.2) should satisfy the following equations: (2.1) From Eqs. (2.1), we find that, if positive equilibria exist, then h 2 < δ and x * is a root of the equation To determine the type of E * (x * , y * ), the Jacobian matrix J evaluated at E * (x * , y * ) is given by .

Thus
Det According to the sign of Det(J(E * )), the equilibrium E * (x * , y * ) can be divided into three types: an anti-saddle equilibrium, a hyperbolic saddle or a degenerated equilibrium if Det(J(E * )) > 0, Det(J(E * )) < 0 or Det(J(E * )) = 0, respectively. Obviously, the type of E * (x * , y * ) is determined by the sign of F (x * ).
Using the root formula of the third-order equation and estimating the sign of F (x * ), we can get the following results.
(a) If > 0, then system (1.2) has a unique positive equilibrium E * (x * , y * ), which is an anti-saddle equilibrium whatever the sign of A is.
and a degenerated equilibrium The detailed proof of Lemma 2.1 is given in the appendix. The phase portraits for the four cases are exhibited in Fig. 2.
Moreover, (i) E * 1 is an unstable multiple focus with multiplicity one. (ii) E * is a codimension 2 BT singularity. The phase diagram is shown in Fig. 3.

Bifurcations
In this section, we study the subcritical Hopf bifurcation in a neighborhood of E * 1 and the Bogdanov-Takens bifurcation in a neighborhood of E * . Now we carry out bifurcation analysis for system (1.2) by choosing h 1 and h 2 as bifurcation parameters. When the conditions of Theorem 2.2 are satisfied, the unfolding system of (1.2) is given bẏ where μ 1 and μ 2 are small parameters. Through analysis, we have the following conclusions. Proof If μ 1 = μ 2 = 0, then (3.1) has two equilibria E * 1 and E * , whose types are discussed in Theorem 2.2.
When (μ 1 , μ 2 ) lies in region I, system (3.1) has an unstable focus E * 1 . By Lemma 1.1 and the Poincaré-Bendixson theorem, it can be proved that there exists a stable limit cycle. The corresponding phase portrait is depicted in Fig. 5. From the biological point of view, we know that the prey and the predators can coexist and will periodically fluctuate along this cycle eventually.
When (μ 1 , μ 2 ) lies on the curve H 1 , the stable limit cycle still exists and the equilibrium E * 1 turns into an unstable multiple focus with multiplicity one. When (μ 1 , μ 2 ) crosses H 1 Figure 5 For (μ 1 , μ 2 ) = (0.01, 0.007) in region I, there exists an unstable focus E * 1 surrounded by a stable limit cycle Figure 6 When (μ 1 , μ 2 ) = (0.01, 0.001) lies in region VI, a stable focus E * 1 is surrounded by an unstable limit cycle, which is circled by a stable cycle and enters into the region VI, E * 1 becomes a stable focus and an unstable limit cycle bifurcates, which indicates a subcritical Hopf bifurcation happens at E * 1 when (μ 1 , μ 2 ) goes through H 1 . Figure 6 shows that there is a bistable state (a stable equilibrium E * 1 and a big stable limit cycle) and the unstable limit cycle bifurcated from Hopf bifurcation acts as the separatrices of attractive domains for different attractors. Eventually, trajectories can either tend to an equilibrium or a limit cycle depending on initial values of the species.
When (μ 1 , μ 2 ) lies on the curve SN -, a new saddle-node E * emerges and the other equilibrium E * 1 is still a stable focus. When (μ 1 , μ 2 ) passes through SNand enters into region V, there are two new equilibria E * 2 and E * 3 bifurcated from the saddle-node bifurcation. E * 2 is a hyperbolic saddle and E * 3 is a stable focus. The type of E * 1 is not changed and there is still an unstable limit cycle surrounding E * 1 . The Poincaré-Bendixson theorem implies that there is a big stable limit cycle surrounding the three equilibria and the small limit cycle, which is displayed in Fig. 7. It is shown that there is a tristable state in system (3.1). The stable and unstable manifolds of the saddle E * 2 and the unstable limit cycle act as separatrices of different attractive domains. With different initial values, the species will eventually tend to a constant state or oscillate along a big limit cycle.
When (μ 1 , μ 2 ) lies on the curve H 1 , system (3.1) has three equilibria E * 1 , E * 2 and E * 3 , where E * 1 becomes an unstable multiple focus with multiplicity one and the types of E * 2 and E * 3 are unchanged. When (μ 1 , μ 2 ) crosses H 1 into region IV, the unstable limit cycle enclosing E * 1 disappears and E * 1 turns into an unstable focus, which shows that system (3.1) undergoes a subcritical Hopf bifurcation when (μ 1 , μ 2 ) passes through H 1 . The corresponding phase portrait is displayed in Fig. 8, in which E * 3 is a unique globally asymptotical attractor. When (μ 1 , μ 2 ) lies on the curve HL, there appears an unstable homoclinic loop to the saddle E * 2 which is displayed in Fig. 9(a), and a large stable limit cycle encloses all the equilibria and the homoclinic loop shown in Fig. 9(b). Moreover, when (μ 1 , μ 2 ) deviates (b) A big stable limit cycle encloses all the equilibria from the curve HL, the homoclinic loop breaks, which means that system (3.1) undergoes a homoclinic bifurcation when (μ 1 , μ 2 ) passes through the curve HL. When (μ 1 , μ 2 ) lies in region III, there occurs an unstable limit cycle surrounding the stable focus E * 3 , which means that system (3.1) undergoes a subcritical Hopf bifurcation when (μ 1 , μ 2 ) passes through the curve H. For (μ 1 , μ 2 ) = (0.01, 0.009952), there are three equilibria, a small limit cycle and a large stable limit cycle enclosing all the equilibria and the small cycle, which is shown in Figs. 10(a) and 10(b). In this case there is also a bistable state (a stable equilibrium E * 3 and a big stable limit cycle). When (μ 1 , μ 2 ) goes through H and enters into region II (the region between H and SN + ), the unstable limit cycle disappears and E * 3 becomes an unstable focus while the types of E * 1 and E * 2 are not changed. When (μ 1 , μ 2 ) lies on the curve SN + , two equilibria E * 2 and E * 3 overlap and become a saddle-node E * , which means that system (3.1) undergoes a saddle- node bifurcation when (μ 1 , μ 2 ) passes through the curve SN + . Whatever (μ 1 , μ 2 ) lies in region II or on SN + , there both exists a stable limit cycle enclosing all the equilibria by the Poincaré-Bendixson theorem. The phase portrait of system (3.1) for (μ 1 , μ 2 ) on the curve SN + is exhibited in Fig. 11. From the strategy of the optimal management of renewable resources, we find that the saddle-node bifurcation curve acts as the feasible upper bound for the rescaled harvesting efforts which guarantees the coexistence and sustainable development of the species. From Fig. 12(a), we find that when (μ 1 , μ 2 ) varies along the saddle-node bifurcation curve SN + as shown in Fig. 4, a cusp point CP can be encountered and a cusp bifurcation will happen when (μ 1 , μ 2 ) varies near the CP point. When (μ 1 , μ 2 ) varies along two Hopf bifurcation curves H and H 1 as indicated in Fig. 4, respectively, two generalized Hopf points GH + and GHcan be encountered, at which both of the first Lyapunov coefficients are zero. The associated coordinates and parameter values for the two points and the second Lyapunov coefficients l ± 2 of the generalized Hopf l -2 = -2300.709.
From Fig. 12(b), we find that the saddle-node bifurcation curve SN c originating from the CP point and the saddle-node bifurcation curve SN + originating from the BT point are nearly coincident or tangent at the point CP. Through the curve SN c three equilibria collide to form an equilibrium or one equilibrium splits into three equilibria. Similarly, the Hopf bifurcation curve H (H 1 ) originating from the BT point and the Hopf bifurcation curve H + (H -) originating from the GH + (GH -) point are almost coincident and tangent at the GH + (GH -) point. In addition, near the Hopf bifurcation curves H + and H -, two fold bifurcation curves of the cycles T + and Talso appear, on which two limit cycles with different stability collide to form a semi-stable cycle. Several typical phase portraits for (μ 1 , μ 2 ) near the generalized Hopf point GH + are displayed in Fig. 13.

Conclusions
The prey-predator model in [20] added with constant-effort harvesting has been investigated in this paper. The Bogdanov-Takens bifurcation, the cusp bifurcation and the generalized Hopf bifurcation are discussed. By analysis, it is shown that this model can generate many novel dynamic behaviors compared with the model with no harvesting. For example, a globally asymptotically stable equilibrium can appear, while the stable equilibria in the original model are all locally asymptotically stable. Through the generalized Hopf bifurcation, the second Lyaponov coefficient can determine the relative position between the stable limit cycle and the unstable one. A cusp point CP has also been detected, from which two saddle-node bifurcation curves of equilibria emanate, through which three equilibria collide to form an equilibrium or one equilibrium splits into three equilibria. These bifurcations have not been discussed in the original model. In addition, compared with the original model, the harvested prey-predator model can exhibit different bifurcation structure in the parameter plane. These phenomena show that the added harvesting terms play an important role in directing the evolutional directions of the species. From the Hopf and Bogdanov-Takens bifurcations, we find that two small limit cycles bifurcating from two different equilibria are both unstable, while the big limit cycle is always stable, which shows that the ultimate numbers of the species circulate periodically along the big cycle rather than the small cycle.