1. Introduction

With the increasing application of ecological models, the host-parasitoid models have been extensively explored. The host-parasitoid models are of great significance among the relationships between the biotic populations. Although this kind of model is investigated by many scholars (see [1] [2] [3] [4] ), the research about discrete systems is relatively few. Compared to the continuous ones, the dynamics of discrete systems are more interesting. In our real life, most practical problems are described by the discrete systems, and it is necessary for us to discretize the continuous systems. And the dynamics of the discrete-time models can present a much richer set of patterns than those observed in continuous-time models, these models can lead to unpredictable dynamics from a biological point of view. So the study of discrete host-parasitoid system is very important. Research on host-parasitoid system, indicates this kind of model can have very complex dynamics (see [5] [6] [7] [8] ).

Xu and Boyce in [9] have investigated a host-parasitoid model which is influenced by the parasitoid interference, and they have found that the density dependence decides the mutual interference of parasitoid. In [10] , Beddington describes the dynamics of a parasitoid-host-parasitoid model, whose populations is non-overlapping generations. He illustrates the model having dynamical behaviors that are closely analogous to those observed in the first-order situation. Using numerical simulation, the authors in [11] also study the chaotic dynamical behavior of the parasitoid-host-parasitoid model, which shows that the advantage coefficient can stabilize the dynamics. The mathematical equation of the model in [10] [11] is proposed as a discrete-time model

$\{\begin{array}{l}H\left(t+1\right)=H\left(t\right)\exp \left[r\left(1-\frac{H\left(t\right)}{K}\right)-a{\left[P\left(t\right)\right]}^{1-m}\right],\\ P\left(t+1\right)=H\left(t\right)\left[1-\exp \left(-a{\left[P\left(t\right)\right]}^{1-m}\right)\right],\end{array}$ (1)

where $H\left(t\right)$ is the host population size in generation t, $P\left(t\right)$ is the parasitoid population size in generation t. The constant a indicates the searching efficiency and m is the interference coefficient. The host grows logistically with the carrying capacity K and the intrinsic growth rate r.

For simplicity, we rewrite the system (1) as the following:

$\{\begin{array}{l}x\to xexp\left[r\left(1-\frac{x}{K}\right)-a{y}^{1-m}\right]\mathrm{,}\\ y\to x\left[1-exp\left(-a{y}^{1-m}\right)\right]\mathrm{.}\end{array}$ (2)

In our paper, the parasitoid-host-parasitoid system (2) is investigated in further details. We mainly focus on its bifurcations and possible chaos qualitatively. Based on the center manifold theorem and bifurcation theory (see [12] [13] ), we can obtain the detailed existence conditions of these bifurcations. Numerical simulations, including bifurcation diagrams, phase portraits, are used to verify theoretical analysis. The results obtained in the paper can be regarded as the beneficial supplement of the work in [10] [11] .

The layout of this paper is organized as follows: The existence and stability criterion of the equilibria of system (2) are presented in Section 2; Section 3 deals with the flip bifurcation and Neimark-Sacker bifurcation, and derives the existence conditions of the bifurcations; Numerical simulations using MATLAB are presented in Section 4 to illustrate the theoretical results; A brief discussion is carried out in Section 5.

2. Stability of equilibria

In order to obtain the equilibria of system (2), we need to use the mathematical software. With the aid of Maple program, we get the following three equilibria $E_0\left(\mathrm{0,0}\right)$ , $E_1\left(K\mathrm{,0}\right)$ , $E_2\left(x_2\mathrm{,}{y}_2\right)$ , where $x_2$ , $y_2$ satisfy the following equation:

$\{\begin{array}{l}{x}_2=x_2\exp \left[r\left(1-\frac{x_2}{K}\right)-a{\left(y_2\right)}^{1-m}\right],\\ y_2=x_2\left[1-\exp \left(-a{\left(y_2\right)}^{1-m}\right)\right].\end{array}$ (3)

The qualitative behavior of the system (2) will be investigated. The local dynamics of the system (2) near a fixed point depends on its Jacobian matrix. The Jacobian matrix at the state variable is given by

$J\left(x,y\right)=\left(\begin{array}{cc}{\text{e}}^{r-\frac{rx}{K}-a{y}^{1-m}}-\frac{rx}{K}{\text{e}}^{r-\frac{rx}{K}-a{y}^{1-m}}& -a\left(1-m\right)x{y}^{-m}{\text{e}}^{r-\frac{rx}{K}-a{y}^{1-m}}\\ 1-{\text{e}}^{-a{y}^{1-m}}& a\left(1-m\right)x{y}^{-m}{\text{e}}^{-a{y}^{1-m}}\end{array}\right),$

According to J, one can obtain two eigenvalues ${\lambda }_1=0,\mathrm{<mi>\; </mi>}{\lambda }_2={\text{e}}^r$ with given $E_0\left(\mathrm{0,0}\right)$ . From which, one can easily check that $E_0\left(\mathrm{0,0}\right)$ is a stable node $\left(\left|{\text{e}}^r\right|<1\right)$ . Two eigenvalues at the equilibrium $E_1\left(K\mathrm{,0}\right)$ are ${\lambda }_1=0,\mathrm{<mi>\; </mi>}{\lambda }_2=1-r$ , then one can get that $E_1\left(K\mathrm{,0}\right)$ is also a stable node $\left(\left|1-r\right|<1\right)$ (see [14] ). Next, we only need to consider the stability of the equilibrium $E_2$ .

The following is Jacobian matrix of (3) at $E_2$ :

$J_2\left(x_2\mathrm{,}{y}_2\right)=\left(\begin{array}{cc}1-rL& -H\\ 1-G& GH\end{array}\right)\mathrm{,}$

The characteristic equation of matrix $J_2$ is

${\lambda }^2+P\left(x_2,y_2\right)\lambda +Q\left(x_2,y_2\right)=0,$ (4)

where

$P\left(x_2,y_2\right)=-\left(1-rL+GH\right),\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}Q\left(x_2,y_2\right)=H-rLGH,$

$L=\frac{x_2}{K},\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}G={\text{e}}^{-a{y}_2^{1-m}},\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}H=a\left(1-m\right)x_2{y}_2^{-m}.$

