BIFURCATIONS ANALYSIS OF LESLIE-GOWER PREDATOR-PREY MODELS WITH NONLINEAR PREDATOR-HARVESTING

. In the present paper the dynamics of a Leslie-Gower predator-prey model with Michaelis-Menten type predator harvesting is studied. We give out all the possible ranges of parameters for which the model has up to ﬁve equilibria. We prove that these equilibria can be topological saddles, nodes, foci, centers, saddle-nodes, cusps of codimension 2 or 3. Numerous kinds of bifurcations also occur, such as the transcritical bifurcation, pitchfork bifurcation, Bogdanov-Takens bifurcation and homoclinic bifurcation. Several numerical simulations are carried out to illustrate the validity of our results.


1.
Introduction. In our real life, along with the increase of people's demand for food and resources, the exploitation of some biological resources is increasing. In fact, the exploitation of biological resources and the harvesting of populations are commonly practiced in fishery, foresty, and wildlife administration [6,8]. Hence it is important to find a sustainable development policy to protect the ecosystem. The predator-prey models are widely used to describe the interaction of two species, one specie feeds on another. The dynamical properties of the predator-prey models not only can be used to analyze the relations between the prey and predator or to predicate whether they can coexist, but also get insight into the optimal management of renewable resources [8,30]. Therefore, in recent years, predator-prey model has becomes one of the most popular areas in biological systems and many classic models have been proposed [34,24,32,2,10,21].
In this paper, we are interested in the following new predator-prey system where x and y represent population densities of prey and predator, respectively. The parameters r 1 , r 2 , K, a, b, q, E, m 1 , m 2 are positives. And, the prey grows with intrinsic growth rate r 1 and carrying capacity K in the absence of predation. r 2 and a stand for the intrinsic growth rate of the predator, and the maximal predator per capita consumption rate, respectively. b is a measure of food quality that the prey provides towards the predator births, bx takes on the role of a prey-dependent carrying capacity for the predator. q is the catchability coefficient, E is the effort applied to harvest individuals and m 1 , m 2 are suitable constants. Moreover, the ratio r 2 /q is known as the biotechnical productivity of the predator species.
As we know, one of important predator-prey models is the Leslie-Gower model [34], which is a modification of the Lotka-Volteera model [31]. The Leslie-Gower type model can be described by the autonomous bi-dimensional, where the prey grows logistically and the interaction between prey and predator is expressed by Holling type-I functional response [16]. While the subsistence of predator depends exclusively on prey population, hence the conventional environmental carrying capacity of predator is taken to be proportional to prey abundance. For the considerations, many researchers [17,10,28,29] investigated the system The stability of the interior equilibrium is studied in [34] by numerical methods. Lindstrom [28] investigated the nonexistence, existence and limit cycles. Hsu and Huang [17] prove that all the solutions are bounded and positive if their initial values are in the first quadrant, and study the globally asymptotical stability of the interior equilibrium using Liapunov function and LaSalle's invariance principle.
For the commercial purpose and the economic interest, people introduced various harvesting in the predator-prey model [6,27,36]. Since the predator-prey model with harvesting has richer dynamics, then the effect of harvesting on the dynamics of predator-prey system and the role of harvesting in the management have attracted great attentions [3,35,25,11]. Generally, there are three basic types of harvesting reported in the literatures: (i) constant-yield harvesting can be described by a constant, which independent of the size of the population under harvest [18,37]; (ii) constant-effort harvesting that means the number of individuals harvested per unit of time is proportional to the current population [27,26]; and (iii) nonlinear harvesting or called Michaelis-Menten type harvesting [17,14], which is a Holling type-II function about the current population. Amongst the three types of harvesting, Michaelis-Menten type harvesting is more realistic from biological and economic point of view. It shows that the nonlinear harvesting function exhibits saturation effect with respect to both the stock abundance and the effort level.
Gupta et al. [13] studied following Leslie-Gower predator-prey model with Michaelis-Menten type prey harvesting The system has up to five equilibrium including the origin. The origin can be attractor and saddle. Other equilibria can be saddles, nodes, focus, centers and saddlenode. Some bifurcations of codimension 1 were discussed, such as saddle-node bifurcations and Hopf bifurcations. One year later, Gupta and Chandra [15] analysed the dynamics of a modified Leslie-Gower predator-prey model with Michaelis-Menten type prey harvesting.
In this paper, we consider Leslie-Gower predator-prey model with Michaelis-Menten type predator harvesting (1). We find that there are many interesting results such as the origin is always a saddle, some equilibria are the cusps of codimension 2 and of codimension 3, and system (1) undergoes very rich bifurcations. For mathematical simplicity, we do some rescalings in system (1). Let s = r 2 t, u = x/K, v = ay/r 1 . After rewriting (u, v, s) back to (x, y, t), system (1) becomes We investigate the positive equilibria and their phase portraits of system (4). Let h 1 = 1 + c + 2α + αc − 2 α(1 + α)(1 + c). When h > h 1 , there are no interior equilibria, which implies that the extinction of prey or predator occurs. When h ≤ h 1 , there are one or two interior equilibria, depending on the ranges of the positive parameters. It is shown that the maximum harvesting rate h MSY of system (4) is depend on the values of 1/α i.e. the food quality that the prey provides towards the predator births. When α < 1/c, i.e. the food is enough for the predator births, then h MSY = h 1 . When α > 1/c, i.e. the food is not enough for the predator births, then h MSY = c < h 1 . When h ≤ h 1 , we obtain the phase portraits of (4) near each equilibria. We shall show that they are topological saddles, nodes, foci, centers, saddle-nodes of codimension 1, nonhyperbolic stable node of codimension 2, cusp of codimension 2 and cusp of codimension 3.
When the equilibria are nonhyperbolic, we give bifurcation analysis for system (4). System (4) has rich bifurcations. It is shown that system (4) can exhibit numerous kinds of bifurcation phenomena in terms of the original parameters in the model, such as the transcritical bifurcation, pitchfork bifurcation, saddle-node bifurcation, Hopf bifurcations, Bogdanov-Takens bifurcation and homoclinic bifurcation. Those dynamical behaviors can be better for us to understand how the nonlinear predator harvesting affects the dynamics of system (4).
The organization of this paper is as following. In Section 2, we study the various conditions for which the existence of the equilibrium of the model. The phase portraits of each equilibria such as saddle, node, focus, saddle-node, cusp of codimension 2 and cusp of codimension 3 are given in Section 3. In section 4, we consider the Hopf bifurcations of the model depending on all parameters. By computing the first Lyapunov number, the subcritical and supercritical Hopf bifurcations are obtained. Some numerical analysis are given to illustrate these Hopf bifurcations. In section 5, the Bogdanov-Takens bifurcation of codimension 2 is studied. This paper ends by a brief discussion in section 6.
2. Equilibria of the model. In the following, we will first prove some basic property of system (4), such as the positivity and boundedness of solutions as well as the permanence of system (4). Recall that system (4) is said to be permanent if there exist positive constants ω 1 and ω 2 , such that each positive solution (x(t), y(t)) of system (4) with initial condition (x 0 , y 0 ) ∈ Int Ω satisfies, min{lim inf t→+∞ x(t), lim inf t→+∞ y(t)} ≥ ω 1 and max{lim sup t→+∞ x(t), lim sup t→+∞ y(t)} ≤ ω 2 .
We will use the following Lemma which is proved in [4] to get the boundedness of solutions of (4).
Proposition 1. (i) The first quadrant is positive invariant for system (4).
(ii) All solutions (x(t), y(t)) of system (4) with the initial condition (x 0 , y 0 ) in the first quadrant are bounded for all t ≥ 0.
Proof. (i) Letx(t) be a solution ofẋ = γx(1 − x) andỹ(t) ≡ 0. Then (x(t),ỹ(t)) is a solution of (4). Hence the x-axis is a positive invariant of (4). By the interpretation of the Leslie-Gower terms when x = 0 described earlier, we have y ≡ 0. We see that any solution of (4) starting at the first quadrant can not cross y-axis to the second quadrant. Then the first quadrant is positive invariant.
Hence, for the sufficiently large t, the first equation of system (4) can writė where ζ 1 = 1 − (1 + ε)/α > 0. By Lemma 2.1, we have lim inf t→+∞ x(t) ≥ ζ 1 . Thence, x(t) ≥ ζ 1 /2 for sufficiently large t. Now using positivity of y and for sufficiently large t, from the second equation of system (4) we havė Then, combine with the proof of (ii), we can choosing ω 1 = min{ζ 1 , ζ 1 ζ 2 /(2α)} and ω 2 = max{1, M 1 /α}, we get the permanence of the system (4). Remark 1. From Proposition 1.(iii), when α > 1 and c > h, i.e. r 1 > abK and 1 < m 1 (r 2 /q), system (4) is permanent. This implies that when the intrinsic growth rate r 1 of prey and the biotechnical productivity of the predator r 2 /q are large, then predator and prey will persist and extinction will not occur. Biologically, the density of prey is large, then the predator has enough food to get. Hence the prey and predator can coexist.
To find the equilibria in IntΩ suffices to solve the following system: Through direct calculations, we have the following results.
The following results give all the equilibria of system (4) in IntΩ.
(3) If α < 1/c and c < h < h 1 , then (4) has two equilibria E 3 = (x 3 , y 3 ) and  From Theorem 2.3, we seen that the maximum harvesting rate h MSY of system (4) is depend on the values of the food quality that the prey provides towards the predator births 1/α. When α < 1/c, i.e., the food quality is good for the predator births, the h MSY = h 1 . When α > 1/c, i.e., the food quality is bad for the predator births, the h MSY = c < h 1 .
3. Phase portraits of the equilibria. Now, we consider the dynamics of system (4) in the neighborhood of each equilibrium. Let (x, y) be an equilibrium of (4). The Jacobian matrix A(x,y) of system (4) at (x,y) is important in discussing the dynamics of (4). Through direct calculations, we have the Jacobian matrix A(x, y) at an equilibrium (x, y) of system (4).
We first consider the dynamical properties of the equilibria on the boundary. At the equilibria E 0 = (0, 0), the Jacobian matrix cannot be calculated directly because the ratio y x is not defined at E 0 . To understand the dynamical behaviors of (0, 0), as in [19], we need a nonlinear transformations to remove the singularity and expand the system to the whole x − axis. Our goal is to introduce transformation x = u, y = uv and expand the system to the whole x − axis.
Proof. Since (4) is singular at (0, 0), we introduce the transformation to consider the dynamics of (0, 0) of (4) such that its the dynamics can be obtained from the ones of the equilibria of the transformed system on the v-axis.
Remark 2. From Theorem 3.2, when we take h as a bifurcation parameter, then system (4) undergoes a transcritical bifurcation around E 1 if h = c = 1/α and a pitchfork bifurcation around From (9) we know when h = c, the matrix A = Df (X 0 , c) has a simple eigenvalue λ = 0 with eigenvector v = (1, −1) T and A T has an eigenvector w = (0, 1) T corresponding to the eigenvalue λ = 0. By computation, when From the ecological view, the transcritical bifurcation and pitchfork bifurcation occurs around E 1 gives the maximum threshold for continuous harvesting without the extinction risk of the predator species.
In fact, from the formulas of x 3 and r(x 3 ), we have r(x 3 ) = √ ∆ > 0. And, through direct calculations, we have the following results.
Remark 3. The so-called 'biological control paradox' states that we cannot have a low and stable prey equilibrium density. By Theorems 3.4 and 3.6, we see that the E 4 is a saddle and E 3 may be a stable node or focus, this contradicts with the 'biological control paradox' because x 4 > x 3 . From the formula of h 0 , we have that h 0 ≥ c since h 0 − c = (1−αc) 2 (1+2α)(2+c) ≥ 0. In fact, we always have that h 0 ≥ c since h 0 −c = (1−αc) 2 (1+2α)(2+c) ≥ 0. By Theorems 3.2(2) and 3.6(1)(ii), (iii), we know that the system will have two stable equilibrium (E 1 and E 3 ) for some parameters, i.e. the bi-stability case occurs. That is there are some regions such that the predator will go extinction, but for some other regions the prey and predator can co-exist.  Figure  5. The bistability occurred.
In order to study the Hopf bifurcation, we need some knowledge of the Hopf bifurcation theory of the autonomous system. To state the result, we consider the following system ẋ = a 10 x + a 01 y + p(x, y) := f (x, y) y = b 10 x + b 01 y + q(x, y) := g(x, y) where p(x, y) = ∞ i+j=2 a ij x i y j and q(x, y) If D = a 10 b 01 − a 01 b 10 > 0, then it follows from the formula (3 ′ ) of section 4.4 in [33] that the first Liapunov number, denoted by σ is given by where ξ 1 = a 10 b 10 (a 2 11 + a 11 b 02 + a 02 b 11 ), ξ 2 = a 10 a 01 (b 2 11 + a 20 b 11 + a 11 b 02 ),  (1) If σ < 0 (or σ > 0), then the equilibrium (0, 0) is a stable (or unstable) center or stable (or unstable) focus with multiplicity one.
Theorem 4.2. Assume that α < 1/c, h 0 < h < h 1 and γ = ω 0 . Then (1) If σ < 0, then a supercritical Hopf bifurcation of (4) occurs at the equilibrium E 3 and an stable limit cycle will occurred near the E 3 as the bifurcation value µ = a 10 + b 01 increases from zero.
(2)If σ = 0, then system (4) has at least two limit cycles for some suitable parameter values.
(3) If σ > 0, then a subcritical Hopf bifurcation of (4) occurs at the equilibrium E 3 and an unstable limit cycle will occurred near the E 3 as the bifurcation value µ = a 10 + b 01 increases from zero.  1), we know that E 1 is a stable node, E 4 = (0.5212, 0.4788) is a saddle and system has a stable limit cycle near E 3 = (0.3323, 0.6677) this means E 3 is a weak focus of multiplicity 1 and unstable, which is shown in Fig.6.
Theorem 5.1. When α < 1/c, h = h 1 , γ = γ 2 and α > α(c) (α < α(c)), then system (4) have a unique interior equilibrium E 2 which is a cusp of codimension 2 i.e. BT singularity. If we choose h and γ as bifurcation parameters, then system (4) undergoes repelling (attracting) BT bifurcation in a small neighborhood of the interior equilibrium E 2 as (h, γ) varies near (h 1 , γ 2 ). Thence, there exist some parameter values such that system (4) has an unstable (stable) limit cycle, and there exist some other parameter values such that system (4) has an unstable (stable) homoclinic loop.
Proof. We choose h and γ as two bifurcation parameters. Consider the following perturbed system of (4).
Remark 4. By Theorems 5.1, we can seen that these bifurcations can lead to a potentially dramatic shift in the system dynamics, hence it's ecologically important. Through a saddle-node bifurcation, the system can have zero, one or two interior equilibrium as the harvesting rate crosses its critical value. Thus, there are some values of parameters such that the prey and the predator co-exist in the form of a positive equilibrium for different initial values. By a Hopf bifurcation, system can have at least a limit cycle. So, there are some values of parameters such that the prey and the predator co-exist in the form of a positive equilibrium for all initial values lying inside the periodic orbit or an periodic solution with a finite period for all initial values on the periodic orbit. Through a homoclinic bifurcation, we know that the system (4) will have a homoclinic loop. Thence, there are some values of parameters such that the prey and the predator co-exist in the form of a positive equilibrium for all initial values lying inside the homoclinic loop or a periodic orbit with infinite period for all initial values on the homoclinic loop.