A Density-Dependent Host-parasitoid Model with Stability, Bifurcation and Chaos Control

: The aim of this article is to study the qualitative behavior of a host-parasitoid system with a Beverton–Holt growth function for a host population and Hassell–Varley framework. Furthermore, the existence and uniqueness of a positive ﬁxed point, permanence of solutions, local asymptotic stability of a positive ﬁxed point and its global stability are investigated. On the other hand, it is demonstrated that the model endures Hopf bifurcation about its positive steady-state when the growth rate of the consumer is selected as a bifurcation parameter. Bifurcating and chaotic behaviors are controlled through the implementation of chaos control strategies. In the end, all mathematical discussion, especially Hopf bifurcation, methods related to the control of chaos and global asymptotic stability for a positive steady-state, is supported with suitable numerical simulations. Formal analysis, Q.D.; Investigation, Q.D.; Resources, X.M. and M.R.; Data curation, N.J. and Y.F.; Writing–original preparation, Q.D. and M.R.; Writing–review and and Funding X.M.


Introduction
The host-parasitoid interaction is the most underlying and substantial procedure in population dynamics. Many species, such as semelparous animals and monocarpic plants, have discrete seasonal (non-overlapping) generations and their births take place in customary reproduction seasons. Their interactions are described by difference equations or formulated as discrete-time mappings.
Discrete-time models are widely used to explore the dynamics of host-parasitoid interactions in an ecosystem. In nature, most of host-parasitoid systems have complexity due to the interaction of species. On the other hand, in most of such cases, the coexistence of both species is mostly stable. We propose a host-parasitoid model with the implementation of a Beverton-Holt growth function for the host population and Hassell-Varley framework to make the coexistence stable. On the other hand, chaotic behavior is also observed for the proposed model. In the case of ecological models related to the breeding of biological species, chaos is considered as irregular and unpredictable behavior; therefore, the presence of chaos is an unfavorable situation for the survival of a species [1]. Thus, chaos control methods can be used by the species to avoid risks of irregular behavior. Qualitative analysis for various classes of host-parasitoid interaction is a topic of great interest. Recently, many researchers have focused their studies on the dynamics of host-parasitoid models. Some of these investigations, which are very close to the present study, are summarized as follows. Taylor [2] investigated the dynamics of some host-parasitoid systems under competition. Kaitala et al. [3] discussed a comprehensive study of the complex dynamics occurring in a basic discrete-time model of host-parasitoid interaction. Tang and Chen [4] applied Holling type II and III functional response functions to a host-parasitoid model. Xu and Boyce [5] reported that mutual interference in host-parasitoid interaction not only can stabilize the dynamics but may strongly destabilize as well. Lv and Zhao [6] studied the complex dynamics in a discrete-time model of host-parasitoid interaction based on a lower bound for the host. Din [7] investigated global dynamics for a class of host-parasitoid systems related to plant-herbivore interaction with strong predator functional response. Din [8] modified the host-parasitoid interaction with implementation of constant refuge effects and reported global dynamics for the proposed model. In [9], the authors modified the host-parasitoid interaction with a Pennycuick growth function for the host population and studied Neimark-Sacker bifurcation and chaos control. In [10], global stability and Hopf bifurcation were carried out for a class of host-parasitoid interaction. Din [11] reported Neimark-Sacker bifurcation and chaos control in a Hassell-Varley model. Global dynamics and Neimark-Sacker bifurcation were investigated for Beddington and generalized Beddington models in [12] and [13], respectively. Global stability, bifurcation analysis and control for host-parasitoid interactions were reported in [14] and [15]. Wu and Zhao [16] studied the qualitative behavior of a discrete host-parasitoid model with the host subjected to refuge and strong Allee effects. Jamieson [17] discussed global dynamics of the May host-parasitoid interaction. Liu et al. [18] explored a complex behavior and bifurcation analysis for a class of host-parasitoid interaction with an application of the Allee effect and Holling type III functional response. Moreover, in [19], chaos control and bifurcation analysis were studied for a host-parasitoid model with a lower bound for the host. In [20] and [21], the influence of a refuge effect was explored for certain classes of host-parasitoid models. In [22], the authors numerically investigated a system of partial differential equations that describe the interactions between populations of predators and prey. Suryanto et al. [23] considered a model of predator-prey interaction at fractional-order with a ratio-dependent functional response. Zhang et al. [24] reported dynamics of a predator-prey system with the weak Allee effect.
Generally, the host-parasitoid interaction is considered as non-overlapping generations. Therefore, difference equations can be used to describe their mathematical framework. One of the pioneer works for mathematical modeling of host-parasitoid interaction was presented by Nicholson and Bailey [25]. The general framework for mathematical modeling of a host-parasitoid system is presented as follows: where H n and P n denote population densities for host and parasitoid, respectively, at the n-th generation; f (H n , P n ) represents the fraction of host population that does not parasitize; r represents the number of eggs that are laid by a host; and c is used for the intrinsic growth rate of a parasitoid population.
Nicholson and Bailey assume that f (H n , P n ) = exp(−aP n ), which is derived from the zero-th term of the Poisson distribution, where a represents per capita searching efficiency of the parasitoid population. With this particular implementation of f (H n , P n ) = exp(−aP n ), the Nicholson and Bailey model is given as follows: It is worthwhile to note that coexistence is unstable in the case of the Nicholson-Bailey model (NBM). On the other hand, in our universe, most of the host-parasitoid interactions are stable. Therefore, many modifications are suggested to NBM to make its resemblance to natural host-parasitoid systems. For this, it is more appropriate to interchange the constant reproduction rate of hosts to some density-dependent function for the host population. Consequently, the implementation of such density-dependent factors can modify NBM to make its resemblance with the natural host-parasitoid interaction.
Keeping in view the single species density-dependent model for the host population, the Ricker model [26] is a classic discrete population model, which gives the expected number H n+1 of individuals in generation n + 1 as a function of the number of individuals in the previous generation. This model is given as follows: where s is interpreted as an intrinsic growth rate and l as the carrying capacity of the environment. The Ricker model was introduced in 1954 by Ricker in the context of stock and recruitment in fisheries. The model can be used to predict the number of fish that will be present in a fishery. Subsequent work has been derived from the model under other assumptions, such as scramble competition, within-year resource-limited competition or even as the outcome of source-sink Malthusian patches linked by density-dependent dispersal. The Ricker model is a limiting case of the Hassell model [27], which takes the form: where r ≥ 1 denotes the growth rate of the population; b > 0 is used to represent intra-specific competition; and k > 0 is a scaling parameter for describing the steady population size. If we take b = 1 in Equation (1), then the Hassell model is simply the Beverton-Holt model [28] given as follows: Biological assumptions related to Equation (2) are given as follows [29]: (i) It is assumed that the first encounter between the parasitoid and host is random. Moreover, one viable egg is laid by a parasitoid on a single host, which is killed by the parasitoid's progeny.
(ii) Keeping in view the law of mass action, the number of encounters H e of resources (hosts) with consumers (parasitoids) in generation n are proportional to the product of hosts and parasitoids present densities, and consequently, one has: H e := aH n P n where a is a positive constant representing the searching efficiency of the parasitoid, and P n , H n denote the population densities of the parasitoid and host, respectively, in generation n. (iii) The next generation of parasitoids is produced due to infection of hosts in the present generation. (iv) The hosts that are not infected produce their own offspring.
Keeping in view the assumptions (i) − (iv) and the Beverton-Holt model in Equation (2), we have the following system for host-parasitoid interaction: where c is used for the intrinsic growth rate of the parasitoid population. Due to the randomness of encounters between hosts and parasitoids, it is more appropriate to represent the probability of N encounters by a distribution, which depends on the average number of encounters per unit time.
Consequently, this condition leads to implementation of the Poisson distribution, and probability mass function is described as follows: where P represents the number of parasitoids, N denotes the number of encounters per unit time, and aP is used for the mean of the distribution, which represents the average number of encounters per unit time. Therefore, the portion of the hosts without parasitism are given as follows: Furthermore, it follows that P n+1 = cH e = cH n (1 − P(0)) = cH n (1 − exp(−aP n )) .
Next, from Equation (4), it follows that Then, from Equations (3), (4) and (5), it follows that: Arguing as in [30], one may consider the influence of parasitoid interference on the interaction of the host-parasitoid model. For this, the searching efficiency a of parasitoids can be taken as a = −qP −m n , where m is a constant for mutual interference, and q represents the quest constant, which denotes the searching efficiency of the parasitoid as P n = 1. Keeping in view the Hassell-Varley modification of the host-parasitoid interaction, System (6) is modified as follows: The remaining discussion of this paper is summarized as follows. In Section 2, we discuss the permanence of solutions of System (7). In Section 3, the existence, uniqueness and local stability analysis of the positive fixed point of System (7) are studied. Section 4 is dedicated to examining the parametric conditions for global asymptotic stability of the positive fixed point of System (7). Neimark-Sacker bifurcation about the positive fixed point of System (7) is studied in Section 5 with the implementation of the theory of normal forms. In Section 6, the pole-placement method and hybrid control strategy are implemented for controlling the chaotic and fluctuating behavior of System (7) about its positive fixed point. At the end, numerical simulations are presented to authenticate the mathematical analysis in Section 7. Moreover, theoretical findings are validated through experimental and field data based on statistical analysis of previous literatures.

