Stability and bifurcation analysis of a discrete-time host-parasitoid model with Holling III functional response

1.   Introduction

Population model, epidemic model and complex network model are important components of biomathematics models, which have been extensively researched by scholars for many centuries. In particular, the stability of an ecological system that contains many species interacting via competition, predation and symbiosis, and subjects to environmental and demographic noise, poses an interesting mathematical and ecological puzzle. Among all these biomathematics models, host-parasitoid models, play an important role in population dynamics. The interactions between parasitoids and their hosts are of great interest to ecologists and mathematicians^[1,2,3,4].

Among the variety of models, the host-parasitoid models with different Holling type functional responses are refined so as to better reflect the specific characteristics of the different populations or economical needs. The functional response represents the intake rate of a species as a function of food density which is one of the primary component in defining the biological system's dynamics. A vast interesting works exhibiting quite rich and complicated dynamics have been performed by the theoretical and mathematical biologists considering complex Holling type functional responses ^[5,6,7,8]. For instance, Xiang et.al ^[5] investigate the Bogdanov-Takens bifurcation and Hopf bifurcation of a host-parasitoid model with Holling II functional response. Wang and Zhang ^[6] show the importance of pesticide application on the pest control in a host-parasitoid model with Holling II functional response. Liu et.al ^[7] analyze the flip and Neimark-Sacker bifurcation in a host-parasitoid model with Allee effect and Holling III functional response. Cabello et.al ^[8] explore the effect of functional response on a host-parasitoid model and give a biological interpretation for this entomophagous species model.

It is worth noting that the numerical simulation is actually based on the discrete form of continuous model and the discrete form is a natural connection between the actual model and the simulation. In fact, many species in real life have no overlap between successive generations and their populations evolve in discrete time-steps. The discrete model governed by a difference equation is more direct, convenient and accurate to formulate and simulate than the continuous one because of the fact that the population has a short life expectancy, non-overlapping generations in the real world. Thus the discrete models are been noticed more in recent years ^[9,10,11,12]. Motivated by the above mentioned discussions, here we consider a discrete-time host-parasitoid model with Holling type III functional response, which will be meaningful work. This model is ^[13]

where $x_{n}$ and $y_{n}$ stand for the densities of host and parasitoid in generation $n\; (n = 0, 1, 2, ...)$, respectively; $r$ denotes the intrinsic growth rate of the host population without parasitoid, $K$ represents the carrying capacity, $T$ is the total time initially available for search, i.e., the total time the hosts are exposed to parasitoids, $T_{h}$ is the handling time, i.e., the time between hosts being encountered and search being resumed, $bx_{n}/(1+cx_{n})$ stands for the instantaneous search rate which depends on the host density, and $b, \; c$ are given constants.

Tang and Chen have investigated the complexities of system (1.1) numerically. So, our main motivation is to algebraically show that the system will experience flip bifurcation and Neimark-Sacker bifurcation which will be derived by the center manifold theorem and bifurcation theory ^[14,15]. The detailed existence conditions of each bifurcation are given in a strict mathematical way. Numerical simulations are performed, including bifurcation diagrams, phase portraits and maximum Lyapunov exponents to validate the theoretical results and exhibit some new and interesting kinetics. The contents studied in this paper have never been considered in reference ^[13].

The remainder of the paper is organized as follows. In Section 2, we establish the boundedness conditions for the solution of model (1.1). In Section 3, we explore the uniqueness and local stability of the positive equilibrium. In Section 4, we show that there exist some values of parameters such that model (1.1) undergoes flip and Neimark-Sacker bifurcation. In Section 5, numerical simulation results are presented for supporting the theoretical findings and exhibiting new and rich dynamical regimes. The final section summarizes results and provides conclusions.

2.   Boundedness of positive equilibrium

Biologically speaking many species obey their own ecological genetic laws, and all have their own bounds. Once the population size is above the upper bound or below the lower bound, the species will explode or die out which can lead to an imbalance in the ecosystem. Thus, it is necessary to consider the upper or lower bound of a biological model.

