Michaelis-Menten-Type Prey Harvesting in Discrete Modified Leslie-Gower Predator-Prey Model

,


Introduction
The basic interaction between various species on this planet is predation and is widely present in nature.The existence of most species in our ecosystem is based on predator-prey relationships, which makes the rich biodiversity of complex ecosystems possible [1].These complexities can aptly be explained using mathematical models and the related dynamics.As Clark has observed in [2], the dynamics of the predator-prey models also provide methods to optimally manage renewable resources, apart from coexistence conditions of predators and preys.For this ecosystem, the most basic model was proposed by Lotka [3] and Volterra [4], in early 20 th century.This model captured the oscillating behavior in populations of a predator and its prey.The simplicity of this model, however, is unable to address most real-world scenarios.To rectify this, many modifications of Lotka-Volterra type model have been proposed by researchers.In [5], Holling proposed various types of functional responses to better understand the components of predator-prey interactions.
An alternative was proposed by Leslie, in [6], in which the number of prey and the carrying capacity of the predator's environment was proportional.Leslie focused on the fact that the increasing capacity of both predator and prey cannot be unbounded.This fact was not incorporated in the Lotka-Volterra model.This model was shown, in [7], to possess globally asymptotically stable and unique positive equilibrium, for any permissible parameters.Then, May in [8] showed its stability with Holling functional responses.L. Chen and F. Chen proved the global stability of the unique interior equilibrium for a Leslie-Gower predator-prey model with feedback controls [9].Uniqueness of limit cycles and the Hopf bifurcation for this model was discussed in [10,11].
The primary drawback of this model is that if the prey population is at low densities, it is not possible for the predator to switch to alternative food source [12].Aziz-Alaoui and Daher [13] rectified this scenario by adding a provisional alternative food source parameter.This model is generally called the modified Leslie-Gower predator-prey model.Caughley [14] modelled the biological control of the prickly-pear cactus by the moth Cactoblastis cactorum, using this system.In [15,16], the authors used this system to model the predator-prey mite outbreak interactions on fruit trees.Much work has been done on this model.For further details, look at [17][18][19][20][21][22][23] and the references therein.In [24], W. C. Allee introduced what was altered termed as the Allee effects, to the Leslie-Gower type models.The stability dynamics and bifurcation analysis of the predator-prey systems subject to Allee effect is discussed in [25][26][27][28].
Harvesting of different biological species is necessary for a variety of reasons.The ever-increasing use of natural resources, either due to the increase in human population or in the name of economic progress, has caused the ecological system of this planet to trip from the equilibrium it was in, for thousands of years.A common outcome of these practices is the extinction of species from the face of this planet.Harvesting of such species is one way to minimize the damage.Harvesting for species is also practices for purely business reasons.But there is a silver lining in this strategy.The more harvesting is done, the more predator species are likely to become extinct [29].Thus, the need to reinforce scientific management of harvesting is much needed.Clark [2] discussed the problem of combined or nonselective harvesting of two ecologically independent populations obeying the logistic law of growth.Multispecies harvesting models are studied in detail by Chaudhuri [30,31], Mesterton-Gibbons [32], Kar and Chaudhuri [33,34], etc. Nonselective harvesting model of a preypredator fishery is studied in detail by Chaudhuri and Ray [35], Kar et al. [36], etc.
Michaelis-Menten-type predator-prey model, alternatively known as the ratio-dependent predator-prey model, uses the idea that the ratio of prey to predator abundance directly affects the per capita predator growth rate.The effectiveness of this model is backed by a plethora of experimental and observational data, for the predator which must compete for prey.For more details, see ( [37][38][39][40][41][42][43]).Michaelis-Menten-type predator-prey model has been rigorously studied in ( [44][45][46][47][48][49]) and references therein.This investigation discovered existence of rich dynamics, like stable limit cycle, multiple attractors, and deterministic extinction, to name a few.Existence of hyperbolic, parabolic, and elliptic orbits was revealed near the origin, as well as combinations of such orbits, for various parametric values.
The continuous predator-prey models have been successfully used to obtain desired results, and the disadvantages of using continuous systems are also quite apparent.By definition, the continuous systems require the subject species have continuous and overlapping generations.This is not generally true, for example, salmon, which have an annual spawning season and are born at the same time each year.For these types of populations which have nonoverlapping generation characteristics, the discrete time models are more descriptive and suitable than the continuous models [50].Also, the discrete models should generate richer and reality-based dynamics compared to the continuous time models [51].In addition, since many continuous models cannot be solved analytically, using difference equations for approximation and finding solution is much practical way to approach the problem.In population biology and complex ecosystems, discrete-time models are used to examine the taxonomic group of organisms and species with the passage of time.These models are best to describe the chaotic behavior of nonlinear dynamics [52,53].
Gupta and Chandra, in [54], proposed and studied modified Leslie-Gower predator-prey model with Michaelis-Menten-type prey harvesting in prey.They proved the permanence and stability and discuss different bifurcations of this model.Their model guaranteed a feasible upper bound of the rate of harvesting for the coexistence of the species.Their model is given in [54] as with the initial conditions xð0Þ > 0, yð0Þ > 0. Here, xðtÞ and yðtÞ represent prey and predator densities, respectively, at time t.The details of the parameters are given in their paper.For our purposes, all parameters are positive.
Another possible way to understand the complex problem of interaction between prey and predator is by using discrete models.We use the forward Euler method to discretize the above system.The discrete counterpart of (1) can be given as In the remainder of this paper, we present the discrete counterpart of the Michaelis-Menten-type prey harvesting in modified Leslie-Gower predator-prey model, ((4)).Our aim is to illustrate the detailed mathematical analysis of the topology of nonnegative interior fixed points, including their existence in Section 2 and stability dynamics in Section 3. The conditions for the existence of flip and Neimark-Sacker bifurcations will be derived in Section 4, using the center manifold theorem and bifurcation theory.The theory will be verified through numerical examples in Section 5. Finally, we will conclude in Section 6.

