Qualitative Analysis of a Leslie-Gower Predator-Prey System with Nonlinear Harvesting in Predator

This paper deals with the study of the stability and the bifurcation analysis of a Leslie-Gower predator-prey model with MichaelisMenten type predator harvesting. It is shown that the proposed model exhibits the bistability for certain parametric conditions. Dulac’s criterion has been adopted to obtain the sufficient conditions for the global stability of the model. Moreover, the model exhibits different kinds of bifurcations (e.g., the saddle-node bifurcation, the subcritical and supercritical Hopf bifurcations, Bogdanov-Takens bifurcation, and the homoclinic bifurcation) whenever the values of parameters of the model vary.The analytical findings and numerical simulations reveal far richer and complex dynamics in comparison to the models with no harvesting and with constant-yield predator harvesting.


Introduction
Marine life is a renewable natural resource that not only provides food to a large population of humans but also is involved in the regulation of the Earth's ecosystem.The growing human needs for more food and more energy have led to increased exploitation of these resources which affects the Earth's ecosystem.Thus, it is imperative to design harvesting strategies which aim at maximizing economic gains giving due consideration to the ecological health of the concerned Earth's ecological system.Mathematical modeling in harvesting of species was started by Clark [1,2].There are mainly three types of harvesting according to Gupta et al. [3]: (i) ℎ() = ℎ, constant rate harvesting (see [4][5][6][7]), (ii) ℎ() = , proportionate harvesting (see [8,9]), and (iii) ℎ() = /( 1 + 2 ) (Holling type II), nonlinear harvesting (see [10][11][12][13]).Nonlinear harvesting is more realistic and exhibits saturation effects with respect to both the stock abundance and the effort level.Notice that ℎ() → / 2 as  → ∞ and ℎ() → / 1 as  → ∞.
The objective of the present work is to study dynamical behaviors of a Leslie-Gower predator-prey model in the presence of nonlinear harvesting in predators depending on parameters of the model.There have been many papers on the Leslie-Gower predator-prey system with harvesting.For example, May et al. [14] proposed a Leslie-Gower predatorprey model in which two different kinds of constant-yield harvesting applied on both prey and predator species have been considered and this model was studied by Beddington and Cooke [15].Beddington and May [16] proposed and studied Leslie-Gower predator-prey model when both the prey and predators were harvested with constant effort.Beddington and Cooke [15] studied a Leslie-Gower predatorprey model in which the preys are harvested at a constantyield rate and predators are harvested with constant-effort rate.Zhu and Lan [17] studied a Leslie-Gower predator-prey model with constant-yield prey harvesting.Gong and Huang [18] studied the Bogdanov-Takens bifurcation for the model and showed that for different parameter values the model has a limit cycle or a homoclinic loop.Gupta et al. [3] discussed the bifurcation analysis of a Leslie-Gower prey-predator model in the presence of nonlinear harvesting in prey.Huang et al. [19] studied the effect of constant-yield predator 2 International Journal of Engineering Mathematics harvesting on the dynamics of a Leslie-Gower type model and showed that the model has Bogdanov-Takens (BT) singularity of codimension 3 or a weak focus of multiplicity two for some parameter values.They have shown that as the parameters change, the model exhibits saddle-node bifurcation, repelling and attracting B-T bifurcations, and supercritical, subcritical, and degenerate Hopf bifurcations.
This article is organized as follows.In Section 2, we describe the mathematical model in detail.In Section 3, we obtain the number and location of equilibrium points and the local and global stability of the equilibria are investigated.In Section 4, the existence of saddle-node, Hopf, and Bogdanov-Takens bifurcations is shown.In Section 5, we present several numerical simulations that support our theoretical results.Finally, we present some ramification of our results in Section 6.