The boundedness of model (1.1) can be found by the following lemma.

Lemma 2.1. (^[16]) Assume that $x_{0} > 0$ for every $x_{t}$ and $x_{t+1} < x_{t}\exp(A(1-Bx_{t}))$ for $t\in[t_{1}, \infty]$, where $B > 0$ is a positive constant. Then,

Using Lemma 2.1, we state the following theorem for the uniform boundedness of model (1.1).

Theorem 2.1. Any positive solution $(x_{n}, y_{n})$ of model (1.1) is uniformly bounded.

Proof. Suppose that $(x_{n}, y_{n})$ is arbitrary positive solution of model (1.1). Then, we have

Let $x_{0} > 0$. Using Lemma 2.1, we acquire

Through the solution of model (1.1), we obtain the following expression:

Let $y_{0} > 0$. We conclude

Thus, it follows that $\lim \limits_{n \rightarrow \infty}Sup (x_{n}, y_{n}) \leq l, \; l = \max\{l_{1}, l_{2}\}$. The proof is completed. □

3.   Existence and stability of positive equilibrium

By solving model (1.1), one can explore the existence of all equilibria. A simple algebraic computation shows that model (1.1) has two boundary fixed points $E_{1}(0, 0)\; and\; E_{2}(K, 0)$. In addition, we also discuss the existence and uniqueness of the positive solution of model (1.1) because the positive equilibria are not in a closed form. For this purpose, the following computation exhibits the existence and uniqueness at the positive solution of model (1.1).

Theorem 3.1. There exists an unique positive equilibrium $(x_{*}, y_{*})\in [0, l_{1}]\times[0, l_{2}]$ of model (1.1).

Proof. We begin our analysis of the model (1.1) by solving the equations

Suppose that

for all $x\in[0, l_{1}]$. Then, it follows that $g(0) = r > 0$. Assume that $\exp[r(1-x/K)] > 1$, we have

for all constants are positive. Thus, there exists at least one positive real root of $g(x) = 0$ in $[0, l_{1}]$. And because

Therefore, the model (1.1) has an unique positive equilibrium $(x_{*}, y_{*})\in [0, l_{1}]\times[0, l_{2}]$. And by direct computations to the model (3.1), it is easy to obtain that the unique equilibrium $E(x_{*}, y_{*})$ satisfies the following equations:

where $R = exp(r(1-x_{*}/K))$. □

Initially, we explore the stability of the boundary points $E_{1}(0, 0)\; and\; E_{2}(K, 0)$. The Jacobian matrix $J_{0}$ at $(0, 0)$ is given as

Accordingly, we can get eigenvalues $\lambda_{1} = e^{r} > 1, \; \lambda_{2} = 0 < 1$. From which we can easily check that $E_{1}$ is an unstable node. Also, two eigenvalues at $E_{2}(K, 0)$ are $\lambda_{1} = 1-r, \; \lambda_{2} = -\frac{bTK^{2}}{1+cK+bT_{h}K^{2}}$. Then, $E_{2}$ is stable when $|1-r| < 1$ and $|bTK^{2}/(1+cK+bT_{h}K^{2})| < 1$. Next, we just need to consider the stability of the equilibrium $E(x_{*}, y_{*})$.

The Jacobian matrix of model (1.1) evaluated at the unique positive equilibrium $E(x_{*}, y_{*})$ is given by

where

Moreover, the characteristic polynomial of $J$ reads as

where

In order to study the local stability and bifurcation phenomenon of the positive equilibrium, the following lemma will be very useful and essential, which can be easily proved by the relationship between roots and coefficients of a quadratic equation ^[17,18].

Lemma 3.1. Let $F(\lambda) = \lambda^{2}+P\lambda+Q$ and $F(1) > 0$. Suppose $\lambda_{1}$ and $\lambda_{2}$ are two roots of $F(\lambda) = 0$. Then, the following results hold true.

