Bifurcation and chaos in a host-parasitoid model with a lower bound for the host

Abstract

In this paper, a discrete-time biological model and its dynamical behaviors are studied in detail. The existence and stability of the equilibria of the model are qualitatively discussed. More precisely, the conditions for the existence of a flip bifurcation and a Neimark-Sacker bifurcation are derived by using the center manifold theorem and bifurcation theory. Numerical simulations are presented not only to validate our results with the theoretical analysis, but also to exhibit the complex dynamical behaviors. We also analyze the dynamic characteristics of the system in a two-dimensional parameter space. Numerical results indicate that we can more clearly and directly observe the chaotic phenomenon, period-doubling and period-adding, and the optimal parameters matching interval can also be found easily.

1 Introduction

A host-parasitoid model with a lower bound for the host is given by [1] as follows:

$$ \textstyle\begin{cases} H(t+1)=H(t)\exp [\frac{r(1-\frac{H(t)}{k})(H(t)-c)}{H(t)+m}-\frac {abP(t)}{1+aH(t)} ], \\ P(t+1)=H(t) [1-\exp (-\frac{abP(t)}{1+aH(t)} ) ], \end{cases} $$

(1.1)

where $H(t)$ is the host population size in generation t, $P(t)$ is the parasitoid population size in generation t, r is the intrinsic growth rate. k is the carrying capacity of the environment and c is the lower bound for the host. Many species abide by their own ecological genetic rules, and all of them have their own lower bound. Once the population size goes below the lower bound, the species would die out. Therefore, it is necessary to investigate the interspecific interaction with a lower bound. In this model, a is a search rate, b is a conversion factor, and m is a constant. These parameters r, k, c, a, b and m are all positive constants.

The two-dimensional parameter-space has attracted considerable interest in recent years. There have been some studies focusing on detailed treatments of this topic [2–8]. For example, in [2] Sun et al. investigate the emergence of quasi-periodic and mode-locked states of an arbitrary period in pairs of coupled maps of a discrete system. Li et al. [3] study the dynamic behaviors of a discrete system with exponential terms, two control parameters are involved in this system. In their paper, the authors show that there are different kinds of fractal structures in both the attractive and the divergent regions. The two-dimensional parameter space has also been reported by Paulo C Rech [4], who investigates a three quadratic discrete-time system in which transitions are followed by periodic, quasi-periodic, chaotic and hyper-chaotic states. The author uses the maximum Lyapunov exponents, parameter planes, phase space portraits and bifurcation diagrams to verify these transitions. In [8] a new method with a discrete physiological control system is introduced and the dynamical behaviors of the discrete model are investigated. The division diagrams are simulated to catch the dynamic behaviors of the system about two parameters. The above-mentioned models in these literatures are parameter spaces of the discrete system. Similar results are observed in [5–7], they demonstrate the complicated behaviors of the continuous system, which is also controlled by two different parameters numerically.

With the increasing application of ecological models in the real world, more and more attention has been paid to the host-parasitoid models, and they have played an important and fundamental role in the relationship among the biological populations. Although this kind of model has been investigated by many scholars, little research has been carried out on discrete systems [9–13]. In this study, the discrete host-parasitoid system (1.1) is further investigated in detail. We mainly pay our attention to deriving the existence of flip and Neimark-Sacker bifurcations. For simplicity, system (1.1) can be rewritten as follows:

$$ \textstyle\begin{cases} x\rightarrow x \exp [\frac{r(1-\frac{x}{k})(x-c)}{x+m}-\frac {aby}{1+ax} ], \\ y\rightarrow x [1-\exp (-\frac{aby}{1+ax} ) ]. \end{cases} $$

(1.2)

Our objective is to study systematically the existence conditions of the flip bifurcation and the Neimark-Sacker bifurcation, which will be derived using the center manifold theorem and the bifurcation theory (see [14, 15]). The effectiveness of these theoretical analyses is determined by bifurcation diagrams with one control parameter. Furthermore, two-dimensional parameter-space diagrams for system (1.2) are presented to clarify the relationship between complex dynamic behaviors and two control parameters. In fact, the problem considered in this paper and the obtained results can be regarded as a beneficial supplement to the work of [1].

The general outline of this paper is as follows. In Section 2, we discuss the stability criterion of the system at different equilibria in the interior of $R^{2}$. We prove that under certain parametric condition system (1.2) admits a bifurcation in Section 3. Numerical simulations using MATLAB are applied in Section 4 to support the theoretical analyses and visualize the newly observed complex dynamics of the system. These findings prove that there are possibilities for periodic and chaotic motions to exist in the parameter space. In addition, the phase diagrams of two control parameters are also presented. Finally, comments and conclusions are summarized.

2 Stability of equilibria

In this section, we first determine the existence of equilibria of system (1.2) and then investigate their stability by calculating the eigenvalues for the Jacobian matrix at each equilibrium. At last, sufficient conditions for the existence of a flip bifurcation and a Neimark-Sacker bifurcation are derived using a qualitative theorem and bifurcation theory [14, 15].

We can get one equilibrium $E_{0}(0, 0)$ and the other positive equilibrium $E_{1}(x_{1}, y_{1})$, where $x_{1}$, $y_{1}$ satisfy the following equation:

$$ \textstyle\begin{cases} x_{1}=-\frac{q\ln\frac{c}{q}}{a [q\ln\frac{c}{q}+bq-b ]}, \\ y_{1}=\frac{(q-1)\ln\frac{c}{q}}{a [q\ln\frac{c}{q}+bq-b ]}, \end{cases} $$

(2.1)

where $q = \exp [\frac {r(1-x_{1}/k)(x_{1}-c)}{x_{1}+m} ]$.

The Jacobian matrix of system (1.2) at $E_{0}$ is

$$J_0(0,0)=\left(\begin{array}{cc}{e}^{-\frac{rc}{m}}& 0\\ 0& 0\end{array}\right).$$

Accordingly, we can get eigenvalues

$$ \lambda_{1}=e^{-\frac{rc}{m}},\qquad \lambda_{2}=0, $$

from which, we find that $E_{0}$ is a stable node $(|\lambda_{1}|<1)$. The rigorous estimate of the result can be found in [1]. Now we are interested in considering the stability of $E_{1}$. The Jacobian matrix of (1.2) at $E_{1}$ can be written as follows:

$$J_1(x_1,y_1)=\left(\begin{array}{cc}1+r{x}_{1}G+H& -\frac{ab{x}_1}{1+a{x}_1}\\ 1-e^{\frac{ab{y}_1}{1+a{x}_1}}-H{e}^{\frac{ab{y}_1}{1+a{x}_1}}& L\end{array}\right),$$

where $G = \frac{-(x_{1}-c)}{k(x_{1}+m)}+ \frac{1-\frac {x_{1}}{k}}{x_{1}+m}-\frac{(1-\frac {x_{1}}{k})(x_{1}-c)}{(x_{1}+m)^{2}}$, $H=\frac {a^{2}bx_{1}y_{1}}{(1+ax_{1})^{2}}, L=\frac{abx_{1}}{1+ax_{1}}e^{\frac {aby_{1}}{1+ax_{1}}}$.

The characteristic equation of the Jacobian matrix $J_{1}$ is given by

$$ F(\lambda)=\lambda^{2} + P(x_{1}, y_{1})\lambda+ Q(x_{1}, y_{1}) = 0, $$

(2.2)

with coefficients