From (4), then we have

$F\left(1\right)=rL-GH-rGHL+H,$

$F\left(-1\right)=2-rL+GH-rGHL+H\mathrm{.}$

In order to disuss the stability of the fixed point $E_2$ , we also need the following Lemma, which can be easily found from the theorem presented in [14] .

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

i) $\left|{\lambda }_1\right|<1$ , $\left|{\lambda }_2\right|<1$ if and only if $F\left(-1\right)>0$ and $Q<1$ ;

ii) $\left|{\lambda }_1\right|<1$ , $\left|{\lambda }_2\right|>1$ (or $\left|{\lambda }_1\right|>1$ and $\left|{\lambda }_2\right|<1$ ) if and only if $F\left(-1\right)<0$ ;

iii) $\left|{\lambda }_1\right|>1$ , $\left|{\lambda }_2\right|>1$ if and only if $F\left(-1\right)>0$ and $Q>1$ ;

iv) ${\lambda }_1=-1$ , $\left|{\lambda }_2\right|\ne 1$ if and only if $F\left(-1\right)=0$ and $P\ne \mathrm{0,2}$ ;

v) ${\lambda }_1$ , ${\lambda }_2$ are complex and $\left|{\lambda }_1\right|=\left|{\lambda }_2\right|=1$ if and only if $P^2-4Q<0$ and $Q=1$ .