i) $|\lambda_{1}| < 1$ and $|\lambda_{2}| < 1 \; \Leftrightarrow \; F(-1) > 0$ and $Q < 1$;

ii) $|\lambda_{1}| > 1$ and $|\lambda_{2}| > 1 \; \Leftrightarrow \; F(-1) > 0$ and $Q > 1$;

iii) $|\lambda_{1}| < 1$ and $|\lambda_{2}| > 1$ (or $|\lambda_{1}| > 1$ and $|\lambda_{2}| < 1$) $\Leftrightarrow \; F(-1) < 0$;

iv) $\lambda_{1} = -1$ and $|\lambda_{2}|\neq1 \; \Leftrightarrow \; F(-1) = 0$ and $P\neq0, 2$;

v) $\lambda_{1}, \lambda_{2}$ are complex and $|\lambda_{1}| = |\lambda_{2}| = 1 \; \Leftrightarrow \; P^{2}-4Q < 0$ and $Q = 1$;

vi) $\lambda_{1} = \lambda_{2} = -1 \; \Leftrightarrow \; F(-1) = 0$ and $P = 2$.

By applying Lemma 3.1 and some definitions in ^[17,18], the following conclusions can be obtained.

Theorem 3.2. For the positive equilibrium $E(x_{*}, y_{*})$, we have the following estimates.

1) When $r < \frac{2+L+M+MN}{G+GMN}$,

a) if $r > \frac{M-1}{GMN}$, $E$ is a sink;

b) if $r < \frac{M-1}{GMN}$, $E$ is a source.

2) When $r > \frac{2+L+M+MN}{G+GMN}$, $E$ is a saddle.

3) $E$ is non-hyperbolic if and only if

or

From the above conclusion, if the first condition of (3) is true then one can easily find that one of the eigenvalues of $E$ is $-1$ and the other is $\lambda_{2} = 2-rG+L+MN$ which is neither 1 nor -1. If the second condition holds true, the eigenvalues of $E$ are a pair of complex conjugate numbers whose modulus is 1.

Consider

The unique positive equilibrium $E(x_{*}, y_{*})$ of model (1.1) undergoes flip bifurcation when parameters vary in a small neighborhood of $F_{B}$.

Let

The unique fixed point $E(x_{*}, y_{*})$ can pass through a Neimark-Sacker bifurcation when the parameters fluctuate in a small neighborhood of $H_{B}$.

4.   Bifurcation analysis

Based on the previous analysis, we choose the parameter $r$ as a bifurcation parameter to explore flip and Neimark-Sacker bifurcation for the point $E(x_{*}, y_{*})$ by using the center manifold theorem and bifurcation theory in ^{[14,15,17,18]}.

Suppose that all parameters $(r_{1}, b, c, K, T, T_{h})$ belong to $F_{B}$. If the parameter $r$ varies in a small neighborhood of $r_{1}$, then, the model (1.1) experiences a flip bifurcation at the fixed point $E$. For simplicity, model (1.1) can be written in the form:

Giving a perturbation $\overline{r}$ of parameter $r_{1}$, the model (4.1) is converted to

Let $H = x-x^{*}, P = y-y^{*}$, and transform the fixed point $E(x^{*}, y^{*})$ to the origin $O(0, 0)$ and model (4.2) into

where

The detailed expressions of the involved coefficients can be found in Appendix A.

The canonical form of (4.3) at $r_{1} = 0$ can be obtained by assuming the following map:

The normal form of model (4.3) under translation (4.4) can be expressed as

where the expressions $g_{1}(u, v, \overline{r})$ and $g_{1}(u, v, \overline{r})$ are provided in Appendix B.

From center manifold theory, we know that there exists a center manifold $W^{C}(0, 0, 0)$, which can be approximated as follows:

where

Consequently, the map restricted to center manifold $W^{C}(0, 0, 0)$ is expressed as follows:

where the detailed expressions of $h_{i}$ are given in Appendix B.