$$\begin{aligned}& P(x_{1}, y_{1})=- (1+rx_{1}G+H+L ), \\& Q(x_{1}, y_{1}) = rx_{1}LG+\frac{abx_{1}}{1+ax_{1}}. \end{aligned}$$

From (2.2) we have

$$\begin{aligned}& F(1)=\frac{abx_{1}}{1+ax_{1}}-H-L-(x_{1}G-x_{1}LG)r, \\& F(-1)=\frac{abx_{1}}{1+ax_{1}}+H+L+2+(x_{1}G+x_{1}LG)r. \end{aligned}$$

In order to discuss the stability of the equilibrium $E_{1}$, we also need the following lemma, which can be easily proved by the relationship between roots and coefficients of a quadratic equation [9].

Lemma 2.1

Let $F(\lambda) = \lambda^{2} + P\lambda+ Q$. Assume that $F(1) > 0$, $\lambda_{1}$ and $\lambda_{2}$ are the roots of $F(\lambda) = 0$. Then we have the following statements:

• (i) $|\lambda_{1}| < 1$, $|\lambda_{2}| < 1$ if and only if $F(-1) > 0$ and $Q < 1$;

• (ii) $|\lambda_{1}| < 1$, $|\lambda_{2}| > 1$ (or$|\lambda_{1}| > 1$ and $|\lambda_{2}| < 1$) if and only if $F(-1) < 0$;

• (iii) $|\lambda_{1}| > 1$, $|\lambda_{2}| > 1$ if and only if $F(-1) > 0$ and $Q > 1$;

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

• (v) $\lambda_{1}$, $\lambda_{2}$ are complex and $|\lambda_{1}| = |\lambda_{2}| = 1$ if and only if $P^{2} - 4Q < 0$ and $Q = 1$.

Let $\lambda_{1}$ and $\lambda_{2}$ be the roots of (1.2), which are called eigenvalues of the fixed point $(x_{1}, y_{1})$. The fixed point $(x_{1}, y_{1})$ is a sink or locally asymptotically stable if $|\lambda_{1}| < 1$, $|\lambda_{2}| < 1$. The fixed point $(x_{1}, y_{1})$ is a source or locally unstable if $|\lambda_{1}| > 1$, $|\lambda_{2}| > 1$. The fixed point $(x_{1}, y_{1})$ is a saddle if $|\lambda_{1}| < 1$ and $|\lambda_{2}| > 1$ (or $|\lambda_{1}| > 1$ and $|\lambda_{2}| < 1$). The fixed point $(x_{1}, y_{1})$ is non-hyperbolic if either $|\lambda_{1}| = 1$ or $|\lambda_{2}| = 1$.

Then we have the following theorem on the stability of a positive fixed point of system (1.2).

Theorem 2.1

For the positive equilibrium $E_{1}$, we have the following estimates:

• (i) it is a sink if the condition holds: $r > -\frac{2+H+L+\frac {abx_{1}}{1+ax_{1}}}{x_{1}G+x_{1}LG}$ and $r < \frac{1-\frac {abx_{1}}{1+ax_{1}}}{x_{1}LG}$;

• (ii) it is a source if the condition holds: $r > -\frac{2+H+L+\frac {abx_{1}}{1+ax_{1}}}{x_{1}G+x_{1}LG}$ and $r > \frac{1-\frac {abx_{1}}{1+ax_{1}}}{x_{1}LG}$;

• (iii) it is a saddle if the condition holds: $r <-\frac{2+H+L+\frac {abx_{1}}{1+ax_{1}}}{x_{1}G+x_{1}LG}$;

• (iv) it is non-hyperbolic if one of the following conditions holds:

  • (iv.1) $r =- \frac{2+H+L+\frac{abx_{1}}{1+ax_{1}}}{x_{1}G+x_{1}LG}$, and $r \neq-\frac{1+H+L}{x_{1}G}, -\frac{3+H+L}{x_{1}G}$;

  • (iv.2) $r = \frac{1-\frac{abx_{1}}{1+ax_{1}}}{x_{1}GL}$, and $-\frac {3+H+L}{x_{1}G} < r < \frac{1-H-L}{x_{1}G}$.

From the above theorem, if the term (iv.1) holds, one can easily obtain that one of the eigenvalues of $E_{1}$ is −1 and the other is $2 + rx_{1}G + H + L$, which is neither 1 nor −1. If the term (iv.2) holds, then the eigenvalues of $E_{1}$ are a pair of complex conjugate eigenvalues with modulus one. By Theorems 4.3 and 4.6 in [16], we can obtain that when the eigenvalues of the positive point are −1 and a pair of conjugate numbers, the system will undergo a flip bifurcation and N-S bifurcation, respectively.

Let

$$F_{B}= \biggl\{ (r, a, m, k, c, b):r = - \frac{2+H+L+\frac {abx_{1}}{1+ax_{1}}}{x_{1}G+x_{1}LG}, m>0 \biggr\} . $$

The equilibrium $E_{1}(x_{1}, y_{1})$ can arise a flip bifurcation when parameters vary in a small neighborhood of $F_{B}$.

Let

$$H_{B}= \biggl\{ (r, a, m, b, c, k):r = \frac{1-\frac {abx_{1}}{1+ax_{1}}}{x_{1}GL}, \frac{-(3+H+L)}{x_{1}G}< r < \frac {1-H-L}{x_{1}G}, m>0 \biggr\} . $$

The equilibrium $E_{1}(x_{1}, y_{1})$ can lead to the bifurcation of Neimark-Sacker when parameters are restricted to a small neighborhood of $H_{B}$. In the following section, we will investigate the flip bifurcation of $E_{1}$ if parameters vary in a small vicinity of $F_{B}$ and the Neimark-Sacker bifurcation of $E_{1}$ if parameters lie in a small scope of $H_{B}$.

3 Bifurcation

Based on the above-mentioned analysis, in this section, we mainly focus on the flip bifurcation and Neimark-Sacker bifurcation of $E_{1}$ because they are correlated with the local and global stability of system (1.2). We choose r as a bifurcation parameter for studying the bifurcations.

3.1 Flip bifurcation

We first discuss the flip bifurcation of system (1.2) at $E_{1}(x_{1},y_{1})$ when parameters vary in a small neighborhood of $F_{B}$ [10, 17, 18]. Select arbitrary parameters $(a, r_{1}, m, k, b, c)$ from $F_{B}$. System (1.2) has a unique positive equilibrium $E_{1}(x_{1}, y_{1})$, the corresponding eigenvalues are $\lambda _{1}=-1, \lambda_{2}= 2 + rx_{1}G + H + L$ with $|\lambda_{2}|\neq1$.

Let $u=x-x_{1}$, $v=y-y_{1}$, $r^{*}=r-r_{1}$, the equilibrium $(x_{1},y_{1})$ is transformed to the origin point $(0,0)$. We consider the parameter $r^{*}$ as a new dependent variable, and we can get

$$\left(\begin{array}{c}u\\ r^{\ast }\\ v\end{array}\right)\to \left(\begin{array}{c}{a}_{11}u+a_{12}{r}^{\ast }+a_{13}v+a_0{u}^2+a_{1}uv+a_{2}u{r}^{\ast }+a_{3}v{r}^{\ast }+a_4{u}^3+a_5{u}^{2}v+a_6{u}^2{r}^{\ast }\\ +a_{7}uv{r}^{\ast }+a_8{v}^2+a_{9}u{v}^2+a_{10}{v}^3+O(4)\\ -r^{\ast }\\ a_{31}u+a_{33}v+b_0{u}^2+b_{1}uv+b_4{u}^3+b_5{u}^{2}v+b_8{v}^2+b_{9}u{v}^2+b_{10}{v}^3+O(4)\end{array}\right),$$

