A CLASS OF DISCRETE PREDATOR–PREY INTERACTION WITH BIFURCATION ANALYSIS AND CHAOS CONTROL

The interaction between prey and predator is well-known within natural ecosystems. Due to their multifariousness and strong link population dynamics, predators contain distinct features of ecological communities. Keeping in view the Nicholson-Bailey framework for host-parasitoid interaction, a discrete-time predator–prey system is formulated and studied with implementation of type-II functional response and logistic prey growth in form of the Beverton-Holt map. Persistence of solutions and existence of equilibria are discussed. Moreover, stability analysis of equilibria is carried out for predator–prey model. With implementation of bifurcation theory of normal forms and center manifold theorem, it is proved that system undergoes transcritical bifurcation around its boundary equilibrium. On the other hand, if growth rate of consumers is taken as bifurcation parameter, then system undergoes Neimark-Sacker bifurcation around its positive equilibrium point. Methods of chaos control are introduced to avoid the populations from unpredictable behavior. Numerical simulation is provided to strengthen our theoretical discussion. Mathematics Subject Classification. 39A30, 40A05, 92D25, 92C50. Received March 10, 2020. Accepted September 11, 2020.

x n+1 = x n G (x n , y n ) F (x n G (x n , y n )) , y n+1 = bx n (1 − G (x n , y n )) , n = 0, 1, 2, . . . , (1.1) where x n and y n are population densities of prey and predator, respectively, at generation n; b is the conversion efficiency from prey to predator, and F (x) is Beverton-Holt growth function defined by where k > 0 is the carrying capacity and λ > 1 is intrinsic per-capita growth. Moreover, G (x, y) is type-II functional response given by where a denotes predation rate and h represents handling time of predators per prey. Putting the values of F and G in system (1.1), we obtain the following predator-prey system: x n+1 = λx n exp − ayn 1+hxn 1 + λ−1 k x n exp − ayn 1+hxn , y n+1 = bx n 1 − exp − ay n 1 + hx n .
(1.2) Conversion efficiency from prey to predator Dimensionless Furthermore, biological meanings and units of state variables and parametric values of system (1.2) are presented in Table 1. We summarize the remaining discussion of this paper as follows. Persistence of solutions for system (1.2) is discussed in Section 2. Section 3 is dedicated to existence of biologically feasible equilibria, and conditions for local asymptotic stability of these equilibria are also investigated. In Section 4, we show that boundary equilibrium of system (1.2) undergoes transcritical bifurcation whenever growth parameter b of predator population is taken as bifurcation parameter. In Section 5, it is proved that system (1.2) undergoes Neimark-Sacker bifurcation around its interior equilibrium point. OGY, hybrid control and exponential-type control methods are introduced in Section 6. Lastly, numerical simulations are provided in Section ?? to illustrate our theoretical discussion.

Boundedness and persistence
In this section, we show that every solution of system (1.2) persists under certain parametric conditions. For this, the following Theorem is presented.
Proof. See Appendix A.
With mathematical induction one can prove the following result.