Straightforward but detailed calculations show that

which are the existence conditions for flip bifurcation.

Combining the above conclusions, and applying the bifurcation theory presented as Section 3.2 in Guckenheimer^[15], we state the following conclusion for flip bifurcation.

Theorem 4.1. If $\alpha_{2}\neq 0$, there exists a flip bifurcation at $E(x_{*}, y_{*})$ of model (1.1). Furthermore, the period-2 orbits bifurcated from $E(x_{*}, y_{*})$ are stable (unstable) if $\alpha_{2} > 0$ ($\alpha_{2} < 0$) which is also named supercritical (subcritical) flip bifurcation.

Next, we investigate the Neimark-Sacker bifurcation of model (1.1), which occurs when the eigenvalues of the characteristic equation are complex conjugates with unit modulus. i.e., $|\lambda_{1}| = 1, \; |\lambda_{2}| = 1$.

Now choose $r$ as a bifurcation parameter and let $H = x-x_{*}, P = y-y_{*}$. Here $r = r_{2}+\widetilde{r}$ and $|\widetilde{r}| < < 1$. The fixed point $E$ is shifted to the origin, the model (4.1) is converted into the new following form:

Where

and $a_{ij}, b_{i}, d_{i}$ are described in Appendix A by replacing $r = r_{2}+\widetilde{r}$. Let

be the characteristic equation of the Jacobian matrix of model (4.6) evaluated at $(0, 0)$ with coefficients

Since $(r_{2}, b, c, K, T, T_{h})\in H_{B}$, the roots of the characteristic equation are

Hence,

There are non degeneracy conditions $\lambda_{1, 2}^{j}\neq1, \; j = 1, 2, 3, 4$ when $\widetilde{r} = 0$ if and only if $m^{2}(\widetilde{r})-4n(\widetilde{r}) < 0$ and $m(\widetilde{r})\neq0, -1$, i.e.,

Now we construct the canonical form of (4.6) at $\widetilde{r} = 0$ by taking $\alpha = \frac{m(0)}{2}, \; \beta = \frac{\sqrt{4n(0)-m^{2}(0)}}{2}$ and assuming

Under transformation (4.7), the map (4.6) yields

where

and $H = a_{12}u, \; P = (\alpha-a_{11})u-\beta v$.

Next, we characterize a nonzero real number

where

Keeping in view the aforementioned calculations and implementing the bifurcation theory for normal forms. Ultimately, we deduce the following theorem for the direction and existence of the N-S bifurcation.

Theorem 4.2. There is a Neimark-Sacker bifurcation at $E(x_{*}, y_{*})$ whenever $r$ deviates within the small neighborhood of $r_{2} = \frac{M-1}{GMN}$. Furthermore, if $\theta < 0$ then the Neimark-Sacker bifurcation is supercritical and the closed invariant curve is stable. If $\theta > 0$, the Neimark-Sacker bifurcation is subcritical and the closed invariant curve is unstable. Also, it implies that both species in the system can coexist under some specific conditions.

5.   Numerical simulations

In this section, we conduct numerical results to validate the above theoretical analyses and to illustrate some complex dynamical behaviors of model (1.1). We take the bifurcation parameters in the following two cases into account.

Case 1: Take parameters $b = 0.0015, \; T = 100, \; K = 3, \; T_{h} = 1, \; c = 0.03$ in model (1.1), and vary $r$ in the range of $1.8 \leq r \leq 2.8$. By some complex calculations, we find that the flip bifurcation emerges from $E(2.8001, \; 0.3737)$ with $\alpha_{1} = -1.07576, \; \alpha_{2} = 0.31807 > 0$ when $r = 2.15$. Hence, the theoretical results in Theorem 4.1 are verified.