The Fixed Points
In order to find the fixed points, we need to solve the following system.
The system has four fixed points on the boundary, namely, For positive equilibrium points, we have the following cases.
(i) The system has no positive stationary point, if (ii) If ðα − βðc + 1ÞÞ 2 = 4β 2 h, then the system has a unique positive stationary point, E * = ðx * , m + x * /β Þ, where (iii) If ðα − βðc + 1ÞÞ 2 > 4β 2 h, the system will have two positive stationary points E 5 = ðx m , m + x m /βÞ and E 6 = ðx p , m + x p /βÞ, where As can be seen in Figure 1, the positive interior points may vanish or merge into one another, depending on the values of parameters.

The Topological Dynamics of the Fixed Points
The respective Jacobian matrices, evaluated at each fixed point, are where In order to find the topology and properties of the aforementioned fixed points, the following lemma will be used, which is omnipresent in the text books related to discrete dynamical systems.
Using the above lemma, we can easily see that then it may undergo transcritical or fold bifurcation, and if h = cð1 + 2/kÞ, it undergoes period-doubling bifurcation.Similarly, the fixed representing the prey free scenario,  2 , the fixed points undergo period-doubling bifurcation.For all these boundary fixed points, Neimark-Sacker bifurcation is not possible.
We are much more interested in the behavior of interior fixed points, which are much more richer in comparison and for practical purposes, usable.Theorem 3.For the fixed points E * , E 5 , and Then, (1) The fixed points are locally asymptotically stable if (4) The fixed points are nonhyperbolic with eigenvalues z 1 = 1 and jz 2 j ≠ 1, if and only if h = ðc + xÞ 2 (5) The fixed points are nonhyperbolic with eigenvalues The fixed points are nonhyperbolic, having complex conjugate eigenvalues, with jz Proof.For any positive interior fixed point, E * , E 5 , or E 6 , the characteristic polynomial of the Jacobian matrix at the given fixed point is Then, we can use Lemma 1, since for any k > 0, Then, using Lemmas 1 and 2, we have the following results.
(1) Each fixed point is a sink if and only if (2) Each fixed point is a source if and only if The Jacobian matrix at the given fixed points has eigenvalues z 1 = 1 and jz 2 j ≠ 1 if and only if h = ðc + xÞ 2 (5) The Jacobian matrix at the given fixed points has eigenvalues The eigenvalues are complex conjugates with jz As used in [55,56], we know that due to the jury condition, a transcritical or fold bifurcation is obtained when p j ð1Þ = 0, j ∈ f * ,5, 6g.As shown above, this is only true if h = ðc + xÞ 2 .For the fixed point E * , this condition is equivalent to h = 1/4ðα/β − c − 1Þ 2 , which we gets automatically, since in this scenario, ðα/β − c − 1Þ 2 − 4h = 0. Thus, we always have p * ð1Þ = 0 for the fixed point E * , i.e., a transcritical or fold bifurcation is always obtained for any parametric values, at the fixed point E * .