Model Equations
2.1.Leslie-Gower Model.In general, the Leslie-Gower predator-prey model [20] is given as follows: ) , (, ) ̸ = (0, 0) ,  (0) > 0,  (0) > 0, where  ≡ () and  ≡ () are the prey density and predator density at time , respectively; , , , , and  are positive parameters and represent the intrinsic growth rate of prey, intrinsic growth rate of predator, carrying capacity of prey in the absence of predator, maximal predator per capita consumption rate, and measure of the food quality that the prey provides for conversion into predator births.Model (1) has been studied by Hsu and Huang [21].They showed that the unique positive interior equilibrium point of system (1) is globally asymptotically stable under all biologically admissible parameters.
The number and location of interior equilibrium points have been depicted in Figure 1.Now, we shall discuss the stability of the equilibrium points.The Jacobian matrix of system (5) at the equilibrium point  0 cannot be calculated as the term / is not defined at (0, 0).We use the blow-up technique to analyze the stability of the equilibrium point  0 as given in [23].Let  =  and  = V; then, system (5) reduces to System (7) has two equilibrium points  0 = (0, 0) and  1 = (0, (1/)( − 1 − ℎ/)) whenever  − 1 − ℎ/ > 0. The Jacobian matrix of system (7) at the equilibrium point  0 is Thus, the equilibrium point  0 of system ( 7) is an unstable point as  − 1 − ℎ/ > 0. The Jacobian matrix of system (7) at the equilibrium point  1 is Thus, the equilibrium point  1 of system ( 7) is always saddle as  − 1 − ℎ/ > 0. Thus, we can conclude the discussion above as follows.
Theorem 2. The trivial equilibrium point  0 of system ( 5) is a saddle point.

Theorem 3. (a)
The equilibrium point  is asymptotically stable whenever  < ℎ and a saddle whenever ℎ < .
(b) The equilibrium point  1 , if existent, is always a saddle point.
(c) The equilibrium point  2 , if existent, is asymptotically stable whenever (d) The equilibrium point  3 , if existent, is a degenerate singular point.
(e) The equilibrium point  4 , if existent, is always asymptotically stable.
Proof.(a) The eigenvalues of the Jacobian matrix   are −1 and (1 − ℎ/), so the equilibrium point  is asymptotically stable whenever  < ℎ and a saddle whenever ℎ < .
(b) The determinant det(  1 ) < 0, so the interior equilibrium point  1 is a saddle point.
From Lemma 1, if ℎ = , system (5) has unique interior equilibrium point  4 , and if ℎ < , system (5) has unique interior equilibrium point  5 .In Theorem 3, it is proved that these interior equilibrium points are always asymptotically stable.Now we show that these equilibrium points are globally asymptotically stable.Theorem 4. The equilibrium points  4 and  5 , if existent, are globally asymptotically stable.
Proof.From Lemma 1, system (5) has one trivial equilibrium point  0 , one axial equilibrium point , and one interior equilibrium point  4 whenever ℎ = .Further, it has one trivial equilibrium point  0 , one axial equilibrium point , and one interior equilibrium point  5 whenever ℎ < .Also, from Theorem 3, the trivial equilibrium point  0 is always a saddle point, axial equilibrium point  is a saddle point whenever ℎ ≤ , and the interior equilibrium points  4 and  5 are asymptotically stable.We define the following function: where After simplification, we obtain (, ) = −(1/ 2 )(1+( 2 + (−ℎ)(+2))/(+) 2 ) < 0 as ℎ =  or ℎ < .Thus, by using Dulac's criterion [24], system (5) will not have any nontrivial periodic orbit in Ω.Note that -axis and -axis are the stable manifolds of the trivial and the axial equilibrium points, respectively.Using this in conjunction with the Poincare-Bendixson theorem [24] gives us that the interior equilibrium points  4 and  5 will be globally stable.
In Theorem 3, it is shown that the equilibrium point  3 is a degenerate singular point.Now we study the property of this point.

Theorem 5. The equilibrium point 𝐸
Proof.First, we shall shift the equilibrium point  3 = ( 3 ,  3 ) to the origin by using the transformations  1 =  −  3 and V 1 =  −  3 ; system (5) reduces into the form where If Both eigenvalues of the Jacobian matrix   3 are zero and system (13) reduces to Now we introduce the affine transformations where Now, we consider the  ∞ change of coordinates in the small vicinity of (0, 0): Then, system (16) reduces to Now, we choose the  ∞ change of coordinates in the small neighbourhood of (0, 0): System (19) reduces to International Journal of Engineering Mathematics Now, we choose the final  ∞ change of coordinates in the small neighbourhood of (0, 0): System ( 21) reduces to If 23) is a cusp of codimension 2; that is, the interior equilibrium point  3 of system ( 5) is a cusp of codimension 2.
In Section 4, we shall study the bifurcations occurring in system (5).