(3.1)

where

$$\begin{aligned}& a_{11} = 1+rx_{1}G+H,\qquad a_{12}=\frac{x_{1}(1-\frac {x_{1}}{k})(x_{1}-c)}{m+x_{1}},\qquad a_{13}=-\frac {abx_{1}}{1+ax_{1}}, \\& a_{0} = r \biggl(G+GHx_{1}-\frac{x_{1}(kG+1)}{k(m+x_{1})} \biggr)+ \frac {x_{1}H^{2}}{2}+H+\frac{aH}{1+ax_{1}}, \\& a_{1} = -\frac{abx_{1}G}{1+ax_{1}}r-\frac{abx_{1}H}{1+x_{1}}-\frac {ab}{(1+ax_{1})^{2}}, \\& a_{2} = Gx_{1}, \qquad a_{3}=0,\qquad a_{31}=1-He^{-\frac {aby_{1}}{1+ax_{1}}}-e^{-\frac{aby_{1}}{1+ax_{1}}},\qquad a_{33}=L, \\& a_{4} = r \biggl(\frac{a^{2}by_{1}G}{(1+ax_{1})^{3}}+\frac {a^{4}b^{2}x_{1}y_{1}^{2}G}{2(1+ax_{1})^{4}}- \frac {m(1+kG)}{k(m+x_{1})^{2}}-\frac{H(kG+1)}{k(m+x_{1})} \biggr)+\frac {a^{4}b^{2}y_{1}^{2}G}{2(1+ax_{1})^{4}} \\& \hphantom{a_{4} ={}}{}-\frac{a^{3}by_{1}}{(1+ax_{1})^{4}}-\frac {a^{5}b^{3}x_{1}y_{1}^{2}}{(1+ax_{1})^{5}}+\frac {a^{6}b^{3}x_{1}y_{1}^{3}}{6(1+ax_{1})^{6}}, \\& a_{5} = r \biggl(\frac{abx_{1}(1+kG)}{k(1+ax_{1})(m+x_{1})}-\frac {abG}{(1+ax_{1})^{2}}- \frac {a^{3}b^{2}x_{1}y_{1}G}{2(1+ax_{1})^{3}} \biggr)+\frac {a^{2}b}{(1+ax_{1})^{3}} \\& \hphantom{a_{5} ={}}{}-\frac{a^{3}b^{2}y_{1}}{(1+ax_{1})^{4}}-\frac {a^{4}b^{2}x_{1}y_{1}}{(1+ax_{1})^{4}}-\frac {a^{5}b^{3}x_{1}y_{1}^{2}}{2(1+ax_{1})^{5}}, \\& a_{6} = G+GHx_{1}-\frac{x_{1}(kG+1)}{k(m+x_{1})}, \qquad a_{7}=\frac {-abx_{1}G}{1+ax_{1}}, \qquad a_{8}=\frac {a^{2}b^{2}x_{1}}{2(1+ax_{1})^{2}}, \\& a_{9} = \frac{a^{2}b^{2}x_{1}G}{2(1+ax_{1})^{2}}r+\frac {a^{2}b^{2}(1+x_{1}H)}{2(1+ax_{1})^{2}}-\frac {a^{3}b^{2}x_{1}}{(1+ax_{1})^{3}}, \qquad a_{10}=\frac {-a^{3}b^{3}x_{1}}{6(1+ax_{1})^{3}}, \\& b_{0} = e^{-\frac{aby_{1}}{1+ax_{1}}} \biggl(\frac {-a^{2}by_{1}}{(1+ax_{1})^{3}}-\frac {a^{4}b^{2}x_{1}y_{1}^{2}}{2(1+ax_{1})^{4}} \biggr),\qquad b_{1}=e^{-\frac {aby_{1}}{1+ax_{1}}} \biggl(\frac {a^{3}b^{2}x_{1}y_{1}}{{(1+ax_{1})^{3}}}+ \frac{ab}{(1+ax_{1})^{2}} \biggr), \\& b_{4} = e^{-\frac{aby_{1}}{1+ax_{1}}} \biggl(\frac {-a^{6}b^{3}x_{1}y_{1}^{3}}{6(1+ax_{1})^{6}}+\frac {a^{5}b^{2}x_{1}y_{1}^{2}}{(1+ax_{1})^{5}}+ \frac {a^{3}by_{1}}{(1+ax_{1})^{4}}-\frac {a^{4}b^{2}y_{1}^{2}}{2(1+ax_{1})^{4}} \biggr), \\& b_{5} = e^{-\frac{aby_{1}}{1+ax_{1}}} \biggl(\frac {a^{3}b^{2}y_{1}}{(1+ax_{1})^{3}}-\frac{a^{2}b}{(1+ax_{1})^{3}}+ \frac {a^{5}b^{3}x_{1}y_{1}^{2}}{2(1+ax_{1})^{5}}-\frac {2a^{4}b^{2}x_{1}y_{1}}{(1+ax_{1})^{4}} \biggr), \\& b_{8} = -e^{-\frac{aby_{1}}{1+ax_{1}}}\frac {a^{2}b^{2}x_{1}}{2(1+ax_{1})^{2}},\qquad b_{10}= \frac {a^{3}b^{3}x_{1}}{6(1+ax_{1})^{3}}e^{-\frac {aby_{1}}{1+ax_{1}}}, \\& b_{9} = e^{-\frac{aby_{1}}{1+ax_{1}}} \biggl(\frac {a^{3}b^{2}x_{1}}{(1+ax_{1})^{3}}-\frac {a^{2}b^{2}}{2(1+ax_{1})^{2}}- \frac {a^{4}b^{3}x_{1}y_{1}}{2(1+ax_{1})^{4}} \biggr). \end{aligned}$$

(3.2)

If $a_{11}+a_{33}+2\neq0$, we make an invertible matrix T as follows:

$$T=\left(\begin{array}{ccc}-1-a_{33}& 1& \frac{{\lambda }_2-a_{33}}{a_{31}}\\ 0& -\frac{a_{11}+a_{33}+2}{a_{12}}& 0\\ a_{31}& 0& 1\end{array}\right).$$

Using translation

$$\left(\begin{array}{c}u\\ r^{\ast }\\ v\end{array}\right)=T\left(\begin{array}{c}\overline{x}\\ \overline{r^{\ast }}\\ \overline{y}\end{array}\right),$$

(3.3)

then system (3.1) reduces to

$$\left(\begin{array}{c}\overline{x}\\ \overline{r^{\ast }}\\ \overline{y}\end{array}\right)\to \left(\begin{array}{ccc}-1& 1& 0\\ 0& -1& 0\\ 0& 0& {\lambda }_2\end{array}\right)\left(\begin{array}{c}\overline{x}\\ \overline{r^{\ast }}\\ \overline{y}\end{array}\right)+\left(\begin{array}{c}f(\overline{x},\overline{y},\overline{r^{\ast }})\\ 0\\ g(\overline{x},\overline{y},\overline{r^{\ast }})\end{array}\right),$$

(3.4)

where