Existence of equilibria and stability analysis
In this section, first we explore the biologically feasible fixed points of system (1.2). The steady-states of system (1.2) satisfy the following two-dimensional algebraic equations:  Then, it is quite easy to see that O = (0, 0) and E = (k, 0) are trivial and boundary equilibria of (1.2), respectively. We are interested in interior equilibrium, which cannot be investigated in closed form. Neglecting the trivial solution of system (3.1), we left with the following algebraic system for coexistence: From first equation of system (3.2), it follows that exp Solving (3.3) for y gives Similarly, from second equation of system (3.2) and (3.3), it follows that Comparison of (3.4) and (3.5) with simplification yields that: From (3.6), it follows that x * is a positive real root for the following transcendental equation: where Simple calculation yields that Ψ 2 (0) = ln(λ) > 0, Ψ 2 (k) = 0 and Ψ 1 (0) = Ψ 1 (k) = 0. On the other hand, Ψ 2 (x) = − λ−1 kλ−λx+x < 0 and Ψ 1 (x) < 0 for all x ∈ (0, k). Therefore, Ψ 2 (x) is monotonically decreasing and Ψ 1 (x) is concave down in (0, k). Hence, transcendental equation (3.7) has a unique solution x * satisfying 0 < x * < k.
For a = 1.1, b = 1.5, λ = 2.9, h = 0.65 and k = 19.3, the intersection of Ψ 1 (x) and Ψ 2 (x) is depicted in Figure 1. Now, the following Lemma gives the dynamics of trivial and boundary equilibria of system (1.2).  Proof. (i) First, it is easy to see that the Jacobian matrix of (1.2) at trivial equilibrium O is given by: Then, eigenvalues of J(O) are µ 1 = λ > 1 and µ 2 = 0 < 1. Therefore, O is a saddle point.
(ii) Secondly, it is easy to see that the Jacobian matrix of (1.2) at boundary equilibrium E is given by: Then, eigenvalues of J(O) are µ 1 = 1 λ < 1 and µ 2 = abk hk+1 . Therefore, it follows that E = (k, 0) is a sink if and only if abk hk+1 < 1, it is a saddle point if and only if abk hk+1 > 1 and E is a non-hyperbolic equilibrium point if and only if abk hk+1 = 1.
Moreover, system (1.2) undergoes transcritical bifurcation at its boundary equilibrium E = (k, 0) whenever the parameter b varies in a small neighborhood of b ≡ b 0 = hk+1 ak . Classification of E with transcritical bifurcation curve (yellow curve) in ab-plane is depicted in Figure 2. Next, we see dynamics of system (1.2) at its positive steady-state, say, (x * , y * ), which satisfies the following algebraic system: The Jacobian matrix J(x * , y * ) of system (1.2) evaluated at (x * , y * ) is given by: Now characteristic polynomial of J(x * , y * ) is given by: The following Lemma gives dynamical behavior of system (1.2) at its positive fixed point.
Next, it is interesting to see the dynamical classification of system (1.2) in x * y * -plane about its positive fixed point. For such classification, various parametric values are selected for system (1.2). Moreover, regions for sink, source and saddle cases are represented by blue, red and green regions, respectively. For a = 3.3, b = 2.4, h = 4.2 and λ = 9.8 dynamical classification is depicted in Figure 3. On the other hand, for comparatively smaller parametric values, that is, a = 0.36, b = 0.5, h = 0.42 and λ = 1.4, the dynamical classification is depicted in Figure 4.

Transcritical bifurcation
In this section, we apply bifurcation theory of normal forms and center manifold theorem for emergence of transcritical bifurcation around boundary equilibrium of system (1.2). For this, first we assume that Moreover, we consider the following curve Next, we suppose that (a, b, k, h, λ) ∈ C T B , and writing system (1.2) in the following equivalent map: whereb is very small perturbation in b 0 . Moreover, taking into account the translations p = u − k and q = v, then map (4.1) is converted into the following map with fixed point at origin: , Considering the following transformation From (4.2) and (4.3), we obtain Due to center manifold theory [10], stability analysis of equilibrium (x, y) = (0, 0) nearb = 0 can be discussed by investigating a one-parameter family of reduced equations on a center manifold W C (0, 0, 0), which can be described as follows: On the other hand, the map restricted to center manifold is described as follows: Then, it is easy to see that The following Theorem gives conditions for emergence of transcritical bifurcation about boundary equilibrium of system (1.2).
Theorem 4.1. Assume that κ = 0 and b = hk+1 ak , then system (1.2) undergoes transcritical bifurcation around its boundary equilibrium when parameter b changes in small neighborhood of b 0 = hk+1 ak . Moreover, two fixed points bifurcate from boundary equilibrium for b < b 0 , and merge as the boundary equilibrium at b = b 0 and disappear at b > b 0 .

