Bifurcation analysis and chaos control in discrete-time modified Leslie–Gower prey harvesting model

We investigate the dynamical behavior of a modified Leslie–Gower prey–predator model with harvesting in prey population. In order to explore rich dynamics of the model, Euler approximation is implemented to obtain a discrete-time modified Leslie–Gower model. Existence of equilibria and their local asymptotic stabilities are carried out. Furthermore, with the help of bifurcation theory and center manifold theorem, existence and directions of period-doubling and Neimark–Sacker bifurcations are investigated at positive steady-state. In order to control chaos and bifurcations, the Ott–Grebogi–Yorke (OGY) method and the hybrid control strategy are introduced. Numerical simulations are also provided to illustrate the theoretical discussions.


Introduction
Discrete-time models administered by distinction conditions are much more suitable than the continuous ones when the population is of non-overlapping generations [1][2][3][4]. Due to large applications in nearly all fields of applied sciences, discrete dynamical systems are an effective point of research. Discrete dynamical systems show rich and complex dynamical properties, for example, transcritical bifurcation, flip bifurcation, Hopf bifurcation, and chaos [5][6][7][8]. To study the nonlinear systems and their ramifications, qualitative analysis of discrete dynamical systems is the main mathematical apparatus to use. Integrity of the discrete-time models adds significance in the research community. For an examination of the qualitative behavior of difference equations, an easy computational and graphical explanation tools of illustration are available, but the analytical analysis always keeps a tremendous task especially in the case of the complicated and amazingly varying difference equations [9][10][11]. The predator-prey co-operations are the foremost complex zones of the population environment. Their universal reality and extraordinary behavior have created the interest of scholars, ecologists, and mathematicians during the final few years [12]. An innovative result, i.e., Lotka-Volterra predator-prey model given by Lotka [13] and Volterra [14] independently, is the first and challenging mathematical model. The Lotka-Volterra predator-prey model has neglected many actual conditions and complications and was modeled in such a way that carries the reality. Leslie and Gower [15] suggested a predator-prey model, called the Leslie-Gower predator-prey model, in which the predator growth activity is clear-cut from the predator predation function. They presumed that the predator growth is interpreted by a function of the rate of predators and their prey. Many authors have employed this model to study the real global complications; for example, Wollkind, Logan, and Wollkind [16,17] applied this system to model the predator-prey speck outbreak synergies on fruit trees in Washington state. Aziz and Okiye [18] set up and examined the modified LG model. Xiao and Huang [19] examined a predator-prey system with Holling type-IV functional response. Pattern formation and bifurcation analysis were discussed for a class of predator-prey system in [20]. Influence of fear factor on predator-prey dynamics with prey refuge was analyzed in [21,22].
In biological mathematics, discrete-time models are used to examine the taxonomic group of organisms and species with the passage of time. Discrete-time population models are appropriate where the populations are non-overlapping and essentially remain constant over a generation. These models are best to describe the chaotic behavior of nonlinear dynamics [23,24]. Controlling chaos is an interesting topic in recent studies of nonlinear dynamical systems. Sometime bifurcation and chaotic behavior are highly unfavorable phenomena in dynamical systems because there may be an extinction of population due to chaos. So, controlling chaos by introducing new measures to the population is very important in the population dynamical systems. We are able to get chaos control by using different strategies, e.g., feedback control strategy, hybrid control technique, and poleplacement method. By using these techniques, one can delay, advance, or even eliminate the chaotic behavior due to emergence of bifurcation in dynamical systems.
Next, arguing as in [25], we consider the following class of predator-prey interaction with nonlinear harvesting in prey population and Holling type-IV functional response: where x and y are prey and predator density, respectively. All the parameters a, b, c, d, α, and β are positive constants. Moreover, α represents tolerance of the prey, β denotes half saturation constant, a is the catchability coefficient, b represents the effort for harvesting, c is intrinsic growth rate of the predator, and d is the amount of prey required to support one predator at equilibrium. The qualitative behavior of continuous system (1) has been investigated recently in [25]. In order to explore rich dynamical behavior of such predator-prey interaction, one may consider some discrete counterpart of system (1). For this, we apply Euler's approximation to system (1) as follows: y n+1 = y n + h cy n 1d y n x n . (2)