$$\begin{aligned}& \begin{aligned} f\bigl(\bar{x}, \bar{y},\bar{r^{*}}\bigr) =&{} \biggl( \frac{b_{0}(\lambda _{2}-a_{33})}{a_{31}(1+\lambda_{2})}-\frac{a_{0}}{1+\lambda_{2}} \biggr)u^{2}+ \biggl( \frac{b_{1}(\lambda_{2}-a_{33})}{a_{31}(1+\lambda _{2})}-\frac{a_{1}}{1+\lambda_{2}} \biggr)uv -\frac{a_{2}}{1+\lambda_{2}}ur^{*} \\ &{} + \biggl(\frac{b_{4}(\lambda_{2}-a_{33})}{a_{31}(1+\lambda_{2})}-\frac {a_{4}}{1+\lambda_{2}} \biggr)u^{3}+ \biggl(\frac{b_{5}(\lambda _{2}-a_{33})}{a_{31}(1+\lambda_{2})}-\frac{a_{5}}{1+\lambda_{2}} \biggr)u^{2}v - \frac{a_{6}}{1+\lambda_{2}}u^{2}r^{*} \\ &{} + \biggl(\frac{b_{8}(\lambda_{2}-a_{33})}{a_{31}(1+\lambda_{2})}-\frac {a_{8}}{1+\lambda_{2}} \biggr)v^{2} + \biggl(\frac{b_{9}(\lambda_{2}-a_{33})}{a_{31}(1+\lambda_{2})}-\frac {a_{9}}{1+\lambda_{2}} \biggr)uv^{2}- \frac{a_{7}}{1+\lambda_{2}}uvr^{*} \\ &{} + \biggl(\frac{b_{10}(\lambda_{2}-a_{33})}{a_{31}(1+\lambda_{2})}-\frac {a_{10}}{1+\lambda_{2}} \biggr)v^{3}+O \bigl(\bigl( \vert u \vert + \vert v \vert + \bigl\vert r^{*} \bigr\vert \bigr)^{4}\bigr), \end{aligned} \\& \begin{aligned} g\bigl(\bar{x}, \bar{y},\bar{r^{*}}\bigr) = {}& \biggl(\frac{a_{31}a_{0}}{1+\lambda _{2}}+\frac{b_{0}(1+a_{33})}{1+\lambda_{2}} \biggr)u^{2}+ \biggl( \frac {a_{31}a_{1}}{1+\lambda_{2}}+\frac{b_{1}(1+a_{33})}{1+\lambda_{2}} \biggr)uv +\frac{a_{31}a_{2}}{1+\lambda_{2}}ur^{*} \\ &{} + \biggl(\frac{a_{31}a_{4}}{1+\lambda_{2}}+\frac {b_{4}(1+a_{33})}{1+\lambda_{2}} \biggr)u^{3}+ \biggl(\frac {a_{31}a_{5}}{1+\lambda_{2}}+\frac{b_{5}(1+a_{33})}{1+\lambda_{2}} \biggr)u^{2}v + \frac{a_{31}a_{6}}{1+\lambda_{2}}u^{2}r^{*} \\ &{} + \biggl(\frac{a_{31}a_{8}}{1+\lambda_{2}}+\frac {b_{8}(1+a_{33})}{1+\lambda_{2}} \biggr)v^{2}+ \biggl(\frac {a_{31}a_{9}}{1+\lambda_{2}}+\frac{b_{9}(1+a_{33})}{1+\lambda_{2}} \biggr)uv^{2} + \frac{a_{31}a_{7}}{1+\lambda_{2}}uvr^{*} \\ &{} + \biggl(\frac{a_{31}a_{10}}{1+\lambda_{2}}+\frac {b_{10}(1+a_{33})}{1+\lambda_{2}} \biggr)v^{3}+O \bigl(\bigl( \vert u \vert + \vert v \vert + \bigl\vert r^{*} \bigr\vert \bigr)^{4}\bigr), \end{aligned} \\& u = (-1-a_{33})\bar{x}+\bar{r^{*}}+\frac{\lambda_{2}-a_{33}}{a_{31}}\bar {y},\qquad r^{*}=-\frac{a_{11}+a_{33}+2}{a_{12}}\bar{r^{*}},\qquad v=a_{31}\bar {x}+\bar{y}. \end{aligned}$$

Next, we search for the center manifold of (3.4) at the point $(0,0)$ in a small neighborhood of $r^{*}=0$. By a center manifold theorem (Theorem 7 in [14]), we know that there exists a center manifold $W^{c}(0,0,0)$, which can be approximately represented as follows:

$$ W^{c}(0,0,0)=\bigl\{ \bigl(\bar{x}, \bar{y}, \bar{r^{*}} \bigr)\in R^{3}: \bar {y}=c_{1}\bar{x}^{2} + c_{2}\bar{x}\bar{r^{*}} + c_{3} \bar{r^{*2}} + O\bigl(\bigl( \vert \bar{x} \vert + \bigl\vert \bar{r^{*}} \bigr\vert \bigr)^{3}\bigr)\bigr\} , $$

where $O((|\bar{x}|+|\bar{r^{*}}|)^{3})$ is the sum of all terms whose order is greater than 2.

After calculations, the following results can be easily obtained:

$$\begin{aligned}& c_{1} = \frac {(1+a_{33})^{2}[a_{31}a_{0}+b_{0}(1+a_{33})]-a_{31}(1+a_{33})[a_{31}a_{1}+b_{1}(1+a_{33})]+a_{31}^{2}[a_{31}a_{8}+b_{8}(1+a_{33})]}{(1-\lambda _{2}^{2})}, \\& c_{2} = \frac {a_{31}[a_{31}a_{1}+b_{1}(1+a_{33})]-2(1+a_{33})[a_{31}a_{0}+b_{0}(1+a_{33})]}{(1-\lambda _{2})^{2}} \\& \hphantom{c_{2} ={}}{}+\frac {a_{31}a_{2}(1+a_{33})(a_{11}+a_{33}+2)}{a_{12}(1-\lambda_{2})^{2}}, \\& c_{3} = \frac{a_{31}a_{0}+b_{0}(1+a_{33})}{(1-\lambda_{2})^{2}}-\frac {a_{31}a_{2}(a_{11}+a_{33}+2)}{a_{12}(1-\lambda_{2})^{2}}. \end{aligned}$$

Therefore, we consider the following map originating from (3.4) restricted to the center manifold $W^{c}(0,0,0)$.

$$\begin{aligned}& F: \bar{x} \rightarrow-\bar{x} + \bar{r^{*}}+h_{1} \bar{x}^{2} + h_{2}\bar{x}\bar{r^{*}} + h_{3}\bar{r^{*}}^{2} + h_{4} \bar{x}^{3} + h_{5}\bar{x}^{2}\bar{r^{*}} \\& \quad {}+ h_{6}\bar{x}\bar{r^{*}}^{2}+h_{7}\bar {r^{*}}^{3}+O \bigl(\bigl( \vert \bar{x} \vert + \bigl\vert \bar{r^{*}} \bigr\vert \bigr)^{4} \bigr) , \end{aligned}$$

(3.5)

where