Let ${\lambda }_1$ and ${\lambda }_2$ be the roots of (2), which are called eigenvalues of the fixed point $E_2\left(x_2\mathrm{,\; <mi>\; </mi>}{y}_2\right)$ . The fixed point $\left(x_2\mathrm{,\; <mi>\; </mi>}{y}_2\right)$ is a sink or locally asymptotically stable if $\left|{\lambda }_1\right|<1,\mathrm{<mi>\; </mi>}\left|{\lambda }_2\right|<1$ . $E_2\left(x_2\mathrm{,\; <mi>\; </mi>}{y}_2\right)$ is a source or locally unstable if $\left|{\lambda }_1\right|>1,$ $\left|{\lambda }_2\right|>1$ . $E_2\left(x_2\mathrm{,\; <mi>\; </mi>}{y}_2\right)$ is a saddle if $\left|{\lambda }_1\right|<1$ and $\left|{\lambda }_2\right|>1$ (or $\left|{\lambda }_1\right|>1$ and $\left|{\lambda }_2\right|<1$ ). The fixed point $\left(x_2\mathrm{,\; <mi>\; </mi>}{y}_2\right)$ is non-hyperbolic if either $\left|{\lambda }_1\right|=1$ or $\left|{\lambda }_2\right|=1$ .

From Lemma 2.1, we state the following theorem:

Theorem 2.1. For the positive equilibrium $E_2$ , we have the following estimates:

i) $E_2\left(x_2\mathrm{,\; <mi>\; </mi>}{y}_2\right)$ is a sink if the condition hods: $r<\frac{2+GH+H}{L+GHL}$ and $r>\frac{H-1}{GHL}$ ;

ii) $E_2\left(x_2\mathrm{,\; <mi>\; </mi>}{y}_2\right)$ is a source if the condition holds: $r<\frac{2+GH+H}{L+GHL}$ and $r<\frac{H-1}{GHL}$ ;

iii) $E_2\left(x_2\mathrm{,\; <mi>\; </mi>}{y}_2\right)$ is a saddle if the condition holds: $r>\frac{2+GH+H}{L+GHL}$ ;

iv) $E_2\left(x_2\mathrm{,\; <mi>\; </mi>}{y}_2\right)$ is non-hyperbolic if either condition (iv.1) or (iv.2) holds:

iv.1) $r=\frac{2+GH+H}{L+GHL}$ and $r\ne \frac{1+GH}{L}\mathrm{,}r\ne \frac{3+GH}{L}$ ;

iv.2) $r=\frac{H-1}{GHL}$ and $\frac{GH-1}{L}<r<\frac{GH+3}{L}$ .

From the above conclusion, if the term (iv.1) of Theorem 2.1 holds, one can easily find that one of the eigenvalues of $E_2\left(x_2\mathrm{,\; <mi>\; </mi>}{y}_2\right)$ is -1 and the other is $2+GH-rL$ , which is neither 1 nor -1. If the term (iv.2) of Theorem 2.1 holds, then the eigenvalues of $E_2$ are a pair of complex conjugate numbers whose modulus is 1.

Let

The equilibrium can arise flip bifurcation when parameters change in a small neighborhood of.

Let

The equilibrium can lead to the bifurcation of Neimark-Sacker (discrete Hopf bifurcation) when parameters are restricted to a small scope of.

3. Bifurcation

Now we mainly center on bifurcations (see [15] [16] [17] ) of system (2). In the following research, r is chosen as a bifurcation parameter.

3.1. Flip bifurcation

Select arbitrary parameters from, the system (2) is replaced by

(5)

The system (5) has a unique positive fixed point, the responding eigenvalues are, with. Since, , choosing as a bifurcation parameter, we reconsider the map (5) as given belowing:

(6)

where, which is a small perturbation parameter of.

Let, we transform the fixed point to the origin, then we have

(7)

where

(8)

We can build an invertible matrix T

Consider the following translation:

then the system (7) becomes

(9)

where

and

By the center manifold theorem in [12] , one can easily determine the center manifold of (9) at the origin, its expression is given as

By simple calculations, one can obtain

Thus, the map is restricted to the center manifold, which is given by

(10)

where

If the system (10) goes through a flip bifurcation, the following conditions must hold:, where

