Global dynamics and bifurcation analysis of a host–parasitoid model with strong Allee effect

ABSTRACT In this paper, we study the global dynamics and bifurcations of a two-dimensional discrete time host–parasitoid model with strong Allee effect. The existence of fixed points and their stability are analysed in all allowed parametric region. The bifurcation analysis shows that the model can undergo fold bifurcation and Neimark–Sacker bifurcation. As the parameters vary in a small neighbourhood of the Neimark–Sacker bifurcation condition, the unique positive fixed point changes its stability and an invariant closed circle bifurcates from the positive fixed point. From the viewpoint of biology, the invariant closed curve corresponds to the periodic or quasi-periodic oscillations between host and parasitoid populations. Furthermore, it is proved that all solutions of this model are bounded, and there exist some values of the parameters such that the model has a global attractor. These theoretical results reveal the complex dynamics of the present model.


Introduction
Ecological models exhibit a broad spectrum of non-fixed point dynamics ranging from characteristic natural cycles to more complex chaotic oscillations [1,2,[15][16][17]. Among various ecological models, the study of host-parasitoid model attracts more and more interest [9,10]. Parasitoids are insect species whose larva develop by feeding on the bodies of other arthropods, usually killing them. Larva emerge from the host and develop into free-living adults. The adults then lay their eggs in a subsequent generation of hosts. Most parasitoid larva require a specific life stage of the host, so parasitoid and host generations are linked to each other. Consequently, the interaction between host and parasitoid is often modelled by using a discrete-time step corresponding to the common generation length of host and parasitoid [3,6,8]. The following classical Nicholson-Bailey model is a famous example for host-parasitoid dynamics [16,18]: H n+1 = rH n e −γ P n , P n+1 = sH n (1 − e −γ P n ), where H n and P n represent the population of the host and parasitoid at year n, respectively, r is the number of offspring of an unparasitized host surviving to the next year. Assume the random search of parasitoids for host with parameter γ , then the probability that each host escapes parasitism from a single parasitoid is e −γ , and from parasitism is e −γ P n . Therefore, the probability of a host to become infected is 1 − e −γ P n . The parameter s describes the number of parasitoids that hatch from a single infected host. However, it was realized then that the model is not realistic, and the parasitoid population will go extinct while the host population will either go extinct or explodes with long time. Beddington et al. [2] observed that the main defect of the model (1) is that the host population is only regulated by the parasitoid. One method to remedy the situation is to consider both parasitism and the strong Allee effect on the host. The Allee effect refers to the phenomenon that the per capita birth rate declines at low population density [4,19], and the strong Allee effect refers to a phenomenon that a population exhibits a 'critical density', below this density the population declines to extinction and above the density it may increase (see [7,11,19]).
Recently, Livadiotis et al. [14] have investigated a discrete time host-parasitoid model which is an extension of the Nicholson-Bailey model by introducing a strong Allee effect in the dynamics of host [5,18]. The model of Livadiotis et al. [14] is where a, r, s, γ are positive parameters. For more details about its construction and biological meanings of the parameters, we refer the reader to [14]. For model (2), Livadiotis et al. [14] have studied the existence and local dynamics of fixed points, and observed that there are two scenarios, the first with no interior fixed point and the second with one interior fixed point. In the first scenario, they showed that either both host and parasitoid go extinction or there are two regions, an extinction region where both species go extinction and an exclusion region in which host survives and tends to carrying capacity. However, in the second scenario, it is possible that both host and parasitoid go extinction or there are two regions, an extinction region where both species go extinction and a coexistence region where both species survive. The aim of this paper is to study the global dynamics and bifurcations of system (2) from the viewpoint of mathematics. First, we make the following re-scaling transformation: x n a , P n = y n γ , then system (2) becomes Moreover, for the simplicity in mathematics, we assume that a = sγ , and the above system reduces into the following form: where α = r > 0, β = 1/a 2 > 0.
As we mentioned, Livadiotis et al. [14] have done substantial work on the existence and local dynamics of fixed points for the original model (2). In the present paper, we will carry forward the investigations of model (3) from three aspects. Firstly, we provide a complete division of the parameter plane (α, β), according to the number of fixed points. Based on the division, the existence of fixed points and their stability are analysed in each parameter region. Although model (3) is a special case with a = sγ , our analysis shows that the allowable number of fixed points and their possible topological types for model (3) are the same as that for the original model (2). Secondly, our bifurcation analysis shows that there are some values of the parameters such that the model can undergo fold bifurcation and Neimark-Sacker bifurcation, respectively. Thirdly, we discuss the boundedness of all solutions of this model as well as the global attractors in certain parameter regions.
The rest of the paper is organized as follows. In Section 2, we study the local dynamics of system (3) such as the existence and local stability of equilibria in R 2 + for system (3) in all allowed parameter region. In Section 3, we focus on the bifurcations of fixed points for system (3) in R 2 + . It is shown that system (3) can undergo fold bifurcation and Neimark-Sacker bifurcation in R 2 + . In Section 4, we prove that the solutions of system (3) are eventually bounded in the first quadrant R 2 + . The global attractors of system (3) are also analysed for certain parameter regions. In Section 5, some numerical simulations and corresponding biological interpretations are carried out to support our theoretical results. A brief conclusion is given in the last section.

Local dynamics of system (3)
In this section, we will study the existence of fixed points of system (3) and their local stability in the first quadrant R 2 + of (x, y). The results about the existence of fixed points are summarized into the following lemma.
, then D is divided into six sub-regions (I)-(VI) by these three curves (see Figure 1), and fixed points of system (3) are shown in the following statements.  (3) has three fixed points if some conditions are satisfied. More precisely,  Proof: To find fixed points of system (3) in R 2 + , we have to solve the following equations: for (x, y) ∈ R 2 + . We first find the boundary fixed points. We divide two cases to consider: (1) If x = 0, then the first equation of system (4) identically satisfied and from the second equation we get y = 0. Hence, system (3) has always the boundary fixed point O(0, 0) for all α > 0 and β > 0. (2) If y = 0, then the second equation of system (4) holds identically and from the first equation of system (4), we have the solutions (x * 1 , 0) and (x * 2 , 0) of (4), where We now classify the existence of positive x * 1 and x * 2 according to the sign of 1 − 4β(1 − α) to obtain the number of boundary fixed points for system (3).
On the other hand, we consider the existence of positive fixed point of system (3) in the interior of R 2 + . For this x = 0, then system (4) takes the following form: By eliminating x from system (6), we get α = − ye y e y − 1 + e y 1 + β ye y e y − 1 2 .
Let us denote F(y) = −ye y /(e y − 1) + e y (1 + β(ye y /(e y − 1)) 2 ). Then the y coordinates of positive fixed points of system (3) are the corresponding y coordinates of intersection points of the curve Z = F(y) and the line Z = α with y > 0.
Since F(y) is continuous and differentiable as y > 0, and for all y > 0 and β > 0, we calculate the derivative of F(y), Hence, if β ≥ α > 0, then there does not exist any intersection point of the curve Z = F(y) and the line Z = α with y > 0, which implies that system (3) has no positive fixed points in the interior of R 2 + . And if 0 < β < α, then there exists a unique intersection point (α, l) of the curve Z = F(y) and the line Z = α with l > 0 (see Figure 2). Therefore, system (3) has positive fixed points if and only if 0 < β < α. And the positive fixed point of system (3) is unique, we denote it by P(le l /e l − 1, l), where l is the positive solution of (7).
Note that the line α = β tangents to the curve α = 1 − 1/4β at a point ( 1 2 , 1 2 ) in the αβplane and this curve is completely under the line α = β (see Figure 1). By summarizing the above analysis, we can combine the cases for boundary fixed points and the positive fixed points (see Figure 1) to divide the parametric region into six sub-regions such that the number of these fixed points corresponding to their parametric region can be indicated in Table 1.
Hence, we obtain that system (3) has at least one and at most four fixed points in the first quadrant R 2 + for all α > 0 and β > 0. The proof is completed.
We now study the local stability of the fixed points (x 0 , y 0 ) of system (3) in R 2 + . The Jacobian matrix of system (3) at (x 0 , y 0 ) is Clearly, the characteristic equation of J G (x 0 , y 0 ) is where In order to study stability of the fixed point, we have to know the modulus of eigenvalues of the fixed point (x 0 , y 0 ), which are the modulus of roots of the characteristic equation (9), For the convenience of ader, we state the following lemma according to the relations between roots and coefficients of the quadratic equation.
We recall some definitions of topological types for a fixed point Depending on the modulus of eigenvalues of Equation (9), we state and prove the topological classification for the boundary fixed points as follows:

Theorem 2.3:
(1) In region I, both O(0, 0) and A 1 (x * 1 , 0) are saddle; Proof: Note that the eigenvalues of Jacobian matrix at O(0, 0) are λ 1 = α, λ 2 = 0. If x = 0 and y = 0, the Jacobian matrix (8) takes the following form: The eigenvalues of this matrix are For the fixed point O(0, 0), the Jacobian matrix (8) reduces into the following form: The eigenvalues for this matrix are , which exists under the condition α ≥ 1 − 1/4β, the eigenvalues are In fact, on the one hand, it is easy to verify that λ 1 < 1, on the other hand, which is obviously true under the condition α > 1 − 1/4β.
The characteristic equation of the Jacobin matrix J(x 0 , y 0 ) of system (3) about the unique positive fixed point P(le l /e l − 1, l) is given by  An illustration of the results of Theorem 2.4 is given in Figure 3, in which the stability of P(le l /(e l − 1), l) are shown in the parameter region {(α, β) : α > 0, β > 0}. Specifically, the red curve corresponds to the condition α = (l + 1)le l /(e l − 1 + 2l − le l ), and the blue curve corresponds to the condition −1 − q = p. And it should be mentioned that the condition p = 1+q does not hold in the first quadrant {(α, β) : α > 0, β > 0}, by the numerical computation.
By summarizing the results of Lemma 2.1, Theorems 2.3 and 2.4, we show that the allowable number of fixed points and their possible topological types for system (3) are the same as the original system (2), although system (3) is a special case with the assumption a = sγ .

Bifurcations of system (3)
The section is dedicated to the bifurcation analysis of system (3). According to the conclusions in Theorem 2.3 and 2.4, the number of fixed points and their topological types change as the parameters (α, β) vary in different regions of D = {(α, β) : α > 0, β > 0}. First of all, the number of boundary equilibrium points changes from one to three when the parameters (α, β) come across the curve C 1 , and we can obtain the following theorem of fold bifurcation [12].

Theorem 3.1: System (3) undergoes a fold bifurcation as the parameters (α, β) go through the curve
Proof: For β > 1 4 , denoteα = 1 − 1/4β, since y = 0 is invariant with respect to the system (3). Thus, we can restrict system (3) on the line y = 0 to determine the bifurcation, where system (3) becomes and the corresponding map is f (x) = (x 2 + αx)/(1 + βx 2 ). Now at α =α andx = 1/2β, we have This is the non-hyperbolic condition for fold bifurcation. Now, we will check the nondegenerate condition for the map, that is, Thus, the map satisfies the non-degenerate condition of fold bifurcation. Therefore, system (3) undergoes a fold bifurcation when the parameters (α, β) come across the curve C 1 .
In order to guarantee the Neimark-Sacker bifurcation for (16), we require that the following discriminatory quantity is not zero (see [12,13] where After calculating, we obtain Based on the above analysis and the Neimark-Sacker bifurcation theorem discussed in [12], we finally arrive at the following theorem.
According to the Neimark-Sacker bifurcation theorem in [12], the bifurcation is called supercritical and an attracting invariant closed curve bifurcates from the fixed point if the discriminatory quantity < 0. In Section 5, our numerical simulation shows that a supercritical Neimark-Sacker bifurcation can occur for system (11). From the viewpoint of biology, the attracting invariant closed curve implies that the host and parasitoid populations will coexist under the periodic or quasi-periodic oscillations with long time.

Boundedness of solutions and global attractors
In this section, we prove the boundedness of solutions of system (3) in the first quadrant R 2 + and discuss if some fixed points of system (3) becomes attractors in the first quadrant R 2 + . We first recall some concepts. A solution {x n , y n } of system (3)   From the proof of lemma 4.1, we find that if 0 < α < 1 − 1/4β and β > 1 4 , then every solution {(x n , y n )} of system (3) in the first quadrant R 2 + satisfies which leads that lim n→+∞ x n = 0 and lim n→+∞ y n = 0. Hence, the fixed point O(0, 0) is a global attractor of system (3) in the first quadrant R 2 + . (ii) On the curve C 1 , O(0, 0) is a sink and A(1/2β, 0) is non-hyperbolic by the conclusion (7) of Theorem 2.3. The solution {(x n , y n )} of system (3) satisfies Thus, we have lim n→+∞ x n = 0 for all solutions of system (3) with 0 < x n < 1/2β for some n, and lim n→+∞ x n = 1/2β for all other solutions. On the other hand, since the series {y n } is bounded by Lemma 4.1 and y = 0 is the unique fixed point for the second equation of system (3), one can obtain that lim n→+∞ y n = 0.
(iii) If (α, β) ∈ IV, then O(0, 0) is a sink, A 1 (x * 1 , 0) is a sink, and A 2 (x * 2 , 0) is a saddle by the conclusion (4) of Theorem 2.3. Moreover, the solution {(x n , y n )} of system (3) in the first quadrant R 2 + satisfies Hence, lim n→+∞ x n = 0 for all solutions of system (3) with 0 < x n < x * 2 for some n, and lim n→+∞ x n = x * 1 for all other solutions. By similar arguments, we have lim n→+∞ y n = 0. The proof is completed.
An illustration of the results of Theorem 4.2 is given in Figure 5, in which the global attractors of system (3) are shown in the parameter regions V I, C 1 and IV.

Numerical simulations
In this section, we give some numerical simulations for system (3) to support our theoretical results. We will simulate the phase portraits of system (3) for different parameter regions.
If (α, β) ∈ V, system (3) possesses three equilibrium points, where O(0, 0) is a sink, A 1 (x * 1 , 0) is a saddle and A 2 (x * 2 , 0) is a source. An example of the phase portraits of system (3) for (α, β) ∈ V is given in Figure 7(a), and all the orbits tend to the point O(0, 0), that is, both host and parasitoid will go extinction with long time.
If (α, β) ∈ IV, system (3) possesses three equilibrium points, where O(0, 0) and A 1 (x * 1 , 0) are sink, and A 2 (x * 2 , 0) is a saddle. Figure 7(b) presents an example of the phase portraits of system (3) for (α, β) ∈ IV, from which one can see that some orbits tend to O(0, 0), while other orbits tend to A 1 (x * 1 , 0). In the context of biology, the population of parasitoid will go extinction with long time, but the long-time behaviour of the host depends on the initial conditions. More specifically, the population of the host will survive if the initial population is in a certain region; otherwise, it will also go extinction.  If (α, β) ∈ II, system (3) possesses two equilibrium points, where O(0, 0) is a saddle and A 1 (x * 1 , 0) is a sink. An example of the phase portraits of system (3) for (α, β) ∈ II is given in Figure 7(c), in which all the orbits tend to the point A 1 (x * 1 , 0). In this case, the population of parasitoid will go extinction while the host will survive with long time.
If (α, β) ∈ I ∪ III ∪ C 21 , system (3) possesses a unique positive equilibrium point P(le l /(e l − 1), l) besides the boundary points. According to the conclusions of Theorem 2.4, the topological classification of P changes with different values of (α, β). For instance, Figure 7(d,e) gives two examples of system (3) for (α, β) ∈ III, where O(0, 0) is a sink and P is a stable focus. From the viewpoint of biology, both the host and the parasitoid will go extinction for some initial condition, or they will coexist for other initial conditions. Figure 7(f) gives an example of system (3) for (α, β) ∈ III, where O(0, 0) is a sink while P is a source.
If (α, β) ∈ I, system (3) possesses equilibrium points O(0, 0), A 1 (x * 1 , 0) and a unique positive equilibrium point P, where both O(0, 0) and A 1 (x * 1 , 0) are saddle. Figure 8(a) gives an example of phase orbits of system (3), where the point P is a sink. And Figure 8(b) presents an example for which P is a stable focus. Moreover, our numerical simulation shows that a supercritical Neimark-Sacker bifurcation can occur for the system with (α, β) ∈ I. Choose β = 0.65, then we can obtain the bifurcation valueα = 5.71755 of α according to the Neimark-Sacker condition H P . If α <α, as in Figure 8(c, d), P is a stable focus. If α >α, as in Figure 8(e)-(l), the equilibrium point P losses its stability, and an attracting invariant closed curve bifurcates from the fixed point. Our simulation indicates that the closed curve is attractive from both outside and inside (see Figure 8(h, i)). From the viewpoint of biology, the attracting invariant closed curve implies that the host and parasitoid populations will coexist under the periodic or quasi-periodic oscillations with long time. More simulations of the phase portraits for system (3) are shown in Figure 9.

Conclusions
In this paper, we studied the global dynamics and bifurcations of a two-dimensional discrete time host-parasitoid model with strong Allee effect. By a re-scaling transformation and an assumption a = sγ , we reduced the model of Livadiotis et al. [14] to the present model (3). Based on the number of fixed points, the parameter region {(α, β) : α > 0, β > 0} is completely divided into six regions (I)-(VI), and in each parameter region the existence of fixed points of system (3) and their stability are analysed. Compared with the original model (2), the derived model (3) has the same allowable number of fixed points and the same possible topological types. In general, system (3) has at least one fixed point and at most four fixed points. From the viewpoint of biology: (i) if (α, β) ∈ V&VI, then the populations of both host and parasitoid will go extinction with long time; (ii) if (α, β) ∈ IV, then the population of parasitoid will go extinction with long time, but the long time behaviour of the host depends on the initial conditions, more specific, the population of the host will tend to a steady state if the initial populations are in a certain region, otherwise, it will also go to extinction; (iii) if (α, β) ∈ II, then the population of parasitoid will go extinction while the host will tend to a steady state with long time; (iv) if (α, β) ∈ I&III, then both the host and parasitoid will go extinction for some initial condition, or they will coexist for other initial conditions. Moreover, our bifurcation analysis shows that there are some values of the parameters such that the model can undergo fold bifurcation and Neimark-Sacker bifurcation, respectively. In the context of biology, the attracting invariant closed curve bifurcated from Neimark-Sacker bifurcation implies that the host and parasitoid populations will coexist under the periodic or quasi-periodic oscillations with long time. At the end, we discuss the boundedness of all solutions of this model as well as the global attractors in certain parameter regions.