Chaos control
Control of chaos and bifurcation is a theme of remarkable interest. In particular, bifurcation and chaos control methods are more applicable for models of discrete nature because these system are of complex behavior as compare to continuous ones. On the other hand, chaos control methods are widely used in almost all branches of applied science and engineering.
This section is dedicated for implementation of chaos control methods to system (1.2). For this, three chaos control methods are implemented to system (1.2). Implementation of chaos control methods for discrete-time systems is topic of great interest. Recently, many dynamical systems have been discussed with implementation of strategies related to chaos and bifurcations control [1,11,12,40,42,44,54,60,61]. In [15], variety of chaos control methods have been implemented to a discrete class of predator-prey interaction. A generalized hybrid control method based on parameter perturbation and state feedback control has been proposed in [16]. A new chaos control strategy of exponential type has been introduced in [17] for controlling fluctuating and chaotic behavior of discrete-time models. We apply two frequently used and one recently developed chaos control methods to system (1.2) as follows.
First, we consider Ott-Grebogi-Yorke (OGY) method (see [41]). An application of OGY feedback control method to system (1.2) yields the following control system: where k 1 and k 2 are control parameters, and (x * , y * ) is interior equilibrium point of system (1.2). In order to see the controllability of system (6.1), the Jacobian matrix for this controlled system at (x * , y * ) is computed as follows: In addition, the characteristic polynomial for J C (x * , y * ) is given by: Taking into account the controllability of system (6.1), the following Lemma is presented.
Lemma 6.1. System (6.1) is controllable if the following condition is satisfied: Secondly, we apply hybrid control method [35] to system (1.2) to get the following control system: where 0 < α < 1 is control parameter for hybrid control method, which is based on parameter perturbation and state feedback control. Again controllability of system (6.2) depends upon stability of this system at its positive fixed point. To see this stability behavior, the Jacobian matrix J H (x * , y * ) of system (6.2) at its unique positive fixed point (x * , y * ) is given as follows: In addition, the characteristic polynomial for J H (x * , y * ) is computed as follows: Taking into account the controllability of system (6.2), the following Lemma is presented.
Lemma 6.2. Positive fixed point (x * , y * ) of system (6.2) is a sink if the following condition is satisfied: Thirdly, our aim is to apply recently proposed exponential chaos control technique to system (1.2) [17]. With application of exponential type chaos control strategy, we obtain the following control system: where s 1 and s 2 are control parameters for exponential control strategy and (x * , y * ) be interior (positive) fixed point of system (1.2). Next, controllability of system (6.3) is directly related to stability of system about its positive fixed point (x * , y * ). For this, Jacobian matrix J E (x * , y * ) of system (6.3) about (x * , y * ) is computed as follows: Moreover, the characteristic polynomial for J E (x * , y * ) is computed as follows: Taking into account the controllability of system (6.3), the following Lemma is presented.
Note that variations of parametric values affect control performance of system (6.3). For example, if we choose a = 11.5, b = 0.68, h = 0.9, λ = 1.9 and k = 3, then controllability (stability) region (blue region) of system (6.3) is depicted in Figure 5.  amplification of Figure 6 is depicted in Figure 7. On the other hand, for b = 29.5, 29.575, 29.58, 30 phase portraits of system (1.2) are depicted in Figure 8.

Numerical simulation and discussion
Next, we check the validity of chaos control strategies discussed in Section 6. For this, we choose (λ, a, h, k, b) = (10.4, 0.5, 5.8, 1.2, 40). At these parametric values system (1.2) has unique positive fixed point (x * , y * ) = (0.360595, 12.5197). Clearly, this fixed point is source because multipliers of system (1.2) are µ 1 = 1.01709 + 0.815335i and µ 2 = 1.01709 − 0.815335i with |µ 1 | = |µ 2 | = 1.30355 > 1. First, we apply OGY method for these parametric values. In order to see the controllability of system (6.1), the Jacobian matrix for this controlled system at (0.360595, 12.5197) reduces to: The characteristic polynomial for J C (0.360595, 12.5197) is given by An application of Jury condition gives that system (6.1) is controllable if |0.312993k 2 − 2.03419| < 2.69925 − 0.0132964k 1 − 0.540299k 2 < 2. Furthermore, the controllable region in k 1 k 2 -plane with OGY method is shown in Figure 9. It must be noticed that the stability region depicted in Figure 9 is a region bounded by the following It is worthwhile to mention that the marginal stability lines L 1 , L 2 and L 3 are consequences of P (1) = 0, P (−1) = 0 and P (0) = 1, respectively. Secondly, we see the effectiveness of hybrid control method. For this, similar parametric values are chosen, that is, (λ, a, h, k, b) = (10.4, 0.5, 5.8, 1.2, 40) are used in system (6.2). Then, simple calculation yields that system (6.2) is controllable if −0.0514052 < α < 0, which is impossible because 0 < α < 1. Therefore, we are unable to control system (1.2) at these parametric values with hybrid control method. The failure of hybrid control method is due to the choice of bifurcation parameter b = 40, which located at most right end of chaotic region [29.575, 45]. On the other hand, if we choose b = 35, then system (6.2) is controllable if 0 < α < 0.257381. Thus    The characteristic polynomial for J E (0.360595, 12.5197) is given by P (ϕ) = ϕ 2 + (0.361s 1 + 12.52s 2 − 2.034)ϕ + 1.699 − 21.6s 2 + s 1 (4.5s 2 − 0.11) . An application of Jury condition gives that system (6.3) is controllable if |0.361s 1 + 12.52s 2 − 2.034| < 2.699 − 21.6s 2 + s 1 (4.5s 2 − 0.11) < 2. Moreover, the controllable region in s 1 s 2 -plane with exponential-type method is shown in Figure 10. It is easy to see that in Figure 10 two red regions are depicted as stability (controllable) regions for system (6.3). Consequently, at λ = 10.4, a = 0.5, h = 5.8, k = 1.2 and b = 40, OGY and exponential-type control methods are partially successful for controlling fluctuating and bifurcating behavior due to appearance of Neimark-Sacker bifurcation. On the other hand, hybrid control method is unable to avoid the population densities from unpredictable situation at extreme right part of chaotic region.