$$\begin{aligned}& h_{1} = \frac{1}{1+\lambda_{2}}\bigl\{ (1+a_{33})^{2} \bigl[a_{31}a_{0}+b_{0}(1+a_{33}) \bigr]-a_{31}(1+a_{33}) \bigl[a_{31}a_{1}+b_{1}(1+a_{33}) \bigr] \\& \hphantom{h_{1} ={}}{} + a_{31}^{2}\bigl[a_{31}a_{8}+b_{8}(1+a_{33}) \bigr]\bigr\} , \\& h_{2} = \frac{1}{1+\lambda_{2}}\bigl\{ a_{31} \bigl[a_{31}a_{1}+b_{1}(1+a_{33}) \bigr]-2(1+a_{33})\bigl[a_{31}a_{0}+b_{0}(1+a_{33}) \bigr] \\& \hphantom{h_{2} ={}}{} + a_{31}a_{2}(1+a_{33}) (a_{11}+a_{33}+2)\bigr\} , \\& h_{3} = \frac{1}{1+\lambda_{2}}\bigl\{ a_{31}a_{0}+b_{0}(1+a_{33})-a_{31}a_{2}(a_{11}+a_{33}+2) \bigr\} , \\& h_{4} = \frac{1}{1+\lambda_{2}}\bigl\{ \bigl[a_{31}a_{1}+b_{1}(1+a_{33}) \bigr](\lambda _{2}-1-2a_{33})c_{1} \\& \hphantom{h_{4} ={}}{}-2(1+a_{33}) (\lambda _{2}-a_{33})\bigl[a_{31}a_{0}+b_{0}(1+a_{33}) \bigr]c_{1} \\& \hphantom{h_{4} ={}}{} + 2a_{31}\bigl[a_{31}a_{8}+b_{8}(1+a_{33}) \bigr]c_{1}-(1+a_{33})^{3}\bigl[a_{31}a_{4}+b_{4}(1+a_{33}) \bigr] \\& \hphantom{h_{4} ={}}{}+a_{31}^{3}\bigl[a_{31}a_{10}+b_{10}(1+a_{33}) \bigr] \\& \hphantom{h_{4} ={}}{} + a_{31}^{2}(1+a_{33}) \bigl[a_{31}a_{9}+b_{9}(1+a_{33}) \bigr]+a_{31}(1+a_{33})^{2}\bigl[a_{31}a_{5}+b_{5}(1+a_{33}) \bigr]\bigr\} , \\& h_{5} = \frac{1}{1+\lambda_{2}}\biggl\{ \bigl[a_{31}a_{1}+b_{1}(1+a_{33}) \bigr](\lambda _{2}-1-2a_{33})c_{2} \\& \hphantom{h_{5} ={}}{}-2(1+a_{33}) (\lambda _{2}-a_{33})\bigl[a_{31}a_{0}+b_{0}(1+a_{33}) \bigr]c_{2} \\& \hphantom{h_{5} ={}}{} + 2a_{31}\bigl[a_{31}a_{8}+b_{8}(1+a_{33}) \bigr]c_{2}+\bigl[a_{31}a_{1}+b_{1}(1+a_{33}) \bigr]c_{1} \\& \hphantom{h_{5} ={}}{}-\frac {a_{2}a_{31}(\lambda_{2}-a_{33})(a_{11}+a_{33}+2)}{a_{12}a_{31}}c_{1} \\& \hphantom{h_{5} ={}}{} + 3(1+a_{33})^{2}\bigl[a_{31}a_{4}+b_{4}(1+a_{33}) \bigr]+a_{31}\bigl[a_{31}a_{9}+b_{9}(1+a_{33}) \bigr] \\& \hphantom{h_{5} ={}}{}-2(1+a_{33})\bigl[a_{31}a_{5}+b_{5}(1+a_{33}) \bigr] \\& \hphantom{h_{5} ={}}{} + \frac{a_{7}a_{31}^{2}(1+a_{33})(a_{11}+a_{33}+2)}{a_{12}}-\frac {a_{6}a_{31}(a_{11}+a_{33}+2)(1+a_{33})^{2}}{a_{12}}\biggr\} , \\& h_{6} = \frac{1}{1+\lambda_{2}}\biggl\{ \bigl[a_{31}a_{1}+b_{1}(1+a_{33}) \bigr](\lambda _{2}-1-2a_{33})c_{3} \\& \hphantom{h_{6} ={}}{}-2(1+a_{33}) (\lambda _{2}-a_{33})\bigl[a_{31}a_{0}+b_{0}(1+a_{33}) \bigr]c_{3} \\& \hphantom{h_{6} ={}}{} + 2a_{31}\bigl[a_{31}a_{8}+b_{8}(1+a_{33}) \bigr]c_{3}+\bigl[a_{31}a_{1}+b_{1}(1+a_{33}) \bigr]c_{2} \\& \hphantom{h_{6} ={}}{}-\frac {a_{2}a_{31}(\lambda_{2}-a_{33})(a_{11}+a_{33}+2)}{a_{12}a_{31}}c_{2}- 3(1+a_{33})\bigl[a_{31}a_{4}+b_{4}(1+a_{33}) \bigr] \\& \hphantom{h_{6} ={}}{}+a_{31}\bigl[a_{31}a_{5}+b_{5}(1+a_{33}) \bigr]-\frac {a_{7}a_{31}^{2}(a_{11}+a_{33}+2)}{a_{12}} \\& \hphantom{h_{6} ={}}{} + \frac{2a_{31}a_{6}(1+a_{33}(a_{11}+a_{33}+2))}{a_{12}}\biggr\} , \\& h_{7} = \frac{1}{1+\lambda_{2}}\biggl\{ \bigl[a_{31}a_{1}+b_{1}(1+a_{33}) \bigr]c_{3}-\frac{a_{31}a_{2}(\lambda _{2}-a_{33})(a_{11}+a_{33}+2)}{a_{12}a_{31}}c_{3} \\& \hphantom{h_{7} ={}}{}+\bigl[a_{31}a_{4}+b_{4}(1+a_{33}) \bigr] - \frac{a_{31}a_{6}(a_{11}+a_{33}+2)}{a_{12}}\biggr\} . \end{aligned}$$

In order to ensure that map (3.5) appears as a flip bifurcation, we require two discriminatory quantities $\alpha_{1}$ and $\alpha_{2}$ to be not zero, where

$$ \alpha_{1}= \biggl(\frac{2\partial^{2}F}{\partial\bar{x}\, \partial r^{\ast }}+\frac{\partial F}{\partial r^{\ast}} \frac{\partial^{2}F}{\partial\bar {x}^{2}} \biggr)_{(0, 0)},\qquad \alpha_{2}= \biggl(\frac{1}{3}\frac{\partial^{3}F}{\partial\bar {x}^{3}}+ \frac{1}{2} \biggl(\frac{\partial^{2} F}{\partial\bar {x}^{2}} \biggr)^{2} \biggr)_{(0, 0)}. $$

(3.6)

According to the above discussions and applying the bifurcation theory presented as Section 3.2 in Guckenheimer [15], we state the following result on a flip bifurcation.

Theorem 3.1

If $\alpha_{2}\neq0$, the parameter $r^{\ast}$ alters in the limited region of the point $(0,0)$, then system (1.2) undergoes a flip bifurcation at $E_{1}(x_{1}, y_{1})$. Moreover, the period-2 orbit that bifurcates from $E_{1}(x_{1}, y_{1})$ is stable (unstable) if $\alpha_{2}> 0$ ($\alpha_{2}< 0$).

3.2 N-S bifurcation

Next, we will give attention to recapitulating the conditions of the existence for a Neimark-Sacker bifurcation (discrete Hopf bifurcation) by using the bifurcation theorem [11, 12, 19].

Considering system (1.2) with arbitrary parameters $(a, r_{2}, m, k, b, c)\in H_{B}$, we write system (1.2) in the form

$$ \textstyle\begin{cases} x\rightarrow x\exp [\frac{r_{2}(1-\frac{x}{k})(x-c)}{x+m}-\frac {aby}{1+ax} ], \\ y\rightarrow x [1-\exp (-\frac{aby}{1+ax} ) ]. \end{cases} $$