Based on the above analyses and using the bifurcation theorems presented in [13] , we obtain the following estimate:

Theorem 3.1. When the parameter varies in a small vicinity of the point, and, the system (2) undergoes a flip bifurcation at. Furthermore, the period-2 orbit bifurcated from this point is stable (unstable) while ().

3.2. Hopf bifurcation

We consider the following system by selecting parameters arbitrarily.

(11)

is a only positive fixed point of the system (11), where is given by (3) and

Choosing as a bifurcation parameter, we reconsider the system (11) as given below:

(12)

where, it is a small perturbation parameter of.

Let. After transforming point to the point, we have

(13)

where are given in (8), and.

The characteristic equation of map (13) at is given by

where

Since parameters, the eigenvalues of are a pair of complex conjugate numbers and with modulus 1, where

then, we have

Moreover, it is required, , that is to say, nondegeneracy condition, which is equivalent to. Notice that, thus,. We just need following conditions to be true, i.e.

(14)

So the eigenvalues do not lie in the intersection of the unit circle with the coordinate axes when and the conditions (14) hold.

Let, and construct an invertible matrix

Using the following translation:

then the system (13) becomes

(15)

where

and

Therefore

To assure that the map (13) passes though Neimark-Sacker bifurcation, we need to let the following discriminatory quantity is not zero:

(16)

where

According to the previous discussions and applying Hopf bifurcation theorems in [12] , we can get

Theorem 3.2. If the parameters, and varies in a limited region of the origin, the system (2) goes through a Neimark-Sacker bifurcation at. Furthermore, there exists a unique attracting (or repelling) invariant closed curve bifurcated from for (or) while (or).

4. Numerical simulation

We have gotten the theoretical results of system (2) based on the qualitative theory. In this section, we outline a numerical methods to validate the previous analysis and provide some numerical results by using MATLAB. We draw the diagrams for bifurcation and phase portraits to show new interesting complex

dynamical behaviors. The bifurcation parameters are considered in the following two cases:

Case 1: Varying r in range, and fixing, , with initial values of. From above data, one can easily obtain that the equilibrium is. By observing the bifurcation diagrams we can see that a flip bifurcation appears at. In this case,. It satisfies the Theorem 3.1.

From Figure 1, we can observe that the fixed point is stable for, loses its stability at and period doubling phenomena lead to chaos for.

The phase portraits in Figure 2 show that there are orbits of period 2, 4, 8 when. And chaotic sets can be seen when.

Case 2: We fix and let the parameter r vary in the range. By calculating, we know that the Neimark-Sacker bifurcation occurs at when, and its eigenvalues are . For, we have

It shows that the Theorem 3.2 holds.

From Figure 3, we observe that there exists a stable equilibrium for, a Hopf circle happens at. And an attracting invariant closed curve bifurcates from with increase of r. The phase portraits Figure 4 for different values of r illustrate that there appears a smooth invariant curve bifurcated from the stable fixed point, and its radius is getting larger with respect to the growth of r.

5. Conclusion

The dynamical behaviors of a parasitoid-host-parasitoid system are investigated. The theoretical analyses demonstrate that the system (2) can appear as flip bifurcation and Neimark-Sacker bifurcation. We present the numerical diagrams to validate analytical effectiveness. We also observe many forms of complexities from these diagrams, such as the cascade of period-doubling bifurcation and Neimark-Sacker bifurcation. Hence we can find that the discrete-time models have far richer dynamical behaviors as compared to continuous-time models. All these results are obtained by a simple system of only two maps and we believe that similar results can be achieved by more general systems.

Acknowledgements

We would like to thank the associate editor and the anonymous referees for their valuable comments and helpful suggestions, which have led to a great improvement of the initial version. We also wish to acknowledge conversations with Yandong Chu of the Department of Mathematics at Lanzhou Jiaotong University regarding this problem. This work is supported by the National Natural Science Foundation of China (No.11161027).