Concluding remarks
This paper is concerned to investigation of some qualitative aspects of a discrete-time prey-predator model. The model is a modification of classical Nicholson-Bailey framework related to host-parasitoid interaction with implementation of Beverton-Holt growth function for prey population and type-II functional response. Persistence of solutions is discussed with implementation of method of comparison. Moreover, the existence of fixed points along their stability analysis are carried out. It is proved that system undergoes transcritical bifurcation around its boundary fixed point. On the other hand, discrete prey-predator model undergoes Neimark-Sacker bifurcation around its interior fixed point. Three chaos control methods are applied for control of bifurcating and chaotic behavior of system. Our investigation reveals that OGY and exponential-type control methods are more effective as compare to hybrid control strategy for stability of corresponding controlled systems. Furthermore, in order to validate our theoretical findings with respect to observed field data related to predator-prey interaction, the following parametric estimations are presented: Taking into account the estimated parametric values from observed field data in Table 3, we have abk 1+hk = 0.033916 < 1. Consequently, boundary fixed point (k, 0) is a sink for these parametric values. Furthermore, if we vary bifurcation parameter b and keep all other parametric values fixed (that is, k = 30, λ = 1.309, a = 2.207 and h = 2.96), then boundary equilibrium is a sink if 0 < b < 1.35629, and it undergoes transcritical bifurcation at b ≈ 1.35629. On the other hand, positive equilibrium point is a sink if 1.35629 < b < 1.51164. At b ≈ 1.51164, the positive fixed point (x * , y * ) = (10.8055, 2.69622) undergoes Neimark-Sacker bifurcation. This whole scenario is depicted in Figure 11. Consequently, conversion efficiency from prey to predator plays vital role for altering the dynamics of predator-prey interaction. Arguing as in [13], the conversion efficiency may increase due to change in temperature (warming) and in a result simulated population dynamics reveal a destabilisation of predator-prey interaction with warming. On the other hand, conversion efficiency must be high when the actual prey biomass is low and must be low when actual prey biomass is high. Furthermore, the impact of changes in edibility and the conversion efficiency on the effective prey biomass must outweigh the impact of changes in actual prey biomass [52]. Then, solution of (A.1) is given by: where c 1 depends upon initial value u 0 of difference equation (A.1). Since λ > 1, therefore obviously one has lim n→∞ u n = k.
Suppose that x 0 = u 0 , then by comparison we have lim n→∞ sup x n ≤ k for all n ≥ 0. It is easy to see that second equation of system (1.2) gives y n+1 ≤ bx n , and thus it follows that lim n→∞ sup y n ≤ bk. Again, we consider first equation of the system (1.2) as follows: Assume that λ > e abk , then again comparison argument yields that lim n→∞ inf x n ≥ m 1 , where Keeping in view the inequality x 1+x ≤ (1 − e −x ) and considering second equation of the system (1.2), we have y n+1 ≥ abm 1 y n 1 + hk + ay n .