Bifurcation analysis and chaos control for a plant–herbivore model with weak predator functional response

Abstract The interaction between plants and herbivores is one of the most fundamental processes in ecology. Discrete-time models are frequently used for describing the dynamics of plants and herbivores interaction with non-overlapping generations, such that a new generation replaces the old at regular time intervals. Keeping in view the interaction of the apple twig borer and the grape vine, the qualitative behaviour of a discrete-time plant–herbivore model is investigated with weak predator functional response. The topological classification of equilibria is investigated. It is proved that the boundary equilibrium undergoes transcritical bifurcation, whereas unique positive steady-state of discrete-time plant–herbivore model undergoes Neimark–Sacker bifurcation. Numerical simulation is provided to strengthen our theoretical discussion.


Introduction
In mathematical biology, plant-herbivore models are basically modifications of preypredator models [8]. The interaction between plants and herbivores has been investigated by many researchers both in differential and difference equations. Kartal [31] investigated the dynamical behaviour of a plant-herbivore model including both differential and difference equations. In [30] Kang et al. discussed bistability, bifurcation and chaos control in a discrete-time plant-herbivore model. In [38] Liu et al. investigated stability, limit cycle, Neimark-Sacker bifurcations and homoclinic bifurcation for a plant-herbivore model with toxin-determined functional response. Li et al. [37] discussed period-doubling and Neimark-Sacker bifurcations for a plant-herbivore model incorporating plant toxicity in the functional response of plant-herbivore interactions. Similarly, for some other discussions related to qualitative behaviour of plant-herbivore models, the interested reader is referred to [2,10,21,[25][26][27]29,33,45,47,52] and references are therein.
Taking into account the interaction of the apple twig borer and the grape vine [3], a mathematical model for plant-herbivore interaction based on weak predator functional response is given by where x n and y n are population densities for grape vine and apple twig borer, respectively. Moreover, α, β and γ are positive parameters. For some earlier investigations related to model (1), we refer to [4][5][6][7]34,35]. Moreover, taking into account strong predator functional response, Din [12] discussed global stability of system (1), and Khan et al. [32] investigated Neimark-Sacker bifurcation for system (1) with strong predator functional response.
The main findings of this paper are summarized as follows: • Existence of equilibria and their stability analysis is investigated.
• It is proved that model (1) undergoes transcritical bifurcation at its boundary equilibrium by implementing bifurcation theory of normal forms and centre manifold theorem. • Direction and existence criteria for Neimark-Sacker bifurcation are investigated at positive steady-state of plant-herbivore model (1). • Pole-placement and hybrid control methods are implemented in order to discuss chaos control in system (1).
Moreover, in Section 2 the existence of steady-states and their local asymptotic behaviour is analysed. Section 3 is dedicated to investigate transcritical bifurcation about boundary steady-state of system (1). In Section 4, we discuss Neimark-Sacker bifurcation at positive steady-state of model (1). We study pole-placement chaos control and hybrid control strategies in Section 5. Finally, in Section 6, numerical simulations are carried out to authenticate the theoretical discussion.

Existence of equilibria and stability
In this section, first we investigate fixed points for system (1). For this, these fixed points must solve the following algebraic system: Simple computation yields the following steady-states for the plant-herbivore model (1): Assume that 0 < α < 1, β > 0, β/(1 − α + β) < γ < 1, then P * is unique positive equilibrium point of the system (1). For β = 0.2, the region in αγ -plane where inequalities 0 < α < 1, β/(1 − α + β) < γ < 1 are satisfied is depicted in Figure 1 by red shaded region. In order to see the dynamical behaviour of system (1) at equilibria, we first compute Jacobian matrix of system (1) at (x, y) as follows: Then, at (x, y) = (0, 0) the Jacobian matrix J(x, y) given in (3) reduces to Then topological classification of P 0 is given as follows: • P 0 is a sink if and only if α > 1 and 0 < γ < 1.
• P 0 is a non-hyperbolic point if and only if α = 1, or γ = 1.

Transcritical bifurcation
In this section, we investigate that boundary equilibrium P 1 undergoes transcritical bifurcation. For this, we assume that Consider the set TB given by Assume that (α, β, γ 0 ) ∈ TB , then system (1) is equivalently represented by the following 2-dimensional map: whereγ is very small perturbation in γ 0 . Furthermore, if we consider x = u − (1 − α)/β and y = v, then map (5) is transformed into the following map: An application of Taylor series expansion about (x, y,γ ) = (0, 0, 0) yields that where Since at γ 0 = β/(1 − α + β) the linear part of map (7) is already in canonical form. Therefore, in order to implement the centre manifold theorem [9], let W C (0, 0, 0) denotes the centre manifold for the map (7) which is evaluated at (0, 0) in a small neighbourhood of γ = 0. Then, W C (0, 0, 0) is computed as follows: Moreover, we define the following mapping which is restricted to the centre manifold W C (0, 0, 0): where Furthermore, it follows that Then, the following theorem gives the parametric conditions for existence and direction of transcritical bifurcation for system (1) at its boundary fixed point P 1 .