Bifurcation Analysis
4.1.Hopf Bifurcation.In Theorem 3, it is shown that the interior equilibrium point  2 is stable whenever In this parametric condition, the equilibrium point  2 is a weak focus or a center.Hence, system (5) may enter a Hopf bifurcation at the point  2 .In this subsection, we consider the parameter  as the Hopf bifurcation parameter and discuss the conditions under which the stability of  2 will change and system (5) exhibits Hopf bifurcations.Theorem 6. System (5) = 0 such that, at  =  * , tr(  ) = 0 and det(  ) > 0. In order to ensure the existence of a Hopf bifurcation, we have to check the transversality condition (/)(tr(  ))| = * ̸ = 0 (see [24]).We have Hence, the transversality condition for a Hopf bifurcation is satisfied.To determine the direction of Hopf bifurcation and stability of  2 , we compute the first Liapunov coefficient of system (5) at the equilibrium point  2 .
Let  =  −  2 and  = V −  2 ; then, the equilibrium point  2 is shifted to the origin (0, 0) and system (5) can be rewritten as where +=4     V  .Hence, using the formula of the first Lyapunov number  at the origin of ( 24), we have where Δ = ( 2 /( +  2 )) √ (1 +  +  − ℎ) 2 − 4ℎ.If  ̸ = 0, then the origin of ( 24) is a weak focus of multiplicity one: also origin is stable when  < 0 and unstable when  > 0. Hence, in parameter space (, , ℎ, ), there exist surfaces , called subcritical and supercritical Hopf bifurcation surface, respectively, such that if the parameter set (, , ℎ, ) is in  sub , the equilibrium point  2 of system ( 5) is a weak focus of multiplicity 1 and is unstable, and if the parameter set (, , ℎ, ) is in  sup , the equilibrium point  2 of system ( 5) is a weak focus of multiplicity 1 and is stable.
From the discussion above, the equilibrium point  2 of system ( 5) is a weak focus of multiplicity 1 and is unstable if (, , ℎ, ) is in  sub .Also from Theorem 3, the equilibrium point  2 is stable whenever ( 2 / 2 )(( 2 −  2 )/( +  2 ) − ) −  2 < 0, that is,  <  * , and is unstable whenever the equilibrium point  2 generates an unstable limit cycle as the bifurcation parameter  passes through the bifurcation value  =  * .From one side of the surface  sub to the other side, system (5) can undergo a subcritical Hopf bifurcation.An unstable limit cycle arises in the small neighbourhood of the equilibrium point  2 whenever (, , ℎ, ) is in  sub , 0 <  <  * and | −  * | ≪ 1.Similarly, a stable limit cycle arises in the small neighbourhood of the equilibrium point  2 whenever (, , ℎ, ) is in  sup ,  >  * and | −  * | ≪ 1.

Bogdanov-Takens Bifurcation.
In Theorem 5, we have proved that the interior equilibrium point  3 is a cusp of codimension 2 whenever ( 3 / 3 )(( 3 −2 3 −)/(+ 3 ))−  3 = 0 and  4  5 ̸ = 0, which implies that there may exist the Bogdanov-Takens bifurcation in system (5).The parameters ℎ and  are taken as the bifurcation parameters as they are important from biological point of view and by means of a series of transformations as given in Xiao and Ruan [5], we derive a normal form.Proof.We consider that the parameters ℎ and  vary in a small neighbourhood of the BT point (ℎ 0 ,  0 ).Let (ℎ 0 + 1 ,  0 +  2 ) be a point of this neighbourhood, where  1 ,  2 are small.System (5) becomes Translating the equilibrium point  3 into the origin by using the transformations  1 =  −  3 and  2 =  −  3 and then using Taylor's series expansion, system (26) reduces to where
The system above is topologically equivalent to the normal form of the Bogdanov-Takens bifurcation which is given by The system undergoes a Bogdanov-Takens bifurcation if =0) ̸ = 0.When system (5) undergoes Bogdanov-Takens bifurcation at  1 =  2 = 0, three bifurcation curves in  1  2 plane through the BT point are given by (1) saddle-node bifurcation curve: (2) Hopf bifurcation curve: (3) Homoclinic bifurcation curve: International Journal of Engineering Mathematics 9
Translating the equilibrium point  3 into the origin by using the transformations  1 =  − 0.343303 and  2 =  − 0.656697, by Taylor's series expansion, system (40) reduces to where where in the small neighbourhood of (0, 0), system (42) reduces to where where (0, ) =  3 and  is a power series in  1 whose coefficients depend on parameters ( 1 ,  2 ).
The local representations of the bifurcation curves are as follows: These bifurcation curves are sketched in a small neighbourhood of the origin in the  1  2 plane (see Figure 5(a)).
The bifurcation curves divide the parameter plane into four parts: I, II, III, and IV.The saddle-node bifurcation curve (SN) is the horizontal axis, that is,  2 = 0 axis, and the homoclinic bifurcation curve (HL) and Hopf bifurcation curve (H) lie in the fourth quadrant.When ( 1 ,  2 ) = (0, 0), then system (5) has a unique interior equilibrium point which is a cusp of codimension 2 (see Figure 5(b)).When the parameters  1 and  2 vary and lie in region I, then system (5) has no interior equilibrium point and all the solution trajectories go to the axial equilibrium point; that is, the axial equilibrium point is globally stable (see Figure 5(c)).Hence, in this case, the predator will go extinct and prey will approach the carrying capacity.When the parameters  1 and  2 pass the SN bifurcation curve and lie in region II, the cusp of codimension 2 breaks into a hyperbolic saddle and an unstable focus (see Figure 5(d)).When the parameters  1 and  2 lie in region III, the cusp of codimension 2 breaks into a stable focus and a hyperbolic saddle.The change of stability of the focus yields an unstable limit cycle (see Figure 5(e)).Further, when the parameters  1 and  2 lie in region IV, the cusp of codimension 2 breaks into a saddle point and a stable focus (see Figure 5(f)).Hence, an open set of initial population densities exists for which both predator and prey approach a stable steady state.