(3.7)

$E_{1}(x_{1}, y_{1})$ is the only positive equilibrium of system (3.7) given by (2.1) and

$$r_{2} = \frac{1}{x_{1}GL} \biggl(1-\frac{abx_{1}}{1+ax_{1}} \biggr). $$

We consider the perturbation of (3.7), as shown below, with $\bar {r}^{*}$ as a bifurcation parameter.

$$ \textstyle\begin{cases} x\rightarrow x\exp [\frac{(\bar{r^{*}}+r_{2})(1-\frac {x}{k})(x-c)}{x+m}-\frac{aby}{1+ax} ], \\ y\rightarrow x [1-\exp (-\frac{aby}{1+ax} ) ], \end{cases} $$

(3.8)

where $|\bar{r^{*}}|\ll 1$, which is a limited perturbation parameter.

Let $u=x-x_{1}$, $v=y-y_{1}$. After shifting the fixed point $(x_{1}, y_{1})$ to the origin, we can get

$$\left(\begin{array}{c}u\\ v\end{array}\right)\to \left(\begin{array}{c}{a}_{11}u+a_{13}v+a_0{u}^2+a_{1}uv+a_4{u}^3+a_5{u}^{2}v+a_8{v}^2+a_{9}u{v}^2+a_{10}{v}^3+O(4)\\ a_{31}u+a_{33}v+b_0{u}^2+b_{1}uv+b_4{u}^3+b_5{u}^{2}v+b_8{v}^2+b_{9}u{v}^2+b_{10}{v}^3+O(4)\end{array}\right),$$

(3.9)

where $a_{11}$, $a_{13}$, $a_{0}$, $a_{1}$, $a_{4}$, $a_{5}$, $a_{8}$, $a_{9}$, $a_{10}$, $a_{31}$, $a_{33}$, $b_{0}$, $b_{1}$, $b_{4}$, $b_{5}$, $b_{8}$, $b_{9}$, $b_{10}$ are given in (3.2), and $r=r_{2}+\bar{r}^{*}$. The characteristic equation associated with the linearization of system (3.9) at $(0,0)$ is given by

$$\lambda^{2} + p\bigl(\bar{r}^{*}\bigr)\lambda+ q\bigl( \bar{r}^{*}\bigr)=0, $$

with coefficients

$$\begin{aligned}& p\bigl(\overline{r}^{*}\bigr) = - \bigl[1+\bigl(r_{2}+ \bar{r^{*}}\bigr)Gx_{1}+H+L \bigr], \\& q\bigl(\overline{r}^{*}\bigr) = \bigl(r_{2}+ \bar{r^{*}}\bigr)x_{1}GL+\frac{abx_{1}}{1+ax_{1}}. \end{aligned}$$

Since the parameters belong to $H_{B}$, the roots of the characteristic equation are

$$ \lambda, \bar{\lambda}= -\frac{p(\bar{r}^{*})}{2} \pm \frac{i}{2}\sqrt { 4q\bigl(\bar{r}^{*}\bigr)-p^{2} \bigl(\bar{r}^{*}\bigr)}. $$

(3.10)

Then we have

$$|\lambda|=\sqrt{q\bigl(\bar{r}^{*}\bigr)},\qquad l=\frac{d|\lambda|}{d\bar{r}^{*}} \bigg| _{\bar{r}^{*}=0} = \frac{x_{1}LG}{2} \neq0. $$

In addition, it is required that if $\bar{r}^{*}=0$, we have $\lambda^{m}, \bar{\lambda}^{m} \neq1$ ($m=1,2,3,4$), i.e., nondegeneracy condition, which is equivalent to $p(0)\neq-2,0,1,2$. Note that $(a, m, k, r_{2}, b, c) \in H _{B}$, we have $p(0)\neq-2,2$. We only need to require $p(0)\neq0,1$ which leads to

$$ -L(3+H+L), -L(1+H+L) \neq1-\frac{abx_{1}}{1+ax_{1}}. $$

(3.11)

Therefore, the eigenvalues λ, λ̄ do not lie in the intersection of the unit circle with the coordinate axes when $\bar {r}^{*} =0$ and condition (3.11) holds.

Letting $\bar{r}^{*} =0$, $\mu= -\frac{p(0)}{2}$, $\omega=\frac{\sqrt{ 4q(0)-p^{2}(0)}}{2}$, we construct an invertible matrix

$$T=\left(\begin{array}{cc}{a}_{12}& 0\\ \mu -a_{11}& -\omega \end{array}\right),$$

(3.12)

and use the following translation:

$$\left(\begin{array}{c}u\\ v\end{array}\right)=T\left(\begin{array}{c}\overline{x}\\ \overline{y}\end{array}\right).$$

(3.13)

Then system (3.9) becomes of the following form:

$$\left(\begin{array}{c}\overline{x}\\ \overline{y}\end{array}\right)\to \left(\begin{array}{cc}\mu & -\omega \\ \omega & \mu \end{array}\right)\left(\begin{array}{c}\overline{x}\\ \overline{y}\end{array}\right)+\left(\begin{array}{c}\overline{f}(\overline{x},\overline{y})\\ \overline{g}(\overline{x},\overline{y})\end{array}\right),$$

(3.14)

where

$$\begin{aligned}& \bar{f}(\bar{x},\bar{y}) = \frac{1}{a_{12}} \bigl(a_{0}u^{2}+a_{1}uv+a_{4}u^{3}+a_{5}u^{2}v+a_{8}v^{2}+a_{9}uv^{2}+a_{10}v^{3} \bigr)+O \bigl(\bigl( \vert u \vert + \vert v \vert \bigr)^{4} \bigr), \\& \bar{g}(\bar{x},\bar{y}) = \frac{1}{a_{12}\omega}\bigl\{ \bigl[a_{0}( \mu-a _{11})-a_{12}b_{0} \bigr]u^{2}+ \bigl[a_{1}(\mu-a _{11})-a_{12}b_{1} \bigr]uv \\& \hphantom{\bar{g}(\bar{x},\bar{y})={}}{}+ \bigl[a_{8}(\mu -a _{11})-a_{12}b_{8} \bigr]v^{2} + \bigl[a_{4}(\mu-a _{11})-a_{12}b_{4} \bigr]u^{3} \\& \hphantom{\bar{g}(\bar{x},\bar{y})={}}{}+ \bigl[a_{5}(\mu-a _{11})-a_{12}b_{5} \bigr]u^{2}v+ \bigl[a_{9}(\mu-a _{11})-a_{12}b_{9} \bigr]uv^{2} \\& \hphantom{\bar{g}(\bar{x},\bar{y})={}}{}+ \bigl[a_{10}(\mu-a _{11})-a_{12}b_{10} \bigr]v^{3}\bigr\} + O \bigl(\bigl( \vert u \vert + \vert v \vert \bigr)^{4} \bigr), \end{aligned}$$

and