Chaos control
Chaos control is a method of stabilization with the help of small perturbations which are applied to unstable periodic orbits for a given system. The purpose of chaos control is to make chaotic behaviour more predictable and stable. For this, a small perturbation is applied as compare to the original size of the system under consideration, and in a result the natural dynamics of the system is prevented from any major modification. Recently, chaos control for irregular complex dynamics has developed as one of the leading topics in nonlinear applied science. The pioneer articles related to chaos control were published in 1990 and after that their number has been steadily increasing [46]. Moreover, the practical methods related to chaos control can be implemented in various areas such as communications, physics laboratories, biochemistry, turbulence, and cardiology [40]. Furthermore, in the case of dynamical systems related to biological breeding of species, chaos control is considered to be an important feature for investigation. On the other hand, population models with non-overlapping generations have more irregular complex behaviour. For some recent investigation related to chaos control in discrete-time models we refer to [1,[13][14][15][16][17][18][19][20][22][23][24] and references are therein.
In this section, first we discuss pole-placement chaos control method based on state feedback control which was introduced by Romeiras et al. [44] (also see [41]). This method is also known as generalized or modified OGY method proposed by Ott et al. [42]. For the application of pole-placement technique to model (1), one can rewrite this system as follows: where γ denotes parameter for chaos control. Furthermore, it is assumed that γ lies in a small interval of type |γ − γ 0 | < δ, where δ > 0 and γ 0 represents the nominal parameter lies in the chaotic region. One can apply pole-placement method in order to move unstable trajectory towards desired stable orbit. For this, suppose that (x * , y * ) be an unstable fixed point for model (1) lies in a chaotic region under the influence of transcritical or Neimark-Sacker bifurcations. Keeping in view the method of linearization, one can approximate system (16) in the neighbourhood of the unstable fixed point (x * , y * ) as follows: where .
In order to check that system (16) is controllable, the following matrix is computed: Assume that (x * , y * ) is positive fixed point, then it follows that y * (1 + x * ) = 0. Therefore, rank of C is 2 at positive equilibrium point. Next, we assume that [γ − γ 0 ] = −K x n −x * y n −y * , where K = [ k 1 k 2 ], then system (17) can be written as Moreover, equilibrium point (x * , y * ) is a sink if and only if both eigenvalues μ 1 and μ 2 for the matrix A−BK lie in an open unit disk. Moreover, μ 1 and μ 2 are called regulator poles. Placing these eigenvalues at appropriate position is called pole-placement method. Furthermore, pole-placement problem has unique solution because rank of matrix C is two. Denote λ 2 + α 1 λ + α 2 = 0 as characteristic equation for matrix A and μ 2 + β 1 μ + β 2 = 0 is characteristic equation of A−BK, then unique solution of the pole-placement problem is given as follows: where T = CM and M = α 1 1 1 0 . Then it follows that where all partial derivatives in (21) are evaluated at (x * , y * , γ 0 ). From (20) and (21), we have the following unique solution of pole-placement problem: Secondly, we apply another comparatively simple chaos control method known as hybrid control strategy based on parameter perturbation and state feedback control [39,50]. An application of hybrid control strategy to model (1) yields the following controlled system: where 0 < κ < 1 is control parameter. Assume that (x * , y * ) be an equilibrium point of system (22), then Jacobian matrix for system (22) computed at (x * , y * ) is given as follows: It is easy to observe that system (22) is controllable as long as the steady-state (x * , y * ) of system (22) is locally asymptotically stable. Keeping in view this fact, one has the following result for controllability of system (22).

Numerical simulation and discussion
Then, variational matrix for system (24) is computed as follows: 0.973684 −0.041407 1.86849 − 2.07035k 1 1 − 2.07035k 2 .   The triangular stability region bounded by these marginal stability lines is depicted in Figure 8. Next, for same parametric values we apply hybrid control method to system (1). In this case system (22) is written as follows: Then, Jacobian matrix of system (25) and its characteristic equation are computed as follows: Keeping in view Jury condition and aforementioned characteristic equation, |2 − 0.0263158κ| < 2 − 0.0263158κ + 0.0773684κ 2 < 2 if and only if 0 < κ < 0.340136. At κ = 0.33 the plots for system (25) are depicted in Figure 9.

Concluding remarks
Local asymptotic stability, bifurcation analysis and chaos control are investigated for the apple twig borer and the grape vine type plant-herbivore model. The discrete-time model is obtained by taking into account the weak predator functional response. Furthermore, existence of equilibria and their stability analysis is investigated. Taking into account the bifurcating behaviour for biological populations, it must be noted that these phenomena are essential for competition between plants and herbivores [11]. It is proved that plant-herbivore model undergoes transcritical bifurcation at its boundary equilibrium by implementing bifurcation theory of normal forms and centre manifold theorem. Direction and existence criteria for Neimark-Sacker bifurcation are investigated at positive steady-state of plant-herbivore model. Bifurcating behaviour and chaos have always been considered as disadvantageous phenomena in biology. On the other hand, due to chaotic behaviour population can undergo a higher risk of extinction due to the unpredictability. Therefore, these are catastrophic for the breeding of biological population [51]. In order to prevent biological populations from this ruinous situation, one might think about applications of chaos and bifurcation control. Pole-placement and hybrid control methods are implemented in order to discuss chaos control in the system. Furthermore, our investigation reveals that the pole-placement method is more effective as compared to the hybrid control method.

Future work
It must be noted that the weak predator functional response in system (1) is similar to Holling type-III function 1. Therefore, it is more interesting to apply Holling type-II functional response. With such implementation of Holling type-II functional response, plant-herbivore model (1) can be rewritten as follows: x n+1 = x n α(1 + y n ) + βx n , y n+1 = γ y n (1 + x n ).
Stability, bifurcation analysis and chaos control for model (26) is our future work for investigation.