Stability analysis of equilibria
The equilibria for system (2) are the solutions of the following system: Then it follows that (x 1 , 0) and (x 2 , 0) are boundary equilibria of system (2), where x 1 and x 2 are real roots of the following quadratic equation: Moreover, assume that 0 < b < 1 and b < a < 1 4 (b + 1) 2 , then both x 1 and x 2 are positive with x 1 < x 2 . Next, we take x 1,2 =x, then the Jacobian matrix of system (2) at boundary equilibrium points (x, 0) is computed as follows: Then λ 1 = 1+h-2xh-abh (b+x) 2 and λ 2 = 1+ch are eigenvalues of J(x, 0) with λ 2 > 1. Therefore, both boundary equilibria (x 1 , 0) and (x 2 , 0) are unstable. Next, we are looking for interior equilibria of system (2). For this, the components of interior equilibrium (x * , y * ) are given as follows: where x * is a positive real root for the following quartic equation: In order to see the existence of positive equilibria graphically, we consider nullclines defined by y = x d and y = (1-x -a b+x )(β + x 2 α + x). Then Fig. 1 shows the existence of two distinct positive equilibria (0.219259, 0.609053) and (0.714336, 1.98427), and Fig. 2 reveals the existence of the unique positive equilibrium point (0.734318, 1.93242) for system (2).
Furthermore, the Jacobian matrix at positive equilibrium (x * , y * ) is given by where Furthermore, the characteristic equation for J(x * , y * ) is given as follows: Keeping in view the relation between roots and coefficients of a quadratic equation [26,27], we have the following lemma for the dynamics of interior equilibrium (x * , y * ).
(ii) (x * , y * ) is a saddle point if

Neimark-Sacker bifurcation
In this section, we study the existence of Neimark-Sacker bifurcation (NSB) for the steady state (x * , y * ) of system (2). We find the parametric conditions for a non-hyperbolic fixed point so that the Jacobian matrix has two complex conjugate eigenvalues with modulus 1. At fixed point E * = (x * , y * ), the characteristic polynomial equation can be written as where A = (1ch + w 11 ) and B = w 11 (1ch) -w 12 ch d .

Chaos control
In this segment, two feedback control techniques are talked about. The first one is the pole-placement methodology, which may be processed as a generalized OGY method, and the second method is known as the hybrid control strategy. In order to use the OGY technique to system (2), we have a tendency to write the system in the following way: y n , h), where h indicates the parameter for chaos control. Furthermore, h is limited to lie in some small interval |hh 0 | < with > 0, and h 0 indicates the nominal value belonging to the chaotic region. Suppose that (x * , y * ) indicates an unstable fixed point for (2) in the chaotic region, which yields under the influence of period-doubling and Neimark-Sacker bifurcation. In this case, (22) is approximated in the neighborhood of (x * , y * ) as follows: where . Moreover, if the following matrix has rank 2, then system (22) is controllable.
Since ch d w 13 = 0, the rank of C is 2. Consequently, system (2) is always controllable with the OGY feedback control method. Next we assume that [hh o ] = -K x n -x * y n -y * , where K = k 1 k 2 , then system (23) can be written as Moreover, both eigenvalues μ 1 and μ 2 of matrix A -BK lie in an open disk if and only if the equilibrium point (x * , y * ) is locally asymptotically stable. These eigenvalues are called regulator poles, and the problem of placing these regulator poles at appropriate place is known as a pole-placement problem. Since the rank of matrix C is 2, hence the poleplacement problem has a unique solution. Next we suppose that γ 2 + α 1 γ + α 2 is the characteristic equation of A and μ 2 + β 1 μ + β 2 is the characteristic equation of A -BK , then the unique solution of pole-placement problem is given by where T = CM and M = α 1 1 1 0 . Secondly, we apply the hybrid control strategy to system (2) for controlling chaos through the effect of both types of bifurcations. Assume that system (2) undergoes PDB and NSB at its equilibrium point (x * , y * ), then the corresponding controlled system using the hybrid control method is written as follows: where ρ ∈ (0, 1) illustrates the control parameter for controlled system (26). Moreover, the Jacobian matrix for (26) at (x * , y * ) is computed as follows: Furthermore, the following lemma gives parametric conditions for local stability of the fixed point (x * , y * ) for controlled system (26).