$$\begin{aligned}& u^{2}=a_{12}^{2}\bar{x}^{2}, \\& uv=a_{12} (\mu-a _{11} )\bar{x}^{2} - a_{12}\omega\bar {x}\bar{y}, \\& v^{2} = (\mu-a _{11} )^{2}\bar{x}^{2} - 2\omega (\mu -a _{11} )\bar{x}\bar{y} + \omega^{2} \bar{y}^{2}, \\& u^{3} = a_{12}^{3}\bar{x}^{3}, \\& u^{2}v = a_{12}^{2} (\mu-a _{11} ) \bar{x}^{3} - a_{12}^{2}\omega\bar{x}^{2} \bar{y}, \\& uv^{2} = a_{12} (\mu-a _{11} )^{2} \bar{x}^{3} - 2a_{12}\omega\bar{x}^{2}\bar{y} + a_{12}\omega^{2}\bar{x}\bar{y}^{2}, \\& v^{3} = (\mu-a _{11} )\bar{x}^{3} - \omega^{3}\bar{y}^{3} - 3\omega (\mu-a _{11} )^{2}\bar{x}^{2}\bar{y} + 3\omega ^{2} (\mu-a _{11} )\bar{x}\bar{y}^{2}. \end{aligned}$$

Therefore

$$\begin{aligned}& \bar{f}_{\bar{x}\bar{x}} = 2a_{12}\bigl[a_{0}+a_{1}( \mu-a _{11})\bigr]+\frac {2a_{8}(\mu-a _{11})^{2}}{a_{12}}, \\& \bar{f}_{\bar{x}\bar{y}} = -\omega a_{1}-\frac{2a_{8}\omega(\mu-a _{11})}{a_{12}}, \qquad \bar{f}_{\bar{y}\bar{y}} =\frac{2a_{8}\omega ^{2}}{a_{12}}, \\& \bar{f}_{\bar{x}\bar{x}\bar{x}} =6a_{12}\bigl[a_{5}(\mu-a _{11})+a_{4}a_{12}\bigr]+\frac{6}{a_{12}} \bigl[a_{10}(\mu-a _{11})+a_{12}a_{9}( \mu-a _{11})^{2}\bigr], \\& \bar{f}_{\bar{x}\bar{x}\bar{y}} = -2\omega a_{12}a_{5}-4a_{9} \omega (\mu-a _{11})-\frac{6a_{10}\omega}{a_{12}}(\mu-a _{11})^{2}, \\& \bar{f}_{\bar{y}\bar{y}\bar{y}}= \frac{-6a_{10}\omega^{3}}{a_{12}},\qquad \bar{f}_{\bar{x}\bar{y}\bar{y}} = 2 \omega^{2}a_{9}+\frac{6a_{1}\omega ^{2}}{a_{12}}(\mu-a _{11}), \\& \bar{g}_{\bar{x}\bar{x}}=\frac{2}{\omega} \bigl\{ a_{12} \bigl[a_{0}(\mu-a _{11})-a_{12}b_{0} \bigr]+(\mu-a _{11})\bigl[a_{1}(\mu-a _{11})-a_{12}b_{1} \bigr] \\& \hphantom{\bar{g}_{\bar{x}\bar{x}}={}}{}-(\mu-a _{11})^{2}\bigl[a_{8}(\mu-a _{11})-a_{12}b_{8}\bigr] \bigr\} , \\& \bar{g}_{\bar{y}\bar{y}}=\frac{2\omega}{a_{12}}\bigl[a_{8}(\mu-a _{11})-a_{12}b_{8}\bigr], \\& \bar{g}_{\bar{x}\bar{y}}=a_{12}b_{1}-(\mu-a _{11}) \bigl[2a_{8}(\mu-a _{11})+2a_{12}b_{8}-a _{1}\bigr], \\& \bar{g}_{\bar{x}\bar{x}\bar{x}} = \frac{6a_{12}}{\omega}\biggl\{ a_{12} \bigl[a_{4}(\mu-a _{11})-a_{12}b_{4} \bigr]+(\mu-a _{11})\bigl[a_{5}(\mu-a _{11})-a_{12b_{5}} \bigr] \\& \hphantom{\bar{g}_{\bar{x}\bar{x}\bar{x}} ={}}{}+\frac{6}{a_{12}\omega}\bigl[a_{9}(\mu-a_{11})^{3}-a_{12}b_{9}(\mu-a _{11})^{2}+a_{10}( \mu-a _{11})^{2}-a_{12}b_{10}(\mu-a _{11})\bigr]\biggr\} , \\& \bar{g}_{\bar{x}\bar{x}\bar{y}}=2a_{12}\bigl[a_{12}b_{5}-a_{5}( \mu-a _{11})\bigr]-\frac{2}{a_{12}}\bigl\{ 2a_{12}(\mu-a _{11})\bigl[a_{9}(\mu-a _{11})-a_{12}b_{9} \bigr] \\& \hphantom{\bar{g}_{\bar{x}\bar{x}\bar{y}}={}}{}+3(\mu-a _{11})^{2}\bigl[a_{10}( \mu-a _{11})-a_{12}b_{10}\bigr]\bigr\} , \\& \bar{g}_{\bar{x}\bar{y}\bar{y}}=2\omega\bigl\{ \bigl[a_{9}(\mu-a _{11})-a_{12}b_{9}\bigr]a_{12}+3(\mu-a _{11})\bigl[a_{10}(\mu -a_{11})-a_{12}b_{10} \bigr]\bigr\} , \\& \bar{g}_{\bar{y}\bar{y}\bar{y}}= \frac{6\omega^{2}}{a_{12}}\bigl[a_{10}(\mu -a_{11})-a_{12}b_{10}\bigr]. \end{aligned}$$

In order to ensure that map (3.9) undergoes a Neimark-Sacker bifurcation, we require the following discriminatory quantity θ to be not zero.

$$ \theta= - \biggl[\operatorname{Re} \biggl(\frac{(1-2\lambda)\bar{\lambda}^{2}}{1-\lambda }\xi_{20} \xi_{11} \biggr)-\frac{1}{2}|\xi_{11}|^{2}-| \xi _{02}|^{2}+\operatorname{Re}(\bar{\lambda} \xi_{21}) \biggr]_{\bar{r}^{*}=0}, $$

where

$$\begin{aligned}& \xi_{20} = \frac{1}{8} \bigl[\bar{f}_{\bar{x}\bar{x}}- \bar{f}_{\bar {y}\bar{y}}+2\bar{g}_{\bar{x}\bar{y}}+i (\bar{g}_{\bar{x}\bar {x}}- \bar{g}_{\bar{y}\bar{y}}-2\bar{f}_{\bar{x}\bar{y}} ) \bigr], \\& \xi_{11} = \frac{1}{4} \bigl[\bar{f}_{\bar{x}\bar{x}}+ \bar{f}_{\bar {y}\bar{y}}+i (\bar{g}_{\bar{x}\bar{x}}+\bar{g}_{\bar{y}\bar {y}} ) \bigr], \\& \xi_{02} = \frac{1}{8} \bigl[\bar{f}_{\bar{x}\bar{x}}- \bar{f}_{\bar {y}\bar{y}}-2\bar{g}_{\bar{x}\bar{y}}+i (\bar{g}_{\bar{x}\bar {x}}- \bar{g}_{\bar{y}\bar{y}}+2\bar{f}_{\bar{x}\bar{y}} ) \bigr], \\& \xi_{21} = \frac{1}{16} \bigl[\bar{f}_{\bar{x}\bar{x}\bar{x}}+\bar {f}_{\bar{x}\bar{y}\bar{y}}+\bar{g}_{\bar{x}\bar{x}\bar{y}}+i (\bar {g}_{\bar{x}\bar{x}\bar{x}}+ \bar{g}_{\bar{x}\bar{y}\bar{y}}-2\bar {f}_{\bar{y}\bar{y}\bar{y}} ) \bigr]. \end{aligned}$$