Figure 1(a) shows two-dimensional bifurcation diagram of model (1.1) in $(r, x)$ space when $r$ varies in $[1.8, 2.8]$. The corresponding Lyapunov exponents are shown in Figure 1(b) to confirm the existence of periodic orbits and chaotic motions. Figure 1(b) shows that some Lyapunov exponents are negative, which implies that there exists stable fixed point or stable periodic window; the others are positive, which means that there exist chaotic regions. From Figure 1, we see that the fixed point $E$ is stable for $r < 2.15$, loses its stability for $r > 2.15$ and a flip bifurcation takes place for $r = 2.15$. With the increase of $r$, period-doubling bifurcation leads to the birth of a chaos. Figure 1 also shows that the loss of stability through increasing $r$ yields drastic swings in host population size. In agricultural scenarios, these host outbreaks could have devastating consequences.

Case 2: Fix $b = 0.003, \; T = 100, \; K = 5, \; T_{h} = 1, \; c = 0.03$, let $r$ vary in the range ^[1,2]. We get an unique equilibrium $E(2.248523, \; 1.150799)$ when $r = 1.303$ and its eigenvalues are $\lambda_{1, 2} = 0.222866 \pm 0.974952i$. For $r = 1.303$, we have $|\lambda_{1, 2}| = 1, \; \theta = -0.36798753 < 0.$ It is clear that the Neimark-Sacker bifurcation occurs at this time, which confirms the effectiveness of theoretical analysis in Theorem 4.2.

Figure 2(a) shows two-dimensional bifurcation diagram of model (1.1) in $(r, x)$ space when $r$ varies in a neighbourhood of $r = 1.303$. The corresponding maximum Lyapunov exponent is shown in Figure 2(b) to confirm the existence of periodic orbits and invariant circle. More specially, the appearance of invariant circle demonstrates that these two populations can coexist in the same environment.

In above sections, we have investigated flip bifurcation and Neimark-Sacker bifurcations analytically and numerically. In fact, we find that model (1.1) is a complex multiparameter system. Therefore, it is natural to ask how the parameters affect the dynamics of the model, what may happen and what match rules about parameters are when two or more parameters change simultaneously, which are what we concern. Here we use some high-definition resolution phase diagrams to display more and richer nonlinear dynamics for this model, which are described by the efficient methods described in ^{[19,20,21,22]}.

Figure 3(a) represents period stability phase of the $(r, \; b)$ parameter-space with $K = 5, \; T = 100, \; T_{h} = 1, \; c = 0.03$. In which period$-1$ solution is plotted in brown, period$-29$ solution is plotted in light gray, the regions for period $2-28$ are codified in other different colors. In which, there exist periodic structures numbered as $1, 3, 4, 5, 6, ...$ embedded in the chaotic regions. Except these periodic solutions, there exist quasiperiodic solutions, which arise after a Neimark-Sacker bifurcation and divergence regions for which the related parameters are corresponding to unbounded orbits. Here, period motions are indicated by numbers, quasiperiodic regions are indicated by the letter Q, and chaotic regions are indicated by the letter C. The variations of complex parameter space reveal that small changes of key parameters could significantly influence the oscillation patterns of both host and parasitoid species. Consequently, this could result in difficulties for the development of biological control theory.

Figure 3(b) shows the magnification of Figure 3(a) for $1 < r < 2$ and $0.0025 < b < 0.01$. From which, we can see Arnold tongue ^[19,20] sequences $5, 10, 20, 6, 12, 24, ...$ due to the period-doubling bifurcation. Figure 3(b) also shows that the Arnold tongues have self-similar properties. Some Arnold tongues with fractal structure are embedded in quasi-periodic regions, that is, the mode-locking regions and the quasi-periodic regions alternately appear.

Unexpectedly, in Figure 3(b), we obtain that the mode-locking structure is formed depending on the broken Farey tree by examining the rotation number of these Arnold tongues, which orders all rations numbers between $0$ and $1$ depending on a strict definition ^[19]. We suppose that $p/q$ and $r/s$ are two adjacent tongues of rotations numbers there exists a tongue with rotation number $(p+r)/(q+s)$ between these two irreducible fractions. For instance, the tongue with the rotation number $3/13$ is situated in tongues with rotation numbers $2/9$ and $1/4$, tongue $2/9$ between the tongues $1/4$ and $1/5$ and the tongue $2/11$ between the tongues $1/5$ and $1/6$, which is only a small and incomplete part of a Farey tree, not an intact one, so we call it the broken Farey tree.