Discussion and Conclusion
In this paper, a Leslie-Gower predator-prey model has been analyzed in the presence of nonlinear predator harvesting.It is shown that the system has at most four equilibrium points in Ω, in which the trivial equilibrium point and the axial equilibrium point always exist while the positive interior equilibrium point changes from two to zero.The qualitative properties of solutions in the vicinity of (0, 0) have been studied by using a blow-up technique, and it is found that the origin will never be stable; ecologically speaking, the two species cannot go to extinction together.The axial equilibrium point is either asymptotically stable or globally stable or a saddle depends on the parametric conditions.Thus, ecologically speaking, the prey species never goes to extinction whatever the choice of initial population density.
If two interior equilibrium points exist, then one is always a saddle point while the other is either asymptotically stable or unstable or the system undergoes a Hopf bifurcation around this point; that is, for certain parametric domain, the system exhibits a bistability phenomenon as well as oscillatory coexistence of both the populations.The stability of limit cycles has been studied and validated through numerical simulations by calculating the first Lyapunov number.It is also found that for certain parametric domain the proposed  International Journal of Engineering Mathematics system has a unique interior equilibrium which is globally asymptotically stable.It is shown that the system can have zero, one, or two interior equilibrium points through saddle-node bifurcation as the bifurcation parameter ℎ crosses a certain threshold value.By means of Sotomayor's theorem, the existence of saddle-node bifurcation has been shown.Ecologically speaking, a maximum threshold value of ℎ exists, below which both species coexist and above which predator species suddenly collapse to extinction.The proposed system undergoes Bogdanov-Takens bifurcation near degenerate equilibria.We have considered the parameters ℎ and  as bifurcation parameters and reduced the system into normal form.Ecologically speaking, a small perturbation in the bifurcation parameters is a cause of extinction, coexistence, and oscillation.
The dynamical analysis of model (5) shows complex and rich dynamics as compared to the model with no harvesting [21] and with constant-yield predator harvesting [19].The model without harvesting has a unique interior equilibrium point which is globally asymptotically stable [21] while model ( 5) has zero to two interior equilibrium points and model ( 5) undergoes a series of bifurcations (Hopf bifurcation, saddlenode bifurcation, and Bogdanov-Takens bifurcation).The model with constant-yield predator harvesting has no axial equilibrium point, and the solutions once touching the -axis will leave the first quadrant [21].But the model with nonlinear predator harvesting has one trivial equilibrium point  0 = (0, 0) and one axial equilibrium point  = (1, 0) and the solutions touching the -axis will go to the axial equilibrium point .Further, if no interior equilibrium point exists, then the axial equilibrium point  is globally asymptotically stable (see Figure 2(e)).Ecologically speaking, in the absence of predator species, the prey density approaches the carrying capacity.Although the model with constant-yield predator harvesting undergoes a series of bifurcations, like Hopf bifurcation, saddle-node bifurcation, and Bogdanov-Takens bifurcation [19], the model with nonlinear harvesting gives more general parametric conditions for the occurrence of these bifurcations.Thus, the results of nonlinear harvesting model can explain the real life situation in a more effective and realistic manner.

Figure 1 :
Figure 1: This diagram shows the number and location of the positive interior equilibrium points of system (5).The green, red, and black color curves are the predator nullclines for different values of ℎ and the straight line is the prey nullcline: (a) ℎ > , for black color parabola ℎ > ℎ, for red color parabola ℎ = ℎ, and for green color parabola ℎ < ℎ; (b) ℎ = ; (c) ℎ < .