Numerical simulation and discussion
This segment is committed to proving the presence of PDB and NSB for system (2) for specific numerical values of its parameters (a, b, c, d, α, β) whereas taking h as bifurcation parameter. The verification of period-doubling and Neimark-Sacker bifurcation is demonstrated using bifurcation diagrams and phase portraits. In addition, hybrid control and pole-placement techniques for chaos control are too supported by numerical simulations.
Example 6.1 First we choose (a, b, c, d, α, β) = (0.4, 2, 5.4, 0.113, 0.012, 6) and the initial values (x 0 , y 0 ) = (0.75, 6.5). Then system (2) experiences PDB as bifurcation parameter h changes in the small neighborhood of h 0 = 0.3805. From bifurcation diagrams, which are shown in Fig. 3, for h < 0.3805, the fixed point E * is stable and unstable at h = 0.3805 and period-doubling bifurcation occurs for h > 0.3805. MLEs related to bifurcation diagrams Fig. 3(a) and (b) are plotted in Fig. 3(c), which confirms the existence of chaotic behavior in system (2). The phase portraits related to Figs. 3(a) and (b) are plotted in Fig. 4. On the other hand, for these parametric values λ 1 = -1, λ 2 = 0.7604011396 = ±1, l 1 = -11.13714613 = 0, and l 2 = 0.3269262762 = 0. Therefore, all the conditions of Theorem 3.1 are satisfied. Taking h = (h ok 1 (x nx * )k 2 (y ny * )), where K = (k 1 , k 2 ) is the gain matrix, the respective controlled system is given as x n y n x 2 n 0.19 + x n + 1 -0.001x n x n + 0.75 , By simple calculations, we get α 1 = -1.967321758, α 2 = 1.000352101, In this case the positive equilibrium of control system (29) is locally asymptotically stable for -69.03215478) < k < -0.4933532528. The bifurcation diagrams for controlled system (29) are shown in Fig. 7 by taking k 1 = 0.1 and -12 < k 2 < 4. Next we take k 1 = 0.1 and k 2 = -10 for controlled system (29), then the plots x n and y n are shown in Figs. 8(a) and (b) and the phase portrait is shown in Fig. 8(c).  (2) is unstable. In this case system (22) can be written as x n y n x 2 n 0.012 + x n + 6 -0.4x n x n + 2 = f (x n , y n , h), y n+1 = y n + h 5.4y n 1 -0.113 y n x n = g(x n , y n , h).
is the gain matrix, the respective controlled system is given as x n y n x 2 n 0.012 + x n + 6 -0.4x n x n + 2 , y n+1 = y n + h ok 1 x nx *k 2 y ny * 5.4y n 1 -0.113 y n x n .
The positive equilibrium of control system (32) is locally asymptotically stable for 0.030673 < k < 0.165560. The bifurcation diagrams for controlled system (32) are shown in Fig. 9 by taking k 1 = 0.1 and -0.02 < k 2 < 0.14. Next we take k 2 = 0.06 and k 1 = 0.1 for controlled system (32), then the plots x n and y n are shown in Figs. 10(a) and (b) and the phase portrait is shown in Fig. 10(c).
Under the Jury condition, the roots of (39) lie in a unit open disk if and only if 0 < ρ < 1. Thus NSB is entirely controlled in the largest admissible interval for ρ ∈ 0.9996888951. For ρ = 0.95, the plots and phase portrait of controlled system (37) are shown in Fig. 12.

Concluding remarks
We studied the qualitative behavior of a modified Leslie-Gower predator-prey model and achieved the results for stability of equilibrium points. In order to support the complexity in system (2), the presence of period-doubling and Neimark-Sacker bifurcation for the fixed point is verified mathematically, further numerically simulations are performed. Using these simulations, we showed that system (2) goes through PDB and NSB for the vast range of bifurcation parameters h. Chaos control is discussed through the implementation of pole-placement and hybrid feedback control methods. It is clear from our numerical observations that, for an extensive range of parameters, stability can be rebuilt through the OGY method and the hybrid control method. The OGY method is based on feedback control methodology, whereas hybrid control depends on feedback control and parameter perturbation. The computation of maximum Lyapunov exponents proved the presence of chaotic behavior of the system.