Next, we will investigate some details about the kinetics of model (1.1) from Figure 3(a). Figure 4 shows the bifurcation diagrams when $r = 1.5$ which is marked by the vertical line in Figure 3(a). By increasing $b$ from 0.0025, the system stays in a period-1 state until $b = 0.003$, being characterized by the maximum Lyapunov exponent (MLE) is less than zero. At $b = 0.003$, MLE reaches zero value and a Neimark-Sacker bifurcation occurs at this point. These two phenomena can be seen in Figure 5(a) and (b). By increasing $b$ until $b = 0.0042$, period-5 orbit appears in the system and flip bifurcation results in period-10 and period-20 orbits which are shown in Figure 5(c)–(e). When $b = 0.005522$, we cross the quasiperiodic region shown in Figure 5(f). By increasing $b$ more, we can find a period-5 orbit again at $b = 0.007578$ which can be seen in Figure 5(g). Finally a chaotic state whose MLE is greater than zero is created at $b = 0.008106$. A typical chaotic attractor appears in Figure 5(h) for $b = 0.008902$. In fact, this shows a classical route to chaos.

Another attractive phenomenon is the coexistence of multi-attractors which shows the host final states. It follows from the bifurcation analysis that multiple attractors can coexist for a wide range of parameters. To further confirm this and address their biological implications, we set all parameters as the same as in Figure 3. The overlapped colors in the neighborhood of two points $P1(1.438, \; 0.00586)$ and $P2(1.492, \; 0.007585)$ in the parameter plane of Figure 3(b) are a forceful implication of multistability. As we know through coloring diverse kinetic behaviors, attraction basins make the separation of attractors possible in kromograms. Figure 6(a) and (b) depict coexistence of period-5 and period-12 attractors, period-5 and period-7 attractors in this discrete model, respectively. Thus, it is shown that the model (1.1) is characterized by a higher degree of unpredictability and complexity of basin's structures.

6.   Conclusions

In summary, we discuss the comprehensive rich dynamics for a discrete host-parasitoid system with Holling type III functional response. The theoretical analysis demonstrates that the model (1.1) undergoes flip bifurcation and Neimark-Sacker bifurcation. In the case of these codimension one bifurcations, the normal forms for maps and their approximations to the flows of the corresponding differential equations are obtained to determine the bifurcation diagrams. All the numerical continuations are in concordance with analytical predictions. Furthermore, the high-quality phase diagrams are employed to illustrate the mode-locking structures in the quasiperiodic rotation motions. Finally, coexisting multi-chaotic attractors such as coexistence of seven-piece, twelve-piece and five-piece chaotic attractors have been found in this model. Applying the new results and techniques in this article to the discrete model is an interesting point for future research work.

The results obtained in the paper show that the complex dynamic behaviors are widespread in some basic discrete models of two host-parasitoid interactions. Some complexities are related to chaotic bands with periodic window, flip bifurcation, Neimark-Sacker bifurcation and mode-locking structure. The others are related to the non-uniqueness of the dynamics, or attractors (meaning that several attractors coexist). As discussed in Section 5, the dynamic behavior of a population may dramatically be affected by small changes in the parameters' values. As a result, in analyzing these biologically realistic models, it is important to choose the suitable values of the parameters. Therefore, identifying complicated, possible chaotic dynamics in population data is still a major challenge in ecological studies.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was partially supported by Tarim University President Fund (TDZKSS201904).

Conflict of interest

All authors declare no conflict of interest.

Appendix A

Appendix B

with $H = a_{12}(u+v), \; P = -(1+a_{11})u+(\lambda_{2}-a_{11})v.$