Bifurcations
For positive fixed points, E 5 and E 6 , define The procedure is identical for finding the normal form of the bifurcations for E 5 and E 6 , as illustrated below.
4.1.Period-Doubling Bifurcation.Let ðα, β, h, m, c, ρ, k 0 Þ ∈ Ω PD and let K denote the perturbation parameter for the mapping (3), i.e., jKj⋘1.Then, variation of parameters in small neighborhood of Ω PD gives emergence of period 5 Journal of Function Spaces doubling bifurcation.Let k 0 be the bifurcation parameter.The perturbed mapping is Define X = x − x and Y = y − y, to translate the fixed point to ð0, 0Þ.Here, ð x, yÞ is the positive interior fixed point.Then, the above mapping is converted to After Taylor series expansion, the above mapping can be written as where a 11 = 1 + kA 1 , a 12 = −kA 2 , a 21 = kA 3 , and a 22 = 1 − kA 4 , and and To convert the system (43) in the normal form, let where T PD is an invertible matrix.Thus, we can convert the above system, such that we have ð0, 0Þ as fixed points, as shown below. and In order to implement the center manifold theorem, we assume that M C be the center manifold of (30), evaluated at ð0, 0Þ in a small neighborhood of K = 0.It can be approximated as follows.
7 Journal of Function Spaces where We need to find G PD ðK, uÞ such that F PD ðK, G PD ðK, uÞÞ = 0. Since we need the conditions in M C to be satisfied, we may assume Thus, we can write By comparing coefficients in the equation F PD ðK, G PD ðK, uÞÞ = 0, we get and The dynamics restricted to M C are given locally by the map Furthermore, define Theorem 5.If Y 1 ≠ 0 and Y 2 ≠ 0, then system undergoes period-doubling bifurcation at the unique positive equilibrium point when parameter k varies in small neighborhood of k 1 or k 2 .Moreover, the period-two orbits that bifurcate from positive equilibrium are stable if Y 2 > 0 and unstable if Y 2 < 0.
Remark 6.Since we have already found the center manifold for our system in (39), we can easily show the existence of transcritical or fold bifurcation, when one of the eigenvalue equals 1 and the other not on the unit circle.
Then, variation of ðα, β, h, m, c, ρ, kÞ in small neighborhood of Ω NS gives emergence to Neimark-Sacker bifurcation.Let K denote the perturbation parameter for the mapping (3), i.e., jKj1.The perturbed mapping is Define X = x − x and Y = y − y, to translate the fixed point to ð0, 0Þ.Here, ð x, yÞ is the positive interior fixed point.Then, the above mapping is converted to After Taylor series expansion, the above mapping can be written as and b i , i ∈ f1, 2, ⋯, 14g are defined as above.Let the characteristic polynomial of matrix in (43) be where Since ðα, β, h, m, c, ρ, kÞ ∈ Ω NS , therefore the roots of (45) are complex conjugate z 1 , z 2 , and jz 1 j = jz 2 j = 1, and it follows that z 1,2 = PðHÞ/2 ± 1/2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4QðHÞ − PðHÞ 2 q and 1 = jz 1 j = jz 2 j = ffiffiffiffiffiffiffiffiffiffiffi QðHÞ p .Note that PðKÞ 2 − 4QðKÞ < 0, which implies that PðKÞ 2 < 4QðKÞ or PðKÞ ∈ ð−2, 2Þ.Also, In order to ensure that the roots of the characteristic polynomial do not lie in the intersection of unit circle of coordinate axis when K = 0, we need to check that z m 1,2 ≠ 1, for m = 1, 2, 3, 4 at K = 0.This is equivalent to checking Pð0Þ ≠ −2, 0, 1, 2. Since ðα, β, h, m, c, ρ, kÞ ∈ Ω NS , we already know that Pð0Þ ≠ −2, 0, 2. Finally, Thus, the roots of ( 45) do not lie in the intersection of unit circle of coordinate axis when K = 0. Now, to convert the system (43) to the normal form, when K = 0, let R = Pð0Þ/2 and S = 1/2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 4Qð0Þ − Pð0Þ 2 q and let Here, T NS is an invertible matrix.Consider the transformation  The system (43) can be written as where and S , e 7 = 0, For the map (52) to undergo Neimark-Sacker bifurcation, the following quantity must not zero [22]:   Using the result in [57] and the above analysis, we have the following result.55) is nonzero, then the model ( 5) undergoes Neimark-Sacker bifurcation at the equilibrium point ð x, yÞ, provided that parameter k changes in the small neighborhood of k 0 = A 1 − A 4 /A 1 A 4 − A 2 A 3 and ðα, β, h, m, c, ρ, k 0 Þ ∈ Ω NS .Moreover, if M NS < 0 (resp., M NS > 0), then an attracting (resp., repelling) invariant closed curve bifurcate from the fixed point ð x, yÞ for k > k 0 (resp., k < k 0 ).