Therefore, by the above analysis and the theorem presented as Section 3.4 in Guckenheimer [15], we obtain the following result.

Theorem 3.2

System (1.2) undergoes a Neimark-Sacker bifurcation at equilibrium $E_{1}$ if the conditions $(a, r, m, k, b, c)\in H_{B}$, $\theta\neq0$ hold and $\bar{r}^{*}$ varies in a small vicinity of the origin. In addition, if $\theta< 0$ (or $\theta> 0$), then an attracting (or repelling) invariant closed curve bifurcates from $E_{1}$ for $\bar{r}^{*} > 0$ (or $\bar{r}^{*} < 0$).

4 Numerical simulation

In this section, we draw the bifurcation diagrams of model (1.2) to validate the previous theoretical analysis and show the new interesting complex dynamical behaviors by using numerical simulations.

The bifurcation parameters are considered in the following two cases.

Case 1: Varying r in range $1.8\leq r \leq2.8$ and fixing $k=5$, $c=0.2$, $m=0.00001$, $b=100$, $a=0.001$. From the above data, we find that model (1.2) has only one positive fixed point $(6.8196,7.3943)$ and $(r_{1},a,m,k,c,b)=(2.079,0.001,0.00001,5,0.2,100)\in F_{B}$. By calculating, we see that the flip bifurcation emerges from this equilibrium at $r_{1}=2.079$ with $\alpha_{1}=-0.463538$ and $\alpha _{2}=0.014119$. This verifies Theorem 3.1.

From Figure 1(a), we see that the stability of a fixed point $E_{1}$ happens for $r<2.079$ and loses its stability at the flip bifurcation parameter value $r=2.079$. We also observe that there is a cascade of period doubling. Moreover, a chaotic set emerges with the increase of r.

Figure 1

Numerical simulations of equation ( 1.2 ). (a) Flip bifurcation of model (1.2) for $k=5$, $c=0.2$, $m=0.00001$, $b=100$, $a=0.001$; (b) NS bifurcation of model (1.2) for $k=5$, $c=0.2$, $m=0.00001$, $b=100$, $a=0.008$. Initial value $(x_{0},y_{0})=(5.0,2.5)$.

Case 2: We fix $k=5$, $c=0.2$, $m=0.00001$, $b=100$, $a=0.008$, and let the parameter r vary in the range $[3.475,3.495]$. By calculating, we see that model (1.2) has a unique endemic equilibrium $(2.5089,2.0006)$ and $(r_{2},a,m,b,c,k)=(3.482,0.008,0.00001,100,0.2,5)\in H_{B}$. By a direct calculation, we know that the Neimark-Sacker bifurcation occurs at $r_{2}=3.482$ with $l=1.353496>0$ and $\theta=-8.147308<0$. This shows the correctness of Theorem 3.2.

The bifurcation diagram shown in Figure 1(b) demonstrates that the stability of a fixed point $(2.5089,2.0006)$ happens for $r<3.482$ and loses its stability at $r=3.482$, and an attracting invariant closed curve appears if $r>3.482$.

In order to accurately catch the dynamic behavior of system (1.2) with respect to two parameters, the division diagrams are drawn in Figures 2-4. These parameters are combined arbitrarily. Lyapunov exponents are calculated in order to perform some numerical experiments. Figure 2 displays the division of this region by a pseudo-colored map, in which different attractors are identified by different colors. Figures 2 and 3 show a global view of the $a\times c$ or $c\times r$ parameter planes, respectively. For each point $(a,c)$ and $(c,r)$, the orbit corresponding initial condition successively converges to a chaotic, antiperiodic, or periodic attractor. In Figure 2, two parameters a and c are varied synchronously, while the other parameters are fixed as follows: $k=5$, $r=0.2$, $m=0.00001$, $b=100$. Then we can see that the periodic motions of the system appear alternatively. In Figure 3, we change two parameters c and r simultaneously, while the other parameters are selected as follows: $k=5$, $a=0.008$, $m=0.00001$, $b=100$. By observing the diagrams of two parameters, we can find that different color intensities represent different periodic values, the color bar on the right-hand side specifies the correspondence between colors and periods. Meanwhile, the periodic solutions are plotted in different colors and marked by the corresponding numbers (e.g., the number 4 represents period-4, painted in yellow; the number 9 represents period-9, painted in red; and the chaos phase appears when the number is equal to 30, painted in black).

Figure 2

Attractor region of different periods. (a) The 2D projection onto the $(a,c)$ plane of equation (1.2); (b) The magnification of (a). In this and further similar plots, different attractors are identified by different colors and numbers. For the sake of beauty, the period-1 is shown in white for (b).

Figure 3

Attractor region of different periods. (a) The $c\times r$ parameter plane of equation (1.2); (b) Local amplification of (a) for $c\in[0,0.2299]$.

In Figure 2, it is an interesting phenomenon that there is a chaotic region embedded in period regions. That is to say, when the searching efficiency a is changed from 0 to 0.001, the lower bound c will appear with abundant dynamic behavior accordingly. The cycle behavior appears alternately in the chaotic area, which can be clearly seen from (b). These cycle windows demonstrate the characteristics of the period-doubling and periodic-adding.

From Figure 3, we can see that there exists a chaotic region filled in period-1 area, which is quite different from what is seen in Figure 2. Compared to Figure 2, there are many periodicities between 2 and 26 embedded in period-1 regions, and the detailed periodic orbits can be found from (b). It can be explained that the hosts’ intrinsic rate r has a very rich behavior when the lower bound c is changed from 0 to 0.2, as shown in Figure 3. Therefore, in this discrete system the lower bound plays an important role. The frequent occurrence of crisis is seen when c increases from 0 to 2. In addition, these cycle windows illuminate the characteristics of the period-doubling and periodic-adding.

Furthermore, we can also observe abundant dynamic behaviors when more parameters vary at the same time, such as a and r, c and k. The corresponding phase diagrams are shown in Figure 4. In these two plots, there also exists a chaotic region embedded in period areas.

Figure 4

Numerical simulations of equation ( 1.2 ) with two parameters. (a) The $a\times r$ parameter plane of system (1.2); (b) The phase diagram for $c\times k$.

From Figure 1, we can observe the flip bifurcation and the Neimark-Sacker bifurcation when the hosts’ intrinsic rate r has different values, which is entirely consistent with the results in Theorems 3.1 and 3.2. Numerical results illustrate that the theoretical analysis is effective. The results obtained in this paper are interesting and useful and may be applied to a variety of the life systems.

5 Conclusion

The dynamic behavior of a host-parasitoid system is examined. We algebraically show that system (1.2) undergoes a bifurcation (flip or Neimark-Sacker) at a unique positive fixed point if r varies around the sets $F_{B}$ and $H_{B}$. Many forms of complexities, such as the cascade of period-doubling bifurcation and Neimark-Sacker bifurcation, are observed by numerical simulation. The numerical simulation results not only show consistency with the theoretical analysis but also exhibit unpredictable behaviors of the system through the bifurcation. Furthermore, we present the system’s parameter spaces which are restricted by two control parameters. Regardless of the combination of parameters, the different phase diagrams indicate that a variety of extraordinary geometries exist in the parameter plane. It can be told from our analysis that there exists a chaotic region in the periodic areas. Now, the new results reported in this paper give us more complete understanding of the discrete host-parasitoid system. However, it is still a challenging problem to explore a multiple parameter bifurcation in the system. We expect to obtain some more analytical results on this issue in the future.