1. Introduction
The host–parasitoid model in discrete dynamical systems is a way to describe the population dynamics of a host species and its parasitoid species over discrete-time intervals. This model is based on a set of equations that describe how the populations of the two species change from one time step to the next. The model assumes that the populations of the host and parasitoid species are discrete, meaning they are represented by integers. A model that has secured appreciable recognition from both theoretical and experimental point of view is the host–parasitoid model. The host–parasitoid collaborations are the most accepted topics in mathematical biology as they are necessary to deal with the natural enemy of a pest or insect [
1,
2,
3,
4,
5,
6]. We observe that parasitoids are insect species whose larvae grow as parasites on further insect species. It has been discovered that larvae of parasitoids generally kill and sometimes paralyze its host. However, usually, mature parasitoid individuals are free-living insects. Moreover, parasitoids are characterized, in general, by insects that show one or more larval stages that parasite other arthropods, developing inside them and killing them before the end of their life cycle. Life cycles of parasitoids and their hosts are usually synched, e.g., hosts and parasitoids are univoltine (one generation per year) species. Discrete-time steps corresponding to generations are usually used in such models [
7,
8,
9,
10]. Parasitoids commonly include species of wasps, flies, beetles, and worms. On the other hand, a special case is where some preys are fully safe from attack of predators within a progressive or spatial refuge. Spatial refuges can take a variety of forms lying between two extremes, which are given as follows:
- (i)
Where the constant ratio of the host populations is safe within refuge.
- (ii)
Where the constant number is protected.
The most common among these are the constant proportion refuges. For example, flour moth caterpillars are protected from attack of parasitoids by the ichneumonid (
Nemeritis canescen) when they are sufficiently deep in the flour medium to be out of reach of the parasitoid’s ovipositor. In this way, only a proportion of the host’s habitat is searched [
11,
12,
13,
14,
15]. Before starting our discussion related to the mathematical modeling of host–parasitoid interaction under refuge effects, it is worthwhile to describe the pioneering contribution proposed by Nicholson and Bailey [
16]. The mathematical model presented by Nicholson and Bailey is acknowledged as the Nicholson–Bailey (NB) model and is directed by the autonomous planar system of nonlinear difference equations as follows:
where
and
are state variables representing population concentrations of host and parasitoid at some instant
, respectively. Moreover,
and
are growth rates of host and parasitoid, respectively, and
denotes the searching efficiency of parasitoids. Further investigation reveals that the NB model has trivial and interior (positive) fixed points. Moreover, a unique positive fixed equilibrium point occurs when
and is always unstable, which contradicts the fact that the NB model is a noble descriptive of regular host–parasitoid interactions because most of these interactions are stable and do not possess any oscillatory behavior. Various attempts were made by numerous mathematicians and ecologists to formulate appropriate modifications to the classical NB model. One of these modifications is to introduce the refuge effect on the classical NB model. For this, we study two modifications proposed by Hassell in [
17]:
Firstly, we consider a constant population of hosts, keeping in view the NB model, to examine the spatial refuge effect. A fraction of host specie
is accessible to the parasitoid specie in every single generation. Consequently,
is the population within spatial refuge. This is represented by the following model.
Relative proportions lead to stability under refuges but large areas of parameter space remain for unstable equilibrium (too few hosts within the refuge) or no equilibrium at all (too many hosts within the refuge).
Secondly, consider constant number refuge, where
hosts are always protected, and we obtain:
The spatial refuge can play an important role in the dynamics of the model, potentially leading to the persistence of the host population even in the presence of high parasitoid attack rates.
The main contribution of this article is to explore the dynamics of systems (2) and (3). It would be beneficial first to clarify certain fundamental concepts and mathematical outcomes pertaining to our primary inquiry, prior to delving into an all-encompassing examination of the qualitative behaviors exhibited by these systems. Since systems (2) and (3) are planar discrete-time models governed by autonomous nonlinear systems of difference equations, therefore we start with the dynamical behavior of the general planar systems. The key investigations of the manuscript are synopsized in the following discussion.
2. Stability Analysis around Equilibria
In order to obtain the equilibrium points of the system (2), we transform the system as follows.
The equilibria
and
are obtained, where
,
. To demonstrate the stability (local) analysis of a biologically suitable equilibrium, the following lemma is the best interpretation [
18,
19].
Lemma 1. Let This is called the characteristic equation of variational matrix and and are the eigenvalues of this equation. Now, a unique positive equilibrium point is termed as a sink if the absolute values of and is one. Consequently, a sink is always a local asymptotic stable equilibrium point. When along with , the equilibrium point is known as a repeller; thus, we always have a source which is unstable. When and or and , then the equilibrium point is called a saddle and it is always unstable. When the absolute value of either of and is one, then we get a non-hyperbolic equilibrium point.
Taking into account the positivity and existence of model (2) at
, one can take
. In
Figure 1, the blue region shows the sink, and the red region shows the source. Moreover, the yellow region shows the Hopf bifurcation. The nonexistence of the model (2) is represented by the white area. Thus, the Jacobian matrix at
is:
and the characteristic polynomial is given as:
By simple calculations, we obtain:
and
Now, suppose be two roots of the system of Equations (2), then the positive equilibrium point is given by According to the lemma 1, system (2) has the following conditions related to the stability analysis.
The equilibrium point is a sink, if and only if and .
The equilibrium point is a source, if and only if and .
The equilibrium point is a saddle point, if and only if .
Now, we consider the dynamics of model (3) by replacing and with and . It is clear that system (3) has two equilibrium points, the trivial equilibrium point (0, 0) and the unique positive equilibrium point , with and
The Jacobian matrix computed at
is:
Clearly, is a sink if and only if and a saddle if .
The Jacobian matrix at the interior steady state
is calculated as:
and the characteristic polynomial is given by:
By simple calculations, we obtain:
From the above mathematical computations, we deduced the following results:
The unique steady-state
is locally asymptotically stable if and only if:
The unique steady-state
is unstable if and only if:
The unique steady-state
is saddle if and only if:
The eigenvalues of the characteristic polynomial are complex conjugate with magnitude 1, if and only if:
3. Neimark–Sacker Bifurcation in Model (2)
Now, we study the presence and trajectory of the Neimark–Sacker bifurcation for the unique positive equilibrium point
of model (2). We must find sufficient conditions for the Neimark–Sacker bifurcation. The interior unique positive equilibrium point
should be non-hyperbolic for the existence of Neimark–Sacker bifurcation, such that the variational matrix calculated at this point has two complex conjugate eigenvalues with an absolute value equal to one. Discrete-time models corresponding to population can also be investigated it the same way [
20,
21,
22]. The conditions that determine the limitations for the existence of the Neimark-Sacker bifurcation are as follows:
or
and
Suppose that
, then there is change of parameters in a small neighborhood of
, where
produces the Neimark–Sacker bifurcation for the positive equilibrium point of system (2). The following two-dimensional map expresses the system of Equation (2) as:
Let
indicate the bifurcating parameter, and then the perturbed mapping of the above system is as follows:
where
and the small perturbation parameter is denoted by
Where
and
. Then from above map, we have
where
The characteristic equation for the Jacobian matrix at equilibrium point (0, 0) is
Since
. Hence, the complex conjugate roots can be written as
where
and
Now, to obtain the non-degeneracy conditions, we observe that
Moreover, we have
because
. Now, consider that
Hence, we observe that
and
, that is consider
Moreover,
Make sure that
, and as a result we have
,
when
. Henceforth, the roots of the system cannot exist within the intersection of the circle with a unit radius with the synchronized axes if
. To obtain the normal form when
, we assume that
at
Thus,
and
Furthermore, considering the following transformation:
Under this transformation, our system can be written as:
Then consider a real number as follow:
where
After conducting the aforementioned detailed investigation, we have the following precise result:
Theorem 1. Assume that (1) holds and then system (2) undergoes the Neimark–Sacker bifurcation at the unique positive equilibrium point , when the parameter varies in a small neighborhood of Furthermore, if then an attracting invariant closed curve bifurcates from the equilibrium point for and if then a repelling invariant closed curve bifurcates from the equilibrium point for
4. Neimark–Sacker Bifurcation in Model (3)
The two-dimensional map for model (3) is defined as:
Let
represent the bifurcation parameter where
, and then corresponding perturbed mapping of above system is stated as follows:
where
and the small perturbation parameter is denoted by
Then afterwards, consider the transformations
x =
X −
and
y =
Y −
, the above map implies that
where
The characteristic equation for the Jacobian matrix at the equilibrium point (0, 0) is
Since
. Hence, the complex conjugate roots can be written as
where
Thus,
Now, for obtaining the non-degeneracy conditions, we have
Moreover, we have
because
. Also,
Now, suppose that
that is,
To obtain the normal form at
, we assume that
at
where
Furthermore, we consider the following transformation:
Under this transformation, our system can be written as
where
and
Now, we consider a real number as follows:
where
After conducting the detailed investigation described above, we have obtained the following precise result:
Theorem 2. Assume that map (5) holds and then system (3) undergoes the Neimark–Sacker bifurcation at the unique positive equilibrium point , when the parameter varies in a small neighborhood of Furthermore, if then an attracting invariant closed curve bifurcates from the equilibrium point for and if then a repelling invariant closed curve bifurcates from the equilibrium point for
5. Chaos Control
For population models, an important feature considered is controlling chaos and bifurcation. Often, we observe that discrete-time models show more complex behavior than continuous models. Chaos control techniques are used to escape population from unpredictable situations. Here, we use the hybrid control feedback methodology to reduce chaos arising because of bifurcation. We use this technique to reduce chaos emerging because of the Neimark–Sacker bifurcation. As we know that system (2) experiences the Neimark–Sacker bifurcation at the equilibrium point
, then a controlled system corresponding to this system can be expressed as [
23,
24,
25,
26]:
where
. Here, the growth rate is increased by
in both populations, the host and parasite, where as
is a harvesting term which is introduced in a correction to the equations. Moreover, the controlled tactic of system (5) is the fusion of both the feedback control and parameter perturbation. Furthermore, by appropriate selection of the controlled parameter
, the Neimark–Sacker bifurcation for the equilibrium point
of the controlled system (5) could be postponed or excluded entirely. Now, the variational matrix of controlled system (5) calculated at equilibrium point
is expressed as:
Furthermore, the characteristic equation of the Jacobian matrix of the controlled system is
The condition for equilibrium point of the controlled system to be locally asymptotically stable is
As we know that system (1) experiences the Neimark–Sacker bifurcation at a unique equilibrium point, then a controlled system corresponding to this system can be expressed as:
where
. Moreover, the controlled strategy of the above controlled system is coma bination of both the feedback control and the parameter perturbation. Furthermore, by an appropriate option for the controlled parameter
the Neimark–Sacker bifurcation for the equilibrium point
of controlled system can be postponed or excluded totally. The variational matrix of the controlled system calculated at equilibrium point is:
Additionally, the characteristic equation associated with variational matrix of the controlled system is as follow:
For the controlled system, the criteria for equilibrium point to be locally asymptotically stable are the following:
6. Numerical Simulations
Here, we discuss the qualitative and complex dynamical behavior of models (2) and (3) with numerical simulations by introducing some examples. In this section, we elaborate our mathematical and theoretical study with the help of numerical simulations by selecting particular parameter values.
Example 1. Let us consider system (2). Assume that , and . Then, the system (2) experiences Neimark-Sacker bifurcation when the bifurcation parameter .
Furthermore, system (2) has unique positive steady-state which is represented by .
For the parametric values the characteristic polynomial becomes Then, the roots of , are and with Therefore, , then the corresponding bifurcation diagram for model (2) can be shown in the following figures:
Figure 2 shows that after the bifurcation point
, the populations described in the model (2) experience quasiperiodic behavior due to emergence of Neimark–Sacker bifurcation. Furthermore,
is the non-chaotic region for the model (2), and
is its chaotic region. From
Figure 2c, it is obvious to observe that the maximum Lyapunov exponents are almost about zero due to the quasiperiodic behavior.
Moreover, it is interesting to see phase portraits of the system (2) for various values of the bifurcation parameter in the chaotic region. For this, we keep fixed
,
and the initial conditions
, and choosing some values of
in the chaotic region, then the phase portraits of the system (2) are portrayed in
Figure 3 as follows:
Next, we want to observe the effectiveness of the hybrid control method in the chaotic region due to the emergence of the Hopf bifurcation in the model (2). For this, we choose two fixed parameters, that is,
,
, and keep on varying the other two parameters in the approximate intervals. For this, it is suitable to choose the whole chaotic region
for the bifurcation parameter
and the whole allowable unit interval
for the control parameter
. The stability region (controllability region) of system (6) can be shown by
Figure 4, which shows the effectiveness of hybrid chaos control methodology:
Example 2. Let ,
,
and ,
then system (3) undergoes the Neimark–
Sacker bifurcation when .
The bifurcation diagrams for model (3) can be shown in Figure 5:
The
Figure 5 shows that at
and
model (3) experiences the Neimark–Sacker bifurcation. Furthermore, at
a unique positive equilibrium point exists. The phase portraits of model (3) for
,
and
are depicted in
Figure 6a–d as:
7. Conclusions
The spatial refuge can have important implications for the population dynamics of host–parasitoid systems. The presence of a refuge can increase the survival of the host population, reduce the parasitoid attack rate, lead to increased competition among parasitoids, and make the system more complex and difficult to predict. Therefore, it is important to consider the effects of spatial structure when modeling or studying host–parasitoid systems.
In this manuscript, two host–parasitoid models are discussed. Both the models are modifications of the Nicholson–Bailey model proposed by Bailey et al., in [
16]. One of these modifications is to introduce the refuge effect on the classical NB model. For this, we study the two modifications proposed by Hassell. At first, some qualitative features of host–parasitoid models are discussed. These features include the existence of equilibrium points, the stability of equilibria, and the numerical simulations. Moreover, the mathematical results of [2, 3] have been reviewed. They investigated the outcomes of the existence of refuge on the local stability and bifurcation of the models on the three different host–parasitoid models with the spatial refuge effect. The Taylor series method has been used in proposed model (2). It is also proved that model (3) undergoes the Neimark–Sacker bifurcation around the interior equilibrium point. Moreover, it is observed that both proposed models undergo the Neimark–Sacker bifurcation, and the maximum Lyapunov exponents also show that the models experience the Neimark–Sacker bifurcation. In addition to this, the hybrid chaos control methodology is used to control the chaotic behavior of the proposed models. Numerical simulations are provided for further explanation.
Author Contributions
Conceptualization, M.S.S. and Q.D.; methodology, M.S.S.; software, M.S.S.; validation, M.S.S.; formal analysis, M.S.S. and Q.D.; investigation, M.S.S.; resources, M.S.S.; data curation, M.S.S.; writing—original draft preparation, M.S.S., W.F.A. and M.D.l.S. All authors have read and agreed to the published version of the manuscript.
Funding
This research is funded by the Basque Government under the Grants IT1555-22 and KK-2022/00090; and to MCIN/AEI 269.10.13039/501100011033 for Grant PID2021-358 1235430B-C21; and to MCIN/AEI 269.10.13039/501100011033 for Grant PID2021-358 1235430B-C22. The APC was funded by M. De la Sen.
Data Availability Statement
There is no data associated with the manuscript.
Acknowledgments
Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2023R 371), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia. The author M. De la Sen is grateful to the Basque Government for its support through Grants IT1555-22 and KK-2022/00090; and to MCIN/AEI 269.10.13039/501100011033 for Grant PID2021-358 1235430B-C21; and to MCIN/AEI 269.10.13039/501100011033 for Grant PID2021-358 1235430B-C22.
Conflicts of Interest
The authors declare no conflict of interest.
References
- Agarwal, T. Bifurcation analysis of a host-parasitoid model with the beddington-deangelis functional response. Glob. J. Math. Anal. 2014, 2, 565–569. [Google Scholar] [CrossRef]
- Bešo, E.; Mujic, N.; Kalabušić, S.; Pilav, E. Neimark–Sacker bifurcation and stability of a certain class of host-parasitoid models with host refuge effect. Int. J. Bifurc. Chaos 2020, 29, 1950169. [Google Scholar] [CrossRef]
- Bešo, E.; Mujic, T.; Kalabušić, S.; Pilav, E. Stability of a certain class of a host-parasitoid model with a spatial refuge effect. J. Biol. Dyn. 2020, 14, 1–31. [Google Scholar] [CrossRef]
- Wu, D.; Zhao, H. Global qualitative analysis of a discrete host-parasitoid model with refuge and strong Allee effects. Math. Methods Appl. Sci. 2018, 41, 2039–2062. [Google Scholar] [CrossRef]
- Din, Q.; Hussain, M. Controlling chaos and Neimark-Sacker bifurcation in a host-parasitoid model. Asian J. Control. 2019, 21, 1202–1215. [Google Scholar] [CrossRef]
- Din, Q.; Saeed, U. Bifurcation analysis and chaos control in a host-parasitoid model. Math. Methods Appl. Sci. 2017, 40, 5391–5406. [Google Scholar] [CrossRef]
- Kapçak, S.; Ufuktepe, U.; Elaydi, S. Stability and invariant manifolds of a generalized Beddington host-parasitoid model. J. Biol. Dyn. 2013, 7, 233–253. [Google Scholar] [CrossRef] [PubMed]
- Livadiotis, G.; Assas, L.; Dennis, B.; Elaydi, S.; Kwessi, E. A discrete-time host–parasitoid model with at Allee effect. J. Biol. Dyn. 2014, 9, 34–51. [Google Scholar] [CrossRef]
- Liu, X.; Chu, Y.; Liu, Y. Bifurcation and chaos in a host-parasitoid model with a lower bound for the host. Adv. Differ. Equ. 2018, 2018, 1–15. [Google Scholar] [CrossRef]
- Ma, X.; Din, Q.; Rafaqat, M.; Javaid, N.; Feng, Y. A Density-dependent host-parasitoid Model with stability, bifurcation and chaos control. Mathematics 2020, 8, 536. [Google Scholar] [CrossRef]
- Schreiber, S.J. Host-parasitoid dynamics of a generalized Thompson model. J. Math. Biol. 2006, 52, 719–732. [Google Scholar] [CrossRef] [PubMed]
- Sebastian, E.; Victor, P. Stability analysis of a host-parasitoid model with logistic growth using Allee effect. Appl. Math. Sci. 2015, 9, 3265–3273. [Google Scholar] [CrossRef]
- Wiggins, S. Introduction to Applied Nonlinear Dynamical Systems and Chaos; Springer: New York, NY, USA, 2003. [Google Scholar]
- Yuan, L.G.; Yang, Q.G. Bifurcation, invariant curve and hybrid control in a discrete-time predator-prey system. Appl. Math. Model. 2015, 39, 2345–2362. [Google Scholar] [CrossRef]
- Luo, X.S.; Chen, G.R.; Wang, B.H.; Fang, J.Q. Hybrid control of period-doubling bifurcation and chaos in discrete nonlinear dynamical systems. Chaos Solitons Fractals 2003, 18, 775–783. [Google Scholar] [CrossRef]
- Bailey, V.A.; Nicholson, J. The balance of animal populations. Proc. Zool. Soc. Lond 1935, 105, 551–598. [Google Scholar] [CrossRef]
- Hassell, M.P. The Dynamics of Arthropod Predator-Prey Systems; Pritnceton Univeristy Press: Princeton, NJ, USA, 1978. [Google Scholar]
- Jury, E.I. Inners and Stability of Dynamic Systems; Research supported by the National Science Foundation; Wiley-Interscience: New York, NY, USA, 1974; Volume 316. [Google Scholar]
- Liu, X.; Xiao, D. Complex dynamic behaviors of a discrete-time predator-prey system. Chaos Solitons Fractals 2007, 32, 80–94. [Google Scholar] [CrossRef]
- Shabbir, M.S.; Din, Q.; Alabdan, R.; Tassaddiq, A.; Ahmad, K. Dynamical complexity in a class of novel discrete-time predator-prey interaction with cannibalism. IEEE Access 2020, 8, 100226–100240. [Google Scholar] [CrossRef]
- Taylor, A.D. Parasitiod competition and the dynamics of host-parasitoid models. Am. Nat. 1988, 132, 417–436. [Google Scholar] [CrossRef]
- Tang, S.Y.; Chen, L.S. Chaos in functional response host-parasitoid ecosystem models. Chaos Solitons Fractals 2002, 13, 875–884. [Google Scholar] [CrossRef]
- Xu, C.L.; Boyce, M.S. Dynamic complexities in a mutual interference host-parasitoid model. Chaos Solitons Fractals 2005, 24, 175–182. [Google Scholar] [CrossRef]
- Jamieson, W.T. On the global behaviour of May’s host-parasitoid model. J. Differ. Equ. Appl. 2019, 25, 583–596. [Google Scholar] [CrossRef]
- Tassaddiq, A.; Shabbir, M.S.; Din, Q. Discretization, bifurcation and control for a class of predator-prey interaction. Fractal Fract. 2022, 6, 31. [Google Scholar] [CrossRef]
- Kaitala, V.; Ylikarjula, J.; Heino, M. Dynamic complexities in host-parasitoid interaction. J. Theor. Biol. 1999, 197, 331–341. [Google Scholar] [CrossRef] [PubMed]
| Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. |
© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).