Numerical Simulations
In this section, we will present some examples which will show the presence of period-doubling and Neimark-Sacker bifurcations for the system (3), using specific numerical values of its parameters ðα, β, h, m, c, ρÞ and taking the step-size, k, as the bifurcation parameter.The illustration will be done using bifurcation diagrams, and it will be ratified by showing that Theorems 5 and 7 are satisfied.
After calculations, we get, Y 1 = 0:330205 ≠ 0 and Y 2 = 336:675 ≠ 0. These values further complement the dynamics of our map, observed above and proves the correctness of Theorem 5.Moreover, Y 2 > 0, which shows that the period-two orbits that bifurcate from positive equilibrium are stable.
After calculations, we get, Y 1 = −0:23876 ≠ 0 and Y 2 = 8:71736 ≠ 0. These values further complement the dynamics of our map, observed above and proved the correctness of Theorem 5.Moreover, Y 2 > 0, which shows that the periodtwo orbits that bifurcate from positive equilibrium are stable.From the phase portrait in Figure 6(a), the attractor becomes a stable spiral point and as the parameter k increases, and the stable spiral enlarges in size.As k continues to increase, the spiral becomes larger before changing in to an attracting invariant closed curve with rough edges due to Neimark-Sacker bifurcation.This spiral and edges are plotted in Figure 6   From the phase portrait in Figure 7(a), we can see that the attractor becomes a stable spiral point and as the parameter k increases, the stable spiral enlarges in size (Figure 7(b)).As k continues to increase, the spiral becomes larger before changing into an attracting invariant closed curve with rough edges due to Neimark-Sacker bifurcation.This spiral and edges are plotted in Figure 7(c) at k = 3:74.At k = 3:79, we see continued invariant closed curves around the fixed point, but their edges begin to vanish as plotted in Figure 7(d).

Conclusion
In this paper, we studied the existence and stability of the positive interior fixed points of system (5) and flip bifurcation and Neimark-Sacker bifurcation, under certain conditions, by using central manifold theorem and bifurcation theory.Our main results are given in Theorem 3 and Theorems 5 and 7 and numerical simulations in Section 5.In the details of the result, when the integral step size k is chosen as a bifurcation parameter, the discrete-time modified Leslie-Gower system with Michaelis-Menten-type prey harvesting displays much richer dynamical behavior.It can be seen from Figures 3-9 that the system displays period-1, 2, 4, and 8 orbits.Moreover, Figures 3 and 5 show that the system exhibits period-doubling bifurcation, and the period-two orbits that bifurcate from positive equilibrium are stable, which shows the correctness of Theorem 5. Similarly, Figures 8 and 9 show that the system exhibits Neimark-Sacker bifurcation, with an attracting invariant closed curve bifurcates from the fixed point, proving the correctness of Theorem 7.
Larger attracting invariant closed curves with edges starting to vanish at k = 3:79
(b) at k = 1:795.At k = 1:81, the attracting invariant closed curves become larger with some rough edges as shown in Figure 6(c).It is clear from the close values of the parameter k how quasi-periodic motions are produced around fixed point due to this type of bifurcation.Increasing k furthermore gives rise to continued invariant closed curves