Permanence
In the case of difference equations, the permanence of solutions is an important tool to discuss further dynamical behavior of the equations. Definition 1. Suppose {(H n , P n )} denotes a positive solution for System (7). Then, we say that System (7) is permanent if the following hold true: where m and M are some positive constants. Theorem 1. Assume that {(H n , P n )} is an arbitrary solution of System (7). Moreover, suppose that exp qM 1−m 2 < r and M m 2 < qcm 1 , then the following results hold true: Proof. Keeping in view the first equation of Model (7), it follows that H n+1 ≤ rH n 1+kH n , then the comparison argument yields that lim n→∞ sup H n ≤ r − 1 k = M 1 for every n = 1, 2, · · · . Next, the second equation of System (7) gives that P n+1 ≤ cH n ; therefore, it is easy to see that lim n→∞ sup P n ≤ c(r − 1) k = M 2 for every n = 1, 2, · · · . Next, focusing on the first equation of System (7), one has that Next, one can implement the transformation S n = 1/H n in Inequality (9) to obtain S n+1 ≤ α + βS n for all n = 0, 1, 2, · · · , where β = e qM 1−m 2 /r and α = k/r. If we suppose that S 0 = 1/H 0 and β < 1, then due to the comparison argument, one has m 1 := inf H n for every n = 1, 2, · · · , and it follows from the fact that lim n→∞ inf(1/S n ) = lim n→∞ sup S n . Next, we again focus on the second equation of System (7) and see that P n+1 ≥ qcm 1 P n M m 2 +qP n for every n = 0, 1, 2, · · · . On the other hand, we consider a transformation of the form P n = 1/Q n , then one has Q n+1 ≤ α 1 + β 1 Q n for every n = 0, 1, 2, · · · , where β 1 = M m 2 /qcm 1 and α 1 = 1/cm 1 . Furthermore, let P 0 = 1/Q 0 and β 1 < 1, then it is quite simple to see that m 2 = (1 − β 1 )/α 1 ≤ lim n→∞ inf P n for all n = 1, 2, · · · . Consequently, one has the desired results as follows: Next, one can apply Theorem 1 along with mathematical induction to prove the following result: (7).

Existence of Positive Fixed Point and Stability Analysis
The equilibria of System (7) satisfy Obviously, System (7) has a trivial equilibrium (0, 0), the second is a boundary equilibrium r−1 k , 0 and a positive equilibrium (H * , P * ), which satisfies 1 = r exp(q(P * ) 1−m )+kH * and P * = cH * 1 − exp −q(P * ) 1−m . The positive equilibrium (H * , P * ) cannot be found in closed form. The following result shows the existence and uniqueness of a positive equilibrium point (H * , P * ) ∈ 0, r−1 (7).
Proof. In order to show the existence and uniqueness of positive equilibrium, we construct a real-valued function F : 0, r−1 k → R defined as follows: Due to the fact r > 1 and f (0) = 1 q ln(r) 1 1−m > 0, one has On the other hand, one has f r−1 k = 0, and lim that 0 < m < 1, then for any ξ ∈ 0, r−1 k one has that r − kξ > 1 and Consequently, one has In order to see the existence of a positive fixed point geometrically, it is easy to see that ordinate of a positive fixed point (H * , P * ) satisfies the following transcendental equation: The intersection of Υ 1 (P) and Υ 2 (P) is depicted in Figure 1. Suppose that (H, P) be any equilibrium point of System (7), then the Jacobian matrix of System (7) evaluated at equilibrium point (H, P) is given as follows: Furthermore, about interior fixed point (H * , (7), the Jacobian matrix J(H, P) given in Equation (10) becomes: The following Lemma gives the conditions for local stability of a positive fixed point. (7) is a sink if the following inequalities are fulfilled: Proof. It is easy to compute the characteristic equation for Jacobian matrix J(H * , P * ) as follows where Arguing as in [31], the roots of Equation (11) lie inside the open unit disk if P(1) > 0, P(−1) > 0 and P(0) < 1, that is, and As Inequality (13) is satisfied automatically, the interior fixed point of System (7) is a sink and, thus, asymptotically stable if the Inequalities (12) and (14)

Global Stability
The global stability of the interior fixed point of System (7) is explored in this section. Arguing as in [32], first of all, the following Lemma is considered.

Lemma 2.
Taking into account some real intervals I = [a, b] and J = [c, d], and assume that f : I × J → I and g : I × J → J are real-valued continuous functions defined on rectangle I × J. Considering the following system: where n = 0, 1, 2, · · · and initial conditions (x 0 , y 0 ) ∈ I × J. Meanwhile, we assume that the following conditions are satisfied: (i) The function f (x, y) is non-decreasing with respect to x, and non-increasing with respect to y.
(ii) The function g(x, y) is non-decreasing in both variables x and y.

Hopf Bifurcation
This section is dedicated to the existence and direction of Hopf bifurcation about an interior fixed point of System (7).
The roots of Equation (11) are complex conjugate numbers if the following inequality holds true: Suppose that Inequality (24) holds true, then System (7) experiences Hopf bifurcation about its interior fixed point whenever det J(H * , P * ) c=c 1 = 1, where Next, we define a Hopf bifurcation curve S NB for System (7) as follows: Assume that (r, q, m, k, c 1 ) ∈ S NB , and consider Suppose thatc indicates a very small perturbation in c 1 such that |c| 1, and interchanging c 1 into c 1 +c in Equation (26), we have a perturbed map of the following form: Further, we suppose that c 1 +c > 1 and 0 < m < 1, then positive fixed point of Equation (27) is given by (H * , P * ) ∈ 0, r−1 k × 0, (c 1 +c)(r−1) k . Focusing on the following translations From Equations (27) and (28), it follows that: where Taking into account Equation (29), then the characteristic equation for its variational matrix about (0, 0) is of the following form: where The roots of Equation (30) are given as follows: On the other hand, we have Assume thatc = 0, then we have A(0) = exp(q(P * ) 1−m ) r + qc 1 (1−m)H * exp(−q(P * ) 1−m ) (P * ) m > 0 and A(0) < 2 because (r, q, m, k, c 1 ) ∈ S NB . Consequently, one has A(0) = ±2, 0, −1 gives λ n 1 , λ n 2 = 1 at c = 0 for all n = 1, 2, 3, 4. Consequently, the roots of Equation (30) do not meet the unit circle atc = 0. Next, we suppose that α = A(0) 2 and β = 1 2 4B(0) − A 2 (0) and considering a transformation of the following form: From Equations (29) and (31), one has where f (u, v) = a 13 a 12 x 2 + a 14 a 12 xy + a 15 a 12 x = a 12 u and y = (α − a 11 )u − βv. Next, we consider the following first Lyapunov exponent: Keeping in view the above calculations and normal form theory of bifurcation [33][34][35][36], we have the following result: Theorem 5. Assuming that L = 0, then Model (7) experiences Hopf bifurcation around its interior fixed point denoted by (H * , P * ) whenever c deviates in small locality of In addition, if L < 0, then a closed invariant curve of attracting nature bifurcates from the positive fixed point towards c > c 1 , on the other hand, if L > 0, then a repelling invariant closed curve bifurcates from the fixed point towards c < c 1 .

Chaos and Bifurcation Control
Control of bifurcation and chaos is a topic of great interest. Particularly, chaos and bifurcation control methods are more applicable for models of discrete nature because such 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 [37]. For further details related to chaos control methods in discrete-time systems, we refer to [38][39][40][41][42][43][44][45][46][47][48].
In this section, we discuss two chaos control methods for System (7). The first method is a modification of the Ott-Grebogi-Yorke (OGY) method [49], known as the pole-placement method [50]. For the application of the pole-placement method to Model (7), we rewrite this system as follows: where c is used for the sake of the control parameter. On the other hand, it is assumed that control parameter c satisfies |c − c 0 | < δ, where δ > 0 and c 0 represents some nominal value, which is located in the chaotic or bifurcating region. Next, we suppose that (H * , P * ) is an interior unstable fixed point of System (7). Furthermore, it is also assumed that (H * , P * ) is located in some chaotic region. Our main target is to move the unstable fixed point towards a stable one. For this, System (33) is linearized about the an unstable fixed point (H * , P * ) as follows: where Next, we define the following controllability matrix for System (33): Then, it is easy to see that rank matrix C is two. Indeed, H * 1 − exp −q(P * ) 1−m = 0; thus, the system is controllable. Furthermore, in order to apply the pole-placement method, we consider that [c − Consequently, System (34) takes the following form: On the other hand, fixed point (H * , P * ) is asymptotically stable provided that the multipliers for the variational matrix A − BK lie inside the open unit disk. Assume that µ 1 and µ 2 are multipliers of variational matrix A − BK. In addition, we take λ 2 + α 1 λ + α 2 = 0 as a characteristic polynomial for matrix A and µ 2 + β 1 µ + β 2 = 0 is taken as a characteristic polynomial for A − BK. Then, a unique solution of the pole-placement method is provided as follows: where T = CM and M = α 1 1 1 0 . Due to some simple computation, one has here, all these partial derivatives are calculated at (H * , P * , c 0 ). Next, from Equations (37) and (38), it follows that: Secondly, the hybrid control method [51] is applied to System (7). For this, the controlled system related to System (7) is written in the following form: Then, a variational matrix for System (39) about a positive fixed point (H * , P * ) is given as follows: For controllability of System (39) about its positive fixed point, the following result is presented: The positive steady-state (H * , P * ) for System (39) is controllable, if the following inequalities are satisfied:

Numerical Simulations and Discussion
This section is dedicated to the numerical verification of our theoretical discussion. For plots related to phase portraits, bifurcation diagrams and maximum Lyapunov exponents, Mathematica and Matlab are used.
It is easy to see that the roots of Equation (41) are given by λ 1,2 = 0.5539551792258897 ± 0.8325464908392879i such that |λ 1,2 | = 1. Consequently, (r, q, m, k, c) = (5.4, 0.9, 0.4, 8.3, 32.49) ∈ S NB . Moreover, some phase portraits for System (7) are depicted in Figure 4.   Example 2. Now, we discuss the affectivity of chaos control methods. First, we implement the pole-placement method. For this, let us consider that r = 12.5, q = 0.6, m = 0.5, k = 15.5, and 80 ≤ c ≤ 130, then Model (7) goes through Hopf bifurcation when c ≈ 93.7 (cf. Figure 5d). The diagrams for bifurcation of System (7) are shown in Figure 5 (see Figure 5a,b). Furthermore, Figure 5c reveals that an invariant closed curve vanishes at c = 125 and, consequently, a chaotic orbit appears. Next, choosing (r, q, m, k, c) = (12.5, 0.6, 0. 5, 15.5, 125). For these selected parameters, System (7) possesses a positive fixed point (0.13491760154988774, 15.244461821003831), which is unstable because, for these parametric values, the characteristic polynomial for variational matrix of Equation (7) about (H * , P * ) = (0.13491760154988774, 15.244461821003831) is of the following form: Then, one can easily compute the roots of Equation (42), which are given by λ 1,2 = 0.478597471 ± 0.92194061i and further calculation yields that |λ 1,2 | = 1.03876 > 1. Here, we select c 0 = 125 as the nominal value in the chaotic region for the following controlled system: where (H * , P * ) is positive fixed point of Equation ( With some simple calculation, one can easily find characteristic polynomial for matrix A − BK in the following: Calculations performing with Mathematica yields that the roots of Equation ( (43) are depicted in Figure 6. Secondly, we apply the hybrid control method. For this, let us consider r = 12.5, q = 0.6, m = 0.5, k = 15.5, and c = 93.7, then System (7) possesses a positive fixed point which is given by (0.17236, 14.50693) and characteristic polynomial for variational matrix of (7) about fixed point (0.17236, 14.50693) is computed as follows: A simple calculation yields that the roots of Equation ( It is easy to see that System (46) possesses a positive fixed point (H * , P * ) = (0.17236, 14.50693). On the other hand, the variational matrix for Equation (46) around (0.17236, 14.50693) is computed as follows: Furthermore, one can simply find characteristic polynomial for Equation (47) as follows: An application of Jury condition gives that roots of Equation (48) lie inside the open unit disk iff 0 < α < 0.999823. Therefore, Hopf bifurcation is almost eliminated completely. Next, choosing α = 0.9995 plots for System (46)    Moreover, for these selected values, the characteristic equation of the System (7) about fixed point (H * , P * ) = (0.3816089836378311, 0.7367710807002495) is given as follows: The real distinct roots of characteristic Equation (49) are given by: Obviously, |λ 1 | < |λ 2 | < 1. Consequently, the positive steady-state for System (7) is a sink. On the other hand, in order to discuss the global asymptotic stability of the system about this fixed point, one can easily verify Condition (16) of Lemma 3 as follows: Condition (50) shows that the unique positive equilibrium point (H * , P * ) = (0.38161, 0.73677) is globally asymptotically stable. In order to verify this fact numerically, we choose r = 4.75, q = 0.99, m = 0.95, k = 5.5 and c = 3.1 with initial conditions (H 0 , P 0 ) = (8 × 10 10 , 5 × 10 10 ) , which are far away from the neighborhood of the equilibrium point (H * , P * ) = (0.38161, 0.73677), but Figure 8 shows that the equilibrium point (H * , P * ) = (0.38161, 0.73677) is globally asymptotically stable.

Concluding Remarks
This paper is concerned with the investigation of some qualitative aspects of a host-parasitoid model. The model is a modification of a generalized Hassell-Varley model of a host-parasitoid system, which can be useful for the description of ecosystems. In addition, the host-parasitoid interactions in which the growth of the host is governed by the Beverton-Holt law have been studied. These findings reveal that the quest constant and mutual interference of parasitoid interaction can be sources for stabilization or destabilization factors for this interaction. On the other hand, restrained interference can assist coexistence between the parasitoid and the host. Meanwhile, this study reveals that the reproductive rate of the parasitoid and the proceeding population growth can also produce a strong destabilization factor, resulting in a variety of complexities and chaotic behavior. Our mathematical findings include the persistence of solutions with the implementation of a method of comparison. Moreover, the existence and uniqueness of an interior fixed point along its local stability analysis are carried out. It is proved that an interior fixed point of the system is globally asymptotically stable. On the other hand, the host-parasitoid model undergoes Neimark-Sacker bifurcation around its interior fixed point. Two chaos control methods are applied for the control of bifurcating and the chaotic behavior of the system. Our investigation reveals that the pole-placement method and hybrid control strategy both are effective for the stability of corresponding controlled systems. In order to validate the proposed generalized Hassell-Varley Model (7) with real data based on statistical analysis using both field and experimental data, the parametric values for System (7) are estimated in Table 1. Keeping in view the estimated parametric value in Table 1, it is easy to see that the positive equilibrium of System (7) is a sink for 20 ≤ c < 54.128. Moreover, System (7) undergoes Neimark-Sacker bifurcation at its positive steady-state as bifurcation parameter c passes through a critical value c ≈ 54.128. The bifurcation diagrams and phase portrait of System (7) are depicted in Figure 9. Consequently, our theoretical investigations show an excellent validation of field and experimental data.

Future Problems
It is interesting to work on bifurcation analysis of the models presented in [52][53][54][55]. Funding: This research is supported by the funds from Higher Education Commission Pakistan."