Original Articles

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

Abdul Qadeer Khan Department of Mathematics, University of Azad Jammu & Kashmir, Muzaffarabad, PakistanView further author information

Jiying Ma College of Science, University of Shanghai for Science and Technology, Shanghai, People's Republic of ChinaCorrespondencemajiying100@126.com
View further author information

Dongmei Xiao Department of Mathematics, Shanghai Jiao Tong University, Shanghai, People's Republic of ChinaView further author information

Original Articles

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

ABSTRACT
Formulae display:?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.

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.

KEYWORDS:

Host–parasitoid model Allee effect fixed points attractor bifurcations

2010 AMS MATHEMATICS SUBJECT CLASSIFICATIONS:

34K20 39A10 40A05 92C50

1. 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–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]: (1) $H_{n + 1} = r H_{n} e^{- γ P_{n}}, P_{n + 1} = s H_{n} (1 - e^{- γ P_{n}}),$ (1) 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 (2) $H_{n + 1} = \frac{a H_{n}^{2} + r H_{n}}{1 + H_{n}^{2}} e^{- γ P_{n}}, P_{n + 1} = s H_{n} (1 - e^{- γ P_{n}}),$ (2) 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: $H_{n} = \frac{x_{n}}{a}, P_{n} = \frac{y_{n}}{γ},$ then system (2) becomes $x_{n + 1} = \frac{x_{n}^{2} + r x_{n}}{1 + \frac{x_{n}^{2}}{a^{2}}} e^{- y_{n}}, y_{n + 1} = \frac{s γ}{a} x_{n} (1 - e^{- y_{n}}) .$ Moreover, for the simplicity in mathematics, we assume that $a = s γ$ , and the above system reduces into the following form: (3) $x_{n + 1} = \frac{x_{n}^{2} + α x_{n}}{1 + β x_{n}^{2}} e^{- y_{n}}, y_{n + 1} = x_{n} (1 - e^{- y_{n}}),$ (3) 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.

2. 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.

Lemma 2.1.

System (3) has at least one fixed point and at most four fixed points in the first quadrant $R_{+}^{2}$ for all $α > 0$ and $β > 0$ . More precisely, in the parameter region $D = {(α, β) : α > 0, β > 0}$ , define three curves $C_{1} = {(α, β) : β = 1 / 4 (1 - α), α < 1},$ $C_{2} = {(α, β) : α = 1, β > 0},$ and $C_{3} = {(α, β) : α = β, α > 0},$ 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.

System (3) has a unique fixed point at $O (0, 0)$ if $(α, β) \in V I,$ here $V I = {(α, β) : 0 < α < 1 - 1 / 4 β, β > \frac{1}{4}};$
System (3) has two fixed points if some conditions are satisfied. More precisely,
- (ii.1) System (3) has two fixed points $A (1 / 2 β, 0)$ and $O (0, 0)$ if $(α, β) \in C_{1};$
- (ii.2) System (3) has two fixed points $B (1 / β, 0)$ and $O (0, 0)$ if $(α, β) \in C_{2}$ with $β \geq 1;$
- (ii.3) System (3) has two fixed points $A_{1} (x_{1}^{*}, 0)$ and $O (0, 0)$ if $(α, β) \in I I,$ here $I I = {(α, β) : 1 < α \leq β},$ $x_{1} = 1 / 2 β + (1 / 2 β) \sqrt{1 - 4 β (1 - α)}$ .
System (3) has three fixed points if some conditions are satisfied. More precisely,
- (iii.1) System (3) has three fixed points $O (0, 0),$ $A_{1} (x_{1}^{*}, 0)$ and $A_{2} (x_{2}^{*}, 0)$ if $(α, β) \in I V,$ here $I V = {(α, β) : α \leq β < 1 / 4 (1 - α), \frac{1}{2} < α < 1},$ $x_{2}^{*} = 1 / 2 β - (1 / 2 β) \sqrt{1 - 4 β (1 - α)};$
- (iii.2) System (3) has three fixed points $O (0, 0),$ $A_{1} (x_{1}^{*}, 0)$ and $A_{2} (x_{2}^{*}, 0)$ if $(α, β) \in V,$ where $V = {(α, β) : α \leq β < 1 / 4 (1 - α), 0 < α < \frac{1}{2}};$
- (iii.3) System (3) has three fixed points $O (0, 0),$ $B (1 / β, 0)$ and $P (l e^{l} / (e^{l} - 1), l)$ if $(α, β) \in C_{2}$ with $0 < β < 1,$ where l is the positive solution of $α = - y e^{y} / (e^{y} - 1) + e^{y} (1 + β (y e^{y} / (e^{y} - 1))^{2});$
- (iii.4) System (3) has three fixed points $O (0, 0),$ $A_{1} (x_{1}^{*}, 0)$ and $P (l e^{l} / e^{l} - 1, l)$ if $(α, β) \in I,$ where $I = {(α, β) : α > β > 0, α > 1};$
System (3) has four fixed points $O (0, 0),$ $A_{1} (x_{1}^{*}, 0),$ $A_{2} (x_{2}^{*}, 0)$ and $P (l e^{l} / (e^{l} - 1), l)$ if $(α, β) \in I I I,$ where $I I I = {(α, β) : 0 < β < α < 1}$ .

Figure 1. The division of the parameter region ${(α, β) : α > 0, β > 0}$ , according to the number of fixed points of system (3).

Display full size

Proof.

To find fixed points of system (3) in $R_{+}^{2}$ , we have to solve the following equations: (4) $x = \frac{x^{2} + α x}{1 + β x^{2}} e^{- y}, y = x (1 - e^{- y}),$ (4) for $(x, y) \in R_{+}^{2}$ .

We first find the boundary fixed points. We divide two cases to consider:

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$ .
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 (5) $x_{1}^{*} = \frac{1}{2 β} + \frac{1}{2 β} \sqrt{1 - 4 β (1 - α)}, x_{2}^{*} = \frac{1}{2 β} - \frac{1}{2 β} \sqrt{1 - 4 β (1 - α)} .$ (5)

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).

If $1 - 4 β (1 - α) < 0$ , that is, $0 < α < 1 - 1 / 4 β$ and $β > \frac{1}{4}$ , then $x_{1}^{*}$ and $x_{2}^{*}$ are not real numbers. Thus, system (3) has a unique boundary fixed point $O (0, 0) .$
If $1 - 4 β (1 - α) = 0$ , that is, $α = 1 - 1 / 4 β > 0$ and $β > \frac{1}{4}$ , then $x_{1}^{*} = x_{2}^{*} = 1 / 2 β > 0$ . Hence, system (3) has two boundary fixed point $A (1 / 2 β, 0)$ and $O (0, 0)$ .
If $1 - 4 β (1 - α) > 0$ , that is, either $0 < β < 1 / 4 (1 - α)$ and $0 < α < 1$ , or $α \geq 1$ and $β > 0$ , then $x_{1}^{*}$ and $x_{2}^{*}$ are real numbers. Moreover,
- (iii.1) $x_{1}^{*} > 0$ and $x_{2}^{*} > 0$ when either $0 < β < 1 / 4 (1 - α)$ , $α \neq β$ and $0 < α < 1$ , or $α = β$ but $α \neq \frac{1}{2}$ . Therefore, system (3) has three boundary fixed points $A_{1} (x_{1}^{*}, 0)$ , $A_{2} (x_{2}^{*}, 0)$ and $O (0, 0)$ if either $0 < β < 1 / 4 (1 - α)$ , $α \neq β$ and $0 < α < 1$ , or $α = β$ but $α \neq \frac{1}{2}$ ;
- (iii.2) $x_{1}^{*} > 0$ and $x_{2}^{*} = 0$ when either $α = 1$ and $β > 0$ , or $α = β = \frac{1}{2}$ . Hence, system (3) has two boundary fixed points $B (1 / β, 0)$ and $O (0, 0)$ ( $A (1 / 2 β, 0)$ and $O (0, 0)$ ) if $α = 1$ and $β > 0$ ( $α = β = \frac{1}{2}$ , resp.);
- (iii.3) $x_{1}^{*} > 0$ and $x_{2}^{*} < 0$ when $α > 1$ and $β > 0$ . Thus, system (3) has two boundary fixed points $A_{1} (x_{1}^{*}, 0)$ and $O (0, 0)$ if $α > 1$ and $β > 0$ .

It is clear that system (3) has at least one and at most three boundary fixed points for all $α > 0$ and $β > 0$ .

On the other hand, we consider the existence of positive fixed point of system (3) in the interior of $R_{+}^{2}$ . For this $x \neq 0$ , then system (4) takes the following form: (6) $1 = \frac{x + α}{1 + β x^{2}} e^{- y}, y = x (1 - e^{- y}) .$ (6) By eliminating x from system (6), we get (7) $α = - \frac{y e^{y}}{e^{y} - 1} + e^{y} (1 + β {(\frac{y e^{y}}{e^{y} - 1})}^{2}) .$ (7) Let us denote $F (y) = - y e^{y} / (e^{y} - 1) + e^{y} (1 + β (y e^{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)$ , $\begin{aligned} F^{'} (y) & = - \frac{e^{y}}{e^{y} - 1} + \frac{y e^{2 y}}{(e^{y} - 1)^{2}} - \frac{y e^{y}}{e^{y} - 1} + e^{y} (1 + \frac{β y^{2} e^{2 y}}{(e^{y} - 1)^{2}}) \\ + e^{y} (\frac{2 β y e^{2 y}}{(e^{y} - 1)^{2}} - \frac{2 β y^{2} e^{3 y}}{(e^{y} - 1)^{3}} + \frac{2 β y^{2} e^{2 y}}{(e^{y} - 1)^{2}}) \\ = \frac{e^{y}}{(e^{y} - 1)^{2}} (e^{2 y} - 3 e^{y} + 2 + y) + \frac{β y e^{3 y}}{(e^{y} - 1)^{3}} (y e^{y} + 2 e^{y} - 2 - 3 y) > 0. \end{aligned}$

Moreover $\begin{aligned} lim_{y \to 0^{+}} F (y) & = lim_{y \to 0^{+}} (- \frac{y e^{y}}{e^{y} - 1} + e^{y} (1 + β {(\frac{y e^{y}}{e^{y} - 1})}^{2})), \\ = β . \end{aligned}$

Hence, if $β \geq α > 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 (l e^{l} / e^{l} - 1, l)$ , where l is the positive solution of (7).

Figure 2. $Z = F (y), Z = α$ : (a) $α < β$ , (b) $α = β$ and (c) $α > β$ .

Display full size

Note that the line $α = β$ tangents to the curve $α = 1 - 1 / 4 β$ at a point $(\frac{1}{2}, \frac{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.

Table 1. Number of fixed points in different parameter regions.

Display Table

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 (8) $J_{G} (x_{0}, y_{0}) = (\begin{array}{cc} (2 \frac{x_{0} + α}{1 + β x_{0}^{2}} - α) \frac{e^{- y_{0}}}{1 + β {x_{0}}^{2}} & - \frac{{x_{0}}^{2} + α x_{0}}{1 + β {x_{0}}^{2}} e^{- y_{0}} \\ [3 p t] 1 - e^{- y_{0}} & x_{0} e^{- y_{0}} \end{array}) .$ (8) Clearly, the characteristic equation of $J_{G} (x_{0}, y_{0})$ is (9) $λ^{2} - p (x_{0}, y_{0}) λ + q (x_{0}, y_{0}) = 0,$ (9) where $\begin{aligned} p (x_{0}, y_{0}) & = e^{- y_{0}} ((2 \frac{x_{0} + α}{1 + β {x_{0}}^{2}} - α) \frac{1}{1 + β {x_{0}}^{2}} + x_{0}), \\ q (x_{0}, y_{0}) & = \frac{e^{- y_{0}}}{1 + β {x_{0}}^{2}} ((2 \frac{x_{0} + α}{1 + β {x_{0}}^{2}} - α) x_{0} e^{- y_{0}} + ({x_{0}}^{2} + α x_{0}) (1 - e^{- y_{0}})) . \end{aligned}$

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.

Lemma 2.2.

Suppose that $p (x_{0}, y_{0})$ and $q (x_{0}, y_{0})$ are real numbers, and $λ_{1}$ and $λ_{2}$ are two roots of characteristic equation (9). Then

$| λ_{1} | < 1$ and $| λ_{2} | < 1$ if and only if $| p (x_{0}, y_{0}) | < 1 + q (x_{0}, y_{0})$ and $q (x_{0}, y_{0}) < 1;$
$| λ_{1} | < 1$ and $| λ_{2} | > 1$ $($ or $| λ_{1} | > 1$ and $| λ_{2} | < 1)$ if and only if $| 1 + q (x_{0}, y_{0}) | > p (x_{0}, y_{0});$
$| λ_{1} | > 1$ and $| λ_{2} | > 1$ if and only if $| p (x_{0}, y_{0}) | < 1 + q (x_{0}, y_{0})$ and $q (x_{0}, y_{0}) > 1;$
$λ_{1} = 1$ and $| λ_{2} | \neq 1$ if and only if $p (x_{0}, y_{0}) = 1 + q (x_{0}, y_{0})$ and $| q (x_{0}, y_{0}) | \neq 1;$
$λ_{1} = - 1$ and $| λ_{2} | \neq 1$ if and only if $p (x_{0}, y_{0}) = - 1 - q (x_{0}, y_{0})$ and $| q (x_{0}, y_{0}) | \neq 1;$
$λ_{1}$ and $λ_{2}$ are conjugate complex and $| λ_{1} | = | λ_{2} | = 1$ if and only if $| p (x_{0}, y_{0}) | < 2,$ and $q (x_{0}, y_{0}) = 1.$

We recall some definitions of topological types for a fixed point $(x_{0}, y_{0})$ . $(x_{0}, y_{0})$ is called a sink if both $| λ_{1} | < 1$ and $| λ_{2} | < 1$ , that is, $(x_{0}, y_{0})$ is locally asymptotic stable. $(x_{0}, y_{0})$ is called a source if both $| λ_{1} | > 1$ and $| λ_{2} | > 1$ , that is, $(x_{0}, y_{0})$ is repelling. $(x_{0}, y_{0})$ is called a saddle if $| λ_{1} | > 1$ and $| λ_{2} | < 1$ (or $| λ_{1} | < 1$ and $| λ_{2} | > 1$ ), that is, $(x_{0}, y_{0})$ is locally unstable. And $(x_{0}, y_{0})$ is called non-hyperbolic if either $| λ_{1} | = 1$ or $| λ_{2} | = 1$ .

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.

In region I, both $O (0, 0)$ and $A_{1} (x_{1}^{*}, 0)$ are saddle;
In region II, $O (0, 0)$ is a saddle, $A_{1} (x_{1}^{*}, 0)$ is a sink;
In region III, $O (0, 0)$ is a sink, both $A_{1} (x_{1}^{*}, 0)$ and $A_{2} (x_{2}^{*}, 0)$ are saddle;
In region IV, $O (0, 0)$ is a sink, $A_{1} (x_{1}^{*}, 0)$ is a sink and $A_{2} (x_{2}^{*}, 0)$ is a saddle;
In region V, $O (0, 0)$ is a sink, $A_{1} (x_{1}^{*}, 0)$ is a saddle and $A_{2} (x_{2}^{*}, 0)$ is a source;
In region VI, $O (0, 0)$ is a sink.
On curve $C_{1},$ $O (0, 0)$ is a sink, and $A (1 / 2 β, 0)$ is non-hyperbolic;
On curve $C_{2},$ $O (0, 0)$ is non-hyperbolic, and $B (1 / β, 0)$ is a sink if $β > 1,$ a saddle if $β < 1$ and non-hyperbolic if $β = 1;$
On curve $C_{3},$ $O (0, 0)$ is a sink, $A_{1} (x_{1}^{*}, 0)$ is a saddle and $A_{2} (x_{2}^{*}, 0)$ is non-hyperbolic if $β < \frac{1}{2};$ $O (0, 0)$ is a sink and $A_{1} (x_{1}^{*}, 0)$ is non-hyperbolic if $β = \frac{1}{2};$ $O (0, 0)$ is a sink, $A_{1} (x_{1}^{*}, 0)$ is non-hyperbolic and $A_{2} (x_{2}^{*}, 0)$ is a saddle if $\frac{1}{2} < β < 1;$ both $O (0, 0)$ and $A_{1} (x_{1}^{*}, 0)$ are non-hyperbolic if $β = 1,$ $O (0, 0)$ is a sink and $A_{1} (x_{1}^{*}, 0)$ is non-hyperbolic if $β > 1$ .

Proof.

Note that the eigenvalues of Jacobian matrix at $O (0, 0)$ are $λ_{1} = α, λ_{2} = 0$ . If $x \neq 0$ and y=0, the Jacobian matrix (8) takes the following form: $J G (x, 0) = (\begin{matrix} \frac{2 - α}{x + α} & - x \\ 0 & x \end{matrix}) .$ The eigenvalues of this matrix are $λ_{1} = \frac{2 - α}{x + α}, λ_{2} = x .$

For the fixed point $O (0, 0)$ , the Jacobian matrix (8) reduces into the following form: $J G (0, 0) = (\begin{matrix} α & 0 \\ 0 & 0 \end{matrix}) .$ The eigenvalues for this matrix are $λ_{1} = α, λ_{2} = 0 < 1$ . Thus, $O (0, 0)$ is a sink if $α < 1$ , a saddle if $α > 1$ and non-hyperbolic if $α = 1$ .

For the boundary point $A_{1} (x_{1}^{*}, 0) = (1 / 2 β + (1 / 2 β) \sqrt{1 - 4 β (1 - α)}, 0)$ , which exists under the condition $α \geq 1 - 1 / 4 β$ , the eigenvalues are $\begin{aligned} λ_{1} & = \frac{2 - α}{x_{1}^{*} + α} = \frac{4 β - 2 α β}{2 α β + 1 + \sqrt{1 - 4 β (1 - α)}}, \\ λ_{2} & = x_{1}^{*} = \frac{1}{2 β} + \frac{1}{2 β} \sqrt{1 - 4 β (1 - α)} . \end{aligned}$

(i $A_{1}$ ) For the eigenvalue $λ_{1}$ , if $α = 1 - 1 / 4 β$ , then $λ_{1} = 1$ , hence $A_{1} (x_{1}^{*}, 0)$ is non-hyperbolic. If $α > 1 - 1 / 4 β$ , we claim that $- 1 < λ_{1} < 1$ . In fact, on the one hand, it is easy to verify that $λ_{1} < 1$ , on the other hand, $\begin{aligned} λ_{1} & = \frac{4 β - 2 α β}{2 α β + 1 + \sqrt{1 - 4 β (1 - α)}} > - 1 \\ \Leftrightarrow \sqrt{1 - 4 β (1 - α)} > - 1 - 4 β, \end{aligned}$ which is obviously true under the condition $α > 1 - 1 / 4 β$ .

(ii $A_{1}$ ) For the eigenvalue $λ_{2} = x_{1}^{*} = 1 / 2 β + (1 / 2 β) \sqrt{1 - 4 β (1 - α)}$ ,

(aa) if $β \leq \frac{1}{2}$ , then $1 / 2 β \geq 1$ , and $λ_{2} > 1$ ;

(bb) if $β > \frac{1}{2}$ , then $0 < λ_{2} < 1$ if $α < β$ , $λ_{2} = 1$ if $α = β$ , and $λ_{2} > 1$ if $α > β$ . In fact, under the condition $β > \frac{1}{2}$ , $\begin{aligned} λ_{2} & = \frac{1}{2 β} + \frac{1}{2 β} \sqrt{1 - 4 β (1 - α)} < 1 \\ \Leftrightarrow \sqrt{1 - 4 β (1 - α)} < 2 β - 1 \\ \Leftrightarrow 4 α β < 4 β^{2} \\ \Leftrightarrow α < β . \end{aligned}$ Based on cases (i $A_{1}$ ) and (ii $A_{1}$ ), we can obtain the topological types for the fixed point $A_{1} (x_{1}^{*}, 0)$ .

For the boundary point $A_{2} (x_{2}^{*}, 0) = (1 / 2 β - (1 / 2 β) \sqrt{1 - 4 β (1 - α)}, 0)$ , which exists under the conditions $α < 1$ and $α > 1 - 1 / 4 β$ , the eigenvalues are $\begin{aligned} λ_{1} & = \frac{2 - α}{x_{2}^{*} + α} = \frac{4 β - 2 α β}{2 α β + 1 - \sqrt{1 - 4 β (1 - α)}}, \\ λ_{2} & = x_{2}^{*} = \frac{1}{2 β} - \frac{1}{2 β} \sqrt{1 - 4 β (1 - α)} . \end{aligned}$ (i $A_{2}$ ) For the eigenvalue $λ_{1}$ , it is easy to verify that $2 α β + 1 - \sqrt{1 - 4 β (1 - α)} > 0$ . We claim that $λ_{1} > 1$ under the conditions $α < 1$ and $α > 1 - 1 / 4 β$ . In fact, $\begin{aligned} λ_{1} & = \frac{4 β - 2 α β}{2 α β + 1 - \sqrt{1 - 4 β (1 - α)}} > 1 \\ \Leftrightarrow 4 β - 2 α β > 2 α β + 1 - \sqrt{1 - 4 β (1 - α)} \\ \Leftrightarrow \sqrt{1 - 4 β (1 - α)} > 1 + 4 α β - 4 β \\ \Leftrightarrow 4 (1 - α)^{2} > 1 - α \\ \Leftrightarrow α > 1 - \frac{1}{4 β} . \end{aligned}$

(ii $A_{2}$ ) For the eigenvalue $λ_{2} = x_{2}^{*} = 1 / 2 β - (1 / 2 β) \sqrt{1 - 4 β (1 - α)}$ , we first claim that $λ_{2} > - 1$ . In fact, $\begin{aligned} λ_{2} & = \frac{1}{2 β} - \frac{1}{2 β} \sqrt{1 - 4 β (1 - α)} > - 1 \\ \Leftrightarrow 2 β + 1 > \sqrt{1 - 4 β (1 - α)} \\ \Leftrightarrow 2 + β > α, \end{aligned}$ which is obviously true under the condition $α < 1$ .

We now claim that (aa) if $β \leq \frac{1}{2}$ , then $λ_{2} < 1$ if $α > β$ , $λ_{2} = 1$ if $α = β$ , and $λ_{2} > 1$ if $α < β$ . In fact, under the condition $β \leq \frac{1}{2}$ , $\begin{aligned} λ_{2} & = \frac{1}{2 β} - \frac{1}{2 β} \sqrt{1 - 4 β (1 - α)} < 1 \\ \Leftrightarrow \sqrt{1 - 4 β (1 - α)} > 1 - 2 β \\ \Leftrightarrow α > β . \end{aligned}$ (bb) if $β > \frac{1}{2}$ , then $1 - 2 β < 0$ , and we have $λ_{2} < 1$ according to the above equivalent conditions. Based on the cases (i $A_{2}$ ) and (ii $A_{2}$ ), we can obtain the topological types for the fixed point $A_{2} (x_{2}^{*}, 0)$ . Moreover, the statements (1)– (9) in this lemma follow. The proof is completed.

From Lemma 2.1, we know that system (3) has a unique positive fixed point $P (x_{0}, y_{0})$ if $α > β > 0$ , that is $(α, β) \in I \cup C_{22} \cup I I I$ . We now consider the stability of the unique positive fixed point $P (x_{0}, y_{0})$ of system (3), where $x_{0} = l e^{l} / e^{l} - 1, y_{0} = l$ .

The characteristic equation of the Jacobin matrix $J (x_{0}, y_{0})$ of system (3) about the unique positive fixed point $P (l e^{l} / e^{l} - 1, l)$ is given by (10) $λ^{2} - p λ + q = 0,$ (10) where $\begin{aligned} p & = \frac{(2 e^{l} - α) (e^{l} - 1)}{l e^{l} + α (e^{l} - 1)} + \frac{l}{e^{l} - 1}, \\ q & = \frac{(2 e^{l} - α) l}{l e^{l} + α (e^{l} - 1)} + l . \end{aligned}$ From Lemma 2.2, we can obtain the local stability and topological classification of the unique positive fixed point by computation.

Theorem 2.4.

If $α > β > 0,$ that is, $(α, β) \in I \cup C_{22} \cup I I I,$ then system (3) has a unique positive fixed point $P (l e^{l} / (e^{l} - 1), l)$ of system (3) in the interior of $R_{+}^{2},$ where $0 < l < l_{*} < 2$ and $l_{*}$ is a unique positive root of $e^{l} - 1 + 2 l - l e^{l} = 0$ . The topological classification of $P (l e^{l} / (e^{l} - 1), l)$ is as follows:

$P (l e^{l} / (e^{l} - 1), l)$ is a sink if and only if $0 < α < (l + 1) l e^{l} / (l + (1 - l) (e^{l} - 1))$ and $- 1 - q < p < 1 + q;$
$P (l e^{l} / (e^{l} - 1), l)$ is a saddle if and only if either p<1+q or $p < - 1 - q;$
$P (l e^{l} / (e^{l} - 1), l)$ is non-hyperbolic, whose eigenvalues are 1 and the other real number, if and only if $α = e^{l} (2 - 4 e^{l} + 2 e^{2 l} + 3 l - 3 e^{l} l + 2 l^{2} - e^{l} l^{2}) / (- 1 + e^{l}) (- 2 + 2 e^{l} - 3 l + e^{l} l);$
$P (l e^{l} / (e^{l} - 1), l)$ is a source if and only if $α > (l + 1) l e^{l} / (l + (1 - l) (e^{l} - 1))$ and $- 1 - q < p < 1 + q;$
$P (l e^{l} / (e^{l} - 1), l)$ is non-hyperbolic, whose eigenvalues are a pair of conjugate complex number with modulus one, if and only if $α = (l + 1) l e^{l} / (e^{l} - 1 + 2 l - l e^{l})$ .

An illustration of the results of Theorem 2.4 is given in Figure 3, in which the stability of $P (l e^{l} / (e^{l} - 1), l)$ are shown in the parameter region ${(α, β) : α > 0, β > 0}$ . Specifically, the red curve corresponds to the condition $α = (l + 1) l e^{l} / (e^{l} - 1 + 2 l - l e^{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 γ$ .

Figure 3. An illustration of the stability of $P (l e^{l} / (e^{l} - 1), l)$ .

Display full size

3. 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 $C_{1} = \{(α, β) : α = 1 - \frac{1}{4 β}, β > \frac{1}{4}\} .$

Proof.

For $β > \frac{1}{4}$ , denote $\tilde{α} = 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 $x_{n + 1} = \frac{x_{n}^{2} + α x_{n}}{1 + β x_{n}^{2}},$ and the corresponding map is $f (x) = (x^{2} + α x) / (1 + β x^{2})$ . Now at $α = \tilde{α}$ and $\tilde{x} = 1 / 2 β$ , we have $f_{x} |_{(\tilde{α}, \tilde{x})} = {[\frac{- α β x^{2} + 2 x + α}{(1 + β x^{2})^{2}}]|}_{(\tilde{α}, \tilde{x})} = 1.$ This is the non-hyperbolic condition for fold bifurcation. Now, we will check the non-degenerate condition for the map, that is, $f_{x x} |_{(\tilde{α}, \tilde{x})} = - 1 / (1 + 4 β) \neq 0$ and $f_{α} |_{(\tilde{α}, \tilde{x})} = x / (1 + β x^{2}) |_{(\tilde{α}, \tilde{x})} = 2 / (4 β + 1) \neq 0$ . 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}$ .

On the other hand, from the conclusion (v) of Theorem 10, we know that the unique positive fixed point $P (l e^{l} / (e^{l} - 1), l)$ of system (3) is non-hyperbolic with a pair of conjugate complex eigenvalues with modulus one if $α = (l + 1) l e^{l} / (l + (1 - l) (e^{l} - 1))$ . By Equation (7), we obtain $β = \frac{e^{- 2 l} (e^{l} - 1) (e^{l} (2 - l) + e^{2 l} (l - 1) + l^{2} - 1)}{l^{2} (e^{l} - 1 + 2 l - l e^{l})} .$

For simplicity, we denote the parameters satisfying $α = (l + 1) l e^{l} / l + (1 - l) (e^{l} - 1)$ and $β = e^{- 2 l} (- 1 + e^{l}) (- 1 - e^{l} (- 2 + l) + e^{2 l} (- 1 + l) + l^{2}) / l^{2} (e^{l} - 1 + 2 l - l e^{l})$ as $\begin{aligned} H_{P} & = \{(α, β) : α = \frac{(l + 1) l e^{l}}{l + (1 - l) (e^{l} - 1)}, \\ β = \frac{e^{- 2 l} (e^{l} - 1) (e^{l} (2 - l) + e^{2 l} (l - 1) + l^{2} - 1)}{l^{2} (e^{l} - 1 + 2 l - l e^{l})}, 0 < l < l_{*}\} . \end{aligned}$

A numerical illustration of the Neimark–Sacker bifurcation curve $H_{P}$ is given in Figure 4, in which the red curve denotes $H_{P}$ on the parameter plane $(α, β)$ .

Figure 4. A numerical illustration of the curve $H_{P}$ , and the stability of $P (l e^{l} / (e^{l} - 1), l)$ .

Display full size

We now consider the Neimark–Sacker bifurcation of system (3) at the unique positive fixed point $P (l e^{l} / e^{l} - 1, l)$ . Choose the parameters $(α_{1}, β_{1}) \in H_{P}$ , then system (3) for $(α, β) = (α_{1}, β_{1})$ becomes (11) $x_{n + 1} = \frac{x_{n}^{2} + α_{1} x_{n}}{1 + β_{1} x_{n}^{2}} e^{- y_{n}}, y_{n + 1} = x_{n} (1 - e^{- y_{n}}) .$ (11) It is clear that $P (x^{*}, y^{*})$ is the unique positive fixed point of system (11), where $x^{*} = l_{1} e^{l_{1}} / e^{l_{1}} - 1, y^{*} = l_{1}$ and $α_{1} = (l_{1} + 1) l_{1} e^{l_{1}} / l_{1} + (1 - l_{1}) (e^{l_{1}} - 1)$ .

Consider the small perturbation of $α_{1}$ by introducing a parameter $α^{*}$ as follows: (12) $x_{n + 1} = \frac{x_{n}^{2} + (α_{1} + α^{*}) x_{n}}{1 + β_{1} x_{n}^{2}} e^{- y_{n}}, y_{n + 1} = x_{n} (1 - e^{- y_{n}}) .$ (12) where $| α^{*} | ≪ 1$ , that is, it is a sufficiently small parameter.

The characteristic equation of the Jacobian matrix of linearized system of (12) about the unique positive fixed point $P (x^{*}, y^{*})$ is given by $λ^{2} - p (α^{*}) λ + q (α^{*}) = 0,$ where $\begin{aligned} p (α^{*}) & = \frac{(2 e^{l_{1}} - (α_{1} + α^{*})) (e^{l_{1}} - 1)}{l_{1} e^{l_{1}} + (α_{1} + α^{*}) (e^{l_{1}} - 1)} + \frac{l_{1}}{e^{l_{1}} - 1}, \\ q (α^{*}) & = \frac{(2 e^{l_{1}} - (α_{1} + α^{*})) l_{1}}{l_{1} e^{l_{1}} + (α_{1} + α^{*}) (e^{l_{1}} - 1)} + l_{1} . \end{aligned}$ Moreover when $α^{*}$ varies in a small neighbourhood of 0, the roots of the characteristic equation are $\begin{aligned} λ_{1, 2} & = \frac{- p (α^{*}) \mp ι \sqrt{4 q (α^{*}) - p^{2} (α^{*})}}{2}, \\ = - \frac{1}{2} \frac{(2 e^{l_{1}} - (α_{1} + α^{*})) (e^{l_{1}} - 1)}{l_{1} e^{l_{1}} + (α_{1} + α^{*}) (e^{l_{1}} - 1)} - \frac{1}{2} \frac{l_{1}}{e^{l_{1}} - 1} \\ \mp \frac{ι}{2} \sqrt{4 l_{1} - {(\frac{(2 e^{l_{1}} - (α_{1} + α^{*})) (e^{l_{1}} - 1)}{l_{1} e^{l_{1}} + (α_{1} + α^{*}) (e^{l_{1}} - 1)} - \frac{l}{e^{l_{1}} - 1})}^{2}}, \end{aligned}$ then we have $\begin{aligned} | λ_{1, 2} | = (q (α^{*}))^{1 / 2}, \frac{d | λ_{1, 2} |}{d α^{*}} |_{α^{*} = 0} = - \frac{(l_{1} + (1 - l_{1}) (e^{l_{1}} - 1))^{2}}{2 \sqrt{q (α^{*})} l_{1} e^{l_{1}} (l_{1} + 2 (e^{l_{1}} - 1))} < 0. \end{aligned}$ Moreover, it is required that when $α^{*} = 0$ , $λ_{1, 2}^{k} \neq 1$ , k=1,2,3,4, which is equivalent to $p (0) \neq - 2, 0, 1, 2$ . Since $p (0)^{2} - 4 q (0) < 0$ and $q (0) = 1$ , therefore $p (0)^{2} < 4$ and hence $p (0) \neq \pm 2$ . So we only need to require that $p (0) \neq 0, 1$ , that is, (13) $(e^{l_{1}} - 1)^{2} \neq \frac{l_{1}^{2}}{l_{1} - 1}, \frac{l_{1}^{2} - l_{1} (e^{l_{1}} - 1)}{l_{1} - 1} .$ (13)

Let $u_{n} = x_{n} - x^{*}, v_{n} = y_{n} - y^{*}$ then we transform the fixed point $P (x^{*}, y^{*})$ of system (12) into the origin. By calculating we obtain (14) $\begin{aligned} u_{n + 1} & = \frac{(u_{n} + x^{*})^{2} + (α_{1} + α^{*}) (u_{n} + x^{*})}{1 + β_{1} (u_{n} + x^{*})^{2}} e^{- (v_{n} + y^{*})} - x^{*}, \\ v_{n + 1} & = (u_{n} + x^{*}) (1 - e^{- (v_{n} + y^{*})}) - y^{*} . \end{aligned}$ (14) In the following, we study the normal form of system (14) when $α^{*} = 0$ . Expanding (14) as a Taylor series at $(u_{n}, v_{n}) = (0, 0)$ to the third order, we obtain (15) $\begin{aligned} u_{n + 1} = & a_{11} u_{n} + a_{12} v_{n} + a_{13} u_{n}^{2} + a_{14} u_{n} v_{n} + a_{15} v_{n}^{2} + a_{16} u_{n}^{3} + a_{17} u_{n}^{2} v_{n} \\ + a_{18} u_{n} v_{n}^{2} + a_{19} v_{n}^{3} + o ((| u_{n} | + | v_{n} |)^{3}), \\ v_{n + 1} = & a_{21} u_{n} + a_{22} v_{n} + a_{23} u_{n} v_{n} + a_{24} v_{n}^{2} + a_{25} u_{n} v_{n}^{2} + a_{26} v_{n}^{3} + o ((| u_{n} | + | v_{n} |)^{3}), \end{aligned}$ (15) where $\begin{aligned} a_{11} & = \frac{2 e^{y^{*}} - α_{1}}{x^{*} + α_{1}}, a_{12} = - x^{*}, \\ a_{13} & = \frac{1}{(1 + β_{1} x^{*^{2}})^{3}} [(1 - 3 α_{1} β_{1} x^{*} - 3 β_{1} x^{*^{2}} + α_{1} β_{1}^{2} x^{*^{3}}) e^{- y^{*}}], \\ a_{14} & = - \frac{2 e^{y^{*}} - α_{1}}{x^{*} + α_{1}}, a_{15} = \frac{1}{2} x^{*}, \\ a_{16} & = - \frac{1}{(1 + β_{1} x^{*^{2}})^{4}} [(α_{1} β_{1} + 4 β_{1} x^{*} - 6 α_{1} β_{1}^{2} x^{*^{2}} - 4 β_{1}^{2} x^{*^{3}} + α_{1} β_{1} x^{*^{4}}) e^{- y^{*}}], \\ a_{17} & = - \frac{1}{(1 + β_{1} x^{*^{2}})^{3}} [(1 - 3 α_{1} β_{1} x^{*} - 3 β_{1} x^{*^{2}} + α_{1} β_{1}^{2} x^{*^{3}}) e^{- y^{*}}], \\ a_{18} & = \frac{1}{2} \frac{2 e^{y^{*}} - α_{1}}{x^{*} + α_{1}}, a_{19} = - \frac{1}{6} x^{*}, \\ a_{21} & = 1 - e^{- y^{*}}, a_{22} = x^{*} e^{- y^{*}}, a_{23} = e^{- y^{*}}, a_{24} = - \frac{1}{2} x^{*} e^{- y^{*}}, \\ a_{25} & = - \frac{1}{2} e^{- y^{*}}, a_{26} = \frac{1}{6} x^{*} e^{- y^{*}} . \end{aligned}$

Now, let $\begin{aligned} η & = - \frac{1}{2} \frac{(2 e^{l_{1}} - α_{1}) (e^{l_{1}} - 1)}{l_{1} e^{l_{1}} + α_{1} (e^{l_{1}} - 1)} - \frac{1}{2} \frac{l_{1}}{e^{l_{1}} - 1}, \\ ζ & = \frac{1}{2} \sqrt{4 l_{1} - {(\frac{(2 e^{l_{1}} - α_{1}) (e^{l_{1}} - 1)}{l_{1} e^{l_{1}} + α_{1} (e^{l_{1}} - 1)} - \frac{l}{e^{l_{1}} - 1})}^{2}}, \end{aligned}$ and $T = (\begin{matrix} a_{12} & 0 \\ η - a_{11} & - ζ \end{matrix}),$ then T is invertible. Using translation $(\begin{matrix} u_{n} \\ v_{n} \end{matrix}) = (\begin{matrix} a_{12} & 0 \\ η - a_{11} & - ζ \end{matrix}) (\begin{matrix} X_{n} \\ Y_{n} \end{matrix}),$ then system (15) becomes of the following form: (16) $\begin{aligned} X_{n + 1} & = η X_{n} - ζ Y_{n} + \bar{F} (X_{n}, Y_{n}), \\ Y_{n + 1} & = ζ X_{n} + η Y_{n} + \bar{G} (X_{n}, Y_{n}), \end{aligned}$ (16) where (17) $\begin{aligned} \bar{F} (X_{n}, Y_{n}) & = c_{11} X_{n}^{2} + c_{12} X_{n} Y_{n} + c_{13} Y_{n}^{2} + c_{14} X_{n}^{3} + c_{15} X_{n}^{2} Y_{n} + c_{16} X_{n} Y_{n}^{2} \\ + c_{17} Y_{n}^{3} + o ((| X_{n} | + | Y_{n} |)^{3}), \\ \bar{G} (X_{n}, Y_{n}) & = c_{21} X_{n}^{2} + c_{22} X_{n} Y_{n} + c_{23} Y_{n}^{2} + c_{24} X_{n}^{3} + c_{25} X_{n}^{2} Y_{n} + c_{26} X_{n} Y_{n}^{2} \\ + c_{27} Y_{n}^{3} + o ((| X_{n} | + | Y_{n} |)^{3}), \end{aligned}$ (17) and $\begin{aligned} c_{11} & = a_{12} a_{13} + (η - a_{11}) a_{14} + \frac{(η - a_{11})^{2} a_{15}}{a_{12}}, c_{12} = - ζ [a_{14} + \frac{2 (η - a_{11}) a_{15}}{a_{12}}], \\ c_{13} & = \frac{ζ^{2} a_{15}}{a_{12}}, c_{14} = a_{12}^{2} a_{16} + (η - a_{11}) a_{12} a_{17} + (η - a_{11})^{2} a_{18} + \frac{(η - a_{11})^{3} a_{19}}{a_{12}}, \\ c_{15} & = - ζ [a_{12} a_{17} + 2 (η - a_{11}) a_{18} + \frac{3 (η - a_{11})^{2} a_{19}}{a_{12}}], \\ c_{16} & = ζ^{2} [a_{18} + \frac{3 (η - a_{11}) a_{19}}{a_{12}}], c_{17} = - \frac{ζ^{3} a_{19}}{a_{12}}, \\ c_{21} & = \frac{η - a_{11}}{ζ} [a_{12} a_{13} + (η - a_{11}) a_{14} - a_{12} a_{23} + \frac{1}{a_{12}} ((η - a_{11}) a_{15} - a_{12} a_{24}) (η - a_{11})], \\ c_{22} & = a_{12} a_{23} - (η - a_{11}) a_{14} - \frac{2}{a_{12}} ((η - a_{11}) a_{15} - a_{12} a_{24}) (η - a_{11}), \end{aligned}$ $\begin{aligned} c_{23} & = \frac{ζ}{a_{12}} [(η - a_{11}) a_{15} - a_{12} a_{24}], \\ c_{24} & = \frac{η - a_{11}}{ζ} [a_{12}^{2} a_{16} + a_{12} a_{17} (η - a_{11}) + ((η - a_{11}) a_{18} - a_{12} a_{25}) (η - a_{11}) \\ + \frac{1}{a_{12}} ((η - a_{11}) a_{19} - a_{12} a_{26}) (η - a_{11})^{2}], \\ c_{25} & = (a_{11} - η) [a_{12} a_{17} + 2 ((η - a_{11}) a_{18} - a_{12} a_{25}) \\ + \frac{3}{a_{12}} ((η - a_{11}) a_{19} - a_{12} a_{26}) (η - a_{11})], \\ c_{26} & = ζ [(η - a_{11}) a_{18} - a_{12} a_{25} + \frac{3}{a_{12}} ((η - a_{11}) a_{19} - a_{12} a_{26}) (η - a_{11})], \\ c_{27} & = - \frac{ζ^{2}}{a_{12}} [(η - a_{11}) a_{19} - a_{12} a_{26}] . \end{aligned}$ Furthermore, $\begin{aligned} {\bar{F}}_{X_{n} X_{n}} |_{(0, 0)} & = 2 c_{11}, {\bar{F}}_{X_{n} Y_{n}} |_{(0, 0)} = c_{12}, {\bar{F}}_{Y_{n} Y_{n}} |_{(0, 0)} = 2 c_{13}, {\bar{F}}_{X_{n} X_{n} X_{n}} |_{(0, 0)} = 6 c_{14}, \\ {\bar{F}}_{X_{n} X_{n} Y_{n}} |_{(0, 0)} & = 2 c_{15}, {\bar{F}}_{X_{n} Y_{n} Y_{n}} |_{(0, 0)} = 2 c_{16}, {\bar{F}}_{Y_{n} Y_{n} Y_{n}} |_{(0, 0)} = 6 c_{17}, \end{aligned}$ and $\begin{aligned} {\bar{G}}_{X_{n} X_{n}} |_{(0, 0)} & = 2 c_{21}, {\bar{G}}_{X_{n} Y_{n}} |_{(0, 0)} = c_{22}, {\bar{G}}_{Y_{n} Y_{n}} |_{(0, 0)} = 2 c_{23}, {\bar{G}}_{X_{n} X_{n} X_{n}} |_{(0, 0)} = 6 c_{24}, \\ {\bar{G}}_{X_{n} X_{n} Y_{n}} |_{(0, 0)} & = 2 c_{25}, {\bar{G}}_{X_{n} Y_{n} Y_{n}} |_{(0, 0)} = 2 c_{26}, {\bar{G}}_{Y_{n} Y_{n} Y_{n}} |_{(0, 0)} = 6 c_{27} . \end{aligned}$

In order to guarantee the Neimark–Sacker bifurcation for (16), we require that the following discriminatory quantity is not zero (see [12, 13]) (18) $Ω = - Re [\frac{(1 - 2 \bar{λ}) {\bar{λ}}^{2}}{1 - λ} ξ_{11} ξ_{20}] - \frac{1}{2} ∥ ξ_{11} ∥^{2} - ∥ ξ_{02} ∥^{2} + Re (\bar{λ} ξ_{21}),$ (18) where (19) $\begin{aligned} ξ_{02} & = \frac{1}{8} [{\bar{F}}_{X_{n} X_{n}} - {\bar{F}}_{Y_{n} Y_{n}} + 2 {\bar{G}}_{X_{n} Y_{n}} + ι ({\bar{G}}_{X_{n} X_{n}} - {\bar{G}}_{Y_{n} Y_{n}} + 2 {\bar{F}}_{X_{n} Y_{n}})] |_{(0, 0)}, \\ ξ_{11} & = \frac{1}{4} [{\bar{F}}_{X_{n} X_{n}} + {\bar{F}}_{Y_{n} Y_{n}} + ι ({\bar{G}}_{X_{n} X_{n}} + {\bar{G}}_{Y_{n} Y_{n}})] |_{(0, 0)}, \\ ξ_{20} & = \frac{1}{8} [{\bar{F}}_{X_{n} X_{n}} - {\bar{F}}_{Y_{n} Y_{n}} + 2 {\bar{G}}_{X_{n} Y_{n}} + ι ({\bar{G}}_{X_{n} X_{n}} - {\bar{G}}_{Y_{n} Y_{n}} - 2 {\bar{F}}_{X_{n} Y_{n}})] |_{(0, 0)}, \\ ξ_{21} & = \frac{1}{16} [{\bar{F}}_{X_{n} X_{n} X_{n}} + {\bar{F}}_{X_{n} Y_{n} Y_{n}} + {\bar{G}}_{X_{n} X_{n} Y_{n}} + {\bar{G}}_{Y_{n} Y_{n} Y_{n}} \\ + ι ({\bar{G}}_{X_{n} X_{n} X_{n}} + {\bar{G}}_{X_{n} Y_{n} Y_{n}} - {\bar{F}}_{X_{n} X_{n} Y_{n}} - {\bar{F}}_{Y_{n} Y_{n} Y_{n}})] |_{(0, 0)} . \end{aligned}$ (19)

After calculating, we obtain (20) $\begin{aligned} ξ_{02} & = \frac{1}{4} [c_{11} - c_{13} + c_{22} + ι (c_{21} - c_{23} + c_{12})], \\ ξ_{11} & = \frac{1}{2} [c_{11} + c_{13} + ι (c_{21} + c_{23})], \\ ξ_{20} & = \frac{1}{4} [c_{11} - c_{13} + c_{22} + ι (c_{21} - c_{23} - c_{12})], \\ ξ_{21} & = \frac{1}{8} [3 c_{14} + c_{16} + c_{25} + 3 c_{27} + ι (3 c_{24} + c_{26} - c_{15} - 3 c_{17})], \end{aligned}$ (20)

Based on the above analysis and the Neimark–Sacker bifurcation theorem discussed in [12], we finally arrive at the following theorem.

Theorem 3.2.

If the condition (13) holds, and $Ω \neq 0,$ then the system (11) undergoes a Neimark–Sacker bifurcation at the unique positive fixed point $P (x^{*}, y^{*})$ as the parameters $(α, β)$ go through the curve $H_{P}$ . Moreover, if $Ω < 0$ (respectively $Ω > 0$ ), then an attracting (respectively repelling) invariant closed curve bifurcates from the fixed point $P (x^{*}, y^{*})$ .

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.

4. 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) is called eventually bounded if there exists a natural number N, and real numbers m and M such that $m \leq x_{n} \leq M$ and $m \leq y_{n} \leq M$ for all $n \geq N$ . A solution ${x_{n}, y_{n}}$ of system (3) is called eventually fixed point if there exists a fixed point $(x^{*}, y^{*})$ of system (3) and a natural number N such that the solution ${x_{N}, y_{N}} = (x^{*}, y^{*})$ . Similarly, A solution ${x_{n}, y_{n}}$ of system (3) is called eventually periodic point with period k if there exists a period point $(x^{*}, y^{*})$ with period k of system (3) and a natural number N such that the solution ${x_{N}, y_{N}} = (x^{*}, y^{*})$ .

Lemma 4.1.

All solutions ${(x_{n}, y_{n})}$ of system (3) in the first quadrant $R_{+}^{2}$ are eventually bounded.

Proof.

Let ${(x_{n}, y_{n})}$ be any positive solution of system (3). It is clear that the first quadrant $R_{+}^{2}$ is positive invariant for system (3). Moreover, any points ${(0, y_{n})}$ in the positive y-axis are eventually fixed point $O (0, 0)$ , that is, it becomes the fixed point $O (0, 0)$ after one step, and the positive x-axis is also positive invariant for system (3). We divide three cases to discuss solutions ${(x_{n}, 0)}$ of system (3) in the positive x-axis.

Case (i): $0 < α < 1 - 1 / 4 β$ and $β > \frac{1}{4}$ . In this case, $0 < x_{n + 1} < x_{n}$ . For any $x_{0} > 0$ , we have $lim_{n \to + \infty} x_{n} = 0.$ This implies that all solutions ${(x_{n}, 0)}$ of system (3) with $x_{0} > 0$ are decreasing to the fixed point $O (0, 0)$ .

Case (ii): $α = 1 - 1 / 4 β$ and $β > \frac{1}{4}$ . In this case, system (3) has two fixed points $O (0, 0)$ and $A (1 / 2 β, 0)$ . All solutions ${(x_{n}, 0)}$ of system (3) with $0 < x_{0} < 1 / 2 β$ are decreasing to the fixed point $O (0, 0)$ , and all solutions ${(x_{n}, 0)}$ of system (3) with $x_{0} > 1 / 2 β$ are decreasing to the fixed point $A (1 / 2 β, 0)$ .

Case (iii): $α > 1 - 1 / 4 β$ and $β > \frac{1}{4}$ . In this case, system (3) has three fixed points $O (0, 0)$ , $A_{1} (x_{1}^{*}, 0)$ and $A_{2} (x_{2}^{*}, 0)$ . It can be checked that all solutions ${(x_{n}, 0)}$ of system (3) with $0 < x_{0} < x_{2}^{*}$ ( $x_{0} > x_{1}^{*}$ ) are decreasing to the fixed point $O (0, 0)$ ( $A_{1} (x_{1}^{*}, 0)$ , resp.), and all solutions ${(x_{n}, 0)}$ of system (3) with $x_{2}^{*} < x_{0} < x_{1}^{*}$ are increasing to the fixed point $A_{1} (x_{1}^{*}, 0)$ .

Therefore, all solutions ${(x_{n}, 0)}$ of system (3) are eventually bounded by the interval $[0, x_{1}^{*} + ϵ]$ in the positive x-axis, where $0 < ϵ ≪ 1$ .

We now consider solutions ${(x_{n}, y_{n})}$ of system (3) with $x_{0} > 0$ and $y_{0} > 0$ . From system (3), we have (21) $0 < x_{n + 1} < \frac{x_{n}^{2} + α x_{n}}{1 + β x_{n}^{2}}, 0 < y_{n + 1} < x_{n}$ (21) for all nonnegative integer n.

Let $f (x) = (x^{2} + α x) / (1 + β x^{2})$ . Then $f (0) = 0$ and $lim_{n \to + \infty} f (x) = 1 / β$ . And $f (x)$ has a unique extreme point at $x^{⋄} = (1 + \sqrt{1 + α^{2} β}) / α β$ for all $x \in (0 + \infty)$ . A straightforward calculation shows that $f (x^{⋄}) = \frac{1 + \sqrt{1 + α^{2} β}}{2 β} = {Max}_{0 < x < + \infty} f (x) .$ Hence, all solutions ${(x_{n}, y_{n})}$ of system (3) are eventually bounded by the rectangle $R$ , $R = \{(x, y) : 0 \leq x \leq \frac{1 + \sqrt{1 + α^{2} β}}{2 β}, 0 \leq y \leq \frac{1 + \sqrt{1 + α^{2} β}}{2 β}\} .$ Thus, every solution of system (3) in the first quadrant $R_{+}^{2}$ is eventually bounded. The proof is completed.

Now, we investigate the global attractor for system (3) in some parameter regions, and arrive at the following theorem.

Theorem 4.2.

The fixed points and global attractors of system (3) are stated as follows:

if $0 < α < 1 - 1 / 4 β$ and $β > \frac{1}{4},$ that is $(α, β) \in V I,$ then system (3) has a unique fixed point $O (0, 0)$ in the first quadrant $R_{+}^{2},$ which is a global attractor of system (3) in $R_{+}^{2}$ ;
if $α = 1 - 1 / 4 β$ and $β > \frac{1}{4},$ that is $(α, β) \in C_{1},$ then system (3) has two fixed points $O (0, 0)$ and $A (1 / 2 β, 0)$ in the first quadrant $R_{+}^{2}$ , and the fixed point $O (0, 0)$ attracts all positive solutions of system (3) with $0 < x_{n} < 1 / 2 β$ for some n, and all the other solutions are attracted to the fixed point $A (1 / 2 β, 0)$ ;
if $(α, β) \in I V,$ then system (3) has three fixed points $O (0, 0),$ $A_{1} (x_{1}^{*}, 0)$ and $A_{2} (x_{2}^{*}, 0),$ and all solutions of system (3) with $0 < x_{n} < x_{2}^{*}$ for some n will be attracted to $O (0, 0)$ in the first quadrant $R_{+}^{2},$ and all other solutions are attracted to $A_{1} (x_{1}^{*}, 0)$ .

Proof.

(i) According to the conclusion (i) in Lemma 2.1, system (3) has a unique fixed point $O (0, 0)$ in the first quadrant $R_{+}^{2}$ , and $O (0, 0)$ is a sink by Theorem 2.3.

From the proof of lemma 4.1, we find that if $0 < α < 1 - 1 / 4 β$ and $β > \frac{1}{4}$ , then every solution ${(x_{n}, y_{n})}$ of system (3) in the first quadrant $R_{+}^{2}$ satisfies $0 < x_{n + 1} \leq x_{n} (\frac{x_{n} + α}{1 + β x_{n}^{2}}) = x_{n} (1 - \frac{β (x_{n} - \frac{1}{2 β})^{2} + 1 - α - \frac{1}{4 β}}{1 + β x_{n}^{2}}) < x_{n}, 0 \leq y_{n + 1} \leq x_{n},$ which leads that $lim_{n \to + \infty} x_{n} = 0$ and $lim_{n \to + \infty} 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 $0 < x_{n + 1} \leq x_{n} (\frac{x_{n} + α}{1 + β x_{n}^{2}}) = x_{n} (1 - \frac{β (x_{n} - \frac{1}{2 β})^{2}}{1 + β x_{n}^{2}}) < x_{n}, 0 \leq y_{n + 1} \leq x_{n} .$ Thus, we have $lim_{n \to + \infty} x_{n} = 0$ for all solutions of system (3) with $0 < x_{n} < 1 / 2 β$ for some n, and $lim_{n \to + \infty} 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 \to + \infty} y_{n} = 0$ .

(iii) If $(α, β) \in I V$ , 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 $0 < x_{n + 1} \leq x_{n} (\frac{x_{n} + α}{1 + β x_{n}^{2}}) = x_{n} (1 - β \frac{(x_{n} - x_{1}^{*}) (x_{n} - x_{2}^{*})}{1 + β x_{n}^{2}}), 0 \leq y_{n + 1} \leq x_{n} .$ Hence, $lim_{n \to + \infty} x_{n} = 0$ for all solutions of system (3) with $0 < x_{n} < x_{2}^{*}$ for some n, and $lim_{n \to + \infty} x_{n} = x_{1}^{*}$ for all other solutions. By similar arguments, we have $lim_{n \to + \infty} 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.

Figure 5. An illustration of the global attractors of system (3) in the parameter regions V I, $C_{1}$ and IV.

Display full size

5. 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 $(α, β) \in V I$ , system (3) possesses the unique equilibrium point $O (0, 0)$ , and it is a global attractor in $R_{+}^{2}$ . Figure 6(a) shows the phase portraits of system (3) for $(α, β) \in V I$ , where the phase orbits of three different initial points are given and all the orbits are eventually attracted to the equilibrium point $O (0, 0)$ . From the viewpoint of biology, if the parameters $(α, β) \in V I$ , then the populations of both host and parasitoid will go extinction with long time. When the parameters $(α, β) \in C_{1}$ , system (3) possesses two equilibrium points $O (0, 0)$ and $A (1 / 2 β, 0)$ , with $O (0, 0)$ a sink and A non-hyperbolic. Figure 6(b) shows an example of the phase portraits of system (3) for $(α, β) \in C_{1}$ , and the simulation coincides with the conclusions of Theorem 4.2.

Figure 6. Phase portraits of system (3): (a) $(α, β) = (0.5, 1.2) \in V I$ , with three initial points $(1.2, 9), (14, 9)$ , $(15, 0.2)$ and (b) $(α, β) = (0.9, 2.5) \in C_{1}$ , with three initial points $(2.2, 2), (4, 0.9)$ , $(2.5, 0.2)$ .

Display full size

If $(α, β) \in 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 $(α, β) \in 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.

Figure 7. Phase portraits of system (3). (a) $(α, β) = (0.15, 0.2) \in V$ , with four initial points $(3.2, 0.02)$ , $(1.25, 0.01)$ , $(2.5, 0.2)$ , $(0.99, 0.02),$ (b) $(α, β) = (0.85, 1) \in I V$ , with four initial points $M_{1} (0.05, 1.2)$ , $M_{2} (0.75, 1)$ , $M_{3} (2.5, 0.2)$ , $M_{4} (3.2, 0.02),$ (c) $(α, β) = (1.2, 1.8) \in I I$ , with initial points $M_{1} - M_{4}$ , (d) $(α, β) = (0.7, 0.5) \in I I I$ , with initial points $M_{1} - M_{4}$ , (e) $(α, β) = (0.9, 0.7) \in I I I$ , with initial points $M_{1} - M_{4}$ and (f) $(α, β) = (0.15, 0.1) \in I I I$ , with initial points $M_{1} - M_{4}$ .

Display full size

If $(α, β) \in I V$ , 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 $(α, β) \in I V$ , 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 $(α, β) \in I I$ , 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 $(α, β) \in I I$ 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 $(α, β) \in I \cup I I I \cup C_{21}$ , system (3) possesses a unique positive equilibrium point $P (l e^{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 $(α, β) \in I I I$ , 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 $(α, β) \in I I I$ , where $O (0, 0)$ is a sink while P is a source.

If $(α, β) \in 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 $(α, β) \in I$ . Choose $β = 0.65$ , then we can obtain the bifurcation value $\hat{α} = 5.71755$ of α according to the Neimark–Sacker condition $H_{P}$ . If $α < \hat{α}$ , as in Figure 8(c, d), P is a stable focus. If $α > \hat{α}$ , 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.

Figure 8. Phase portraits of system (3) with $(α, β) \in I$ , and $β = 0.65$ for (c) to (l). (a) $α = 1.15, β = 0.8$ , and initial points $M_{1} - M_{4}$ , (b) $α = 4, β = 0.8$ , and initial points $M_{1} - M_{4}$ , (c) $α = 5$ and initial point $(3.5, 0.3)$ , (d) $α = 5.7$ and initial point $(3.5, 0.3)$ , (e) $α = 5.73$ and initial point $(3.5, 0.3)$ , (f) $α = 5.799$ and initial point $(3.5, 0.3)$ , (g) $α = 6.2$ and initial point $(3, 0.3)$ , (h) $α = 7.5$ and initial point $(3, 0.3)$ , (i) $α = 7.5$ and initial point $(2, 1.4)$ , (j) $α = 7.5$ and initial point $(3, 0.3)$ , and the first 10,000 points are neglected. (k) $α = 13$ and initial point $(3, 0.3)$ (l) $α = 28$ and initial point $(3, 0.3)$ .

Display full size

Figure 9. Phase portraits of system (3) with initial points $M_{1} - M_{4}$ . (a) $(α, β) = (0.95, 0.1) \in I I I$ , (b) $(α, β) = (0.95, 0.9) \in I I I$ , (c) $(α, β) = (0.5, 0.5) \in C_{1}$ , (d) $(α, β) = (1, 1) \in C_{2}$ , (e) $(α, β) = (2, 2) \in C_{3}$ and (f) $(α, β) = (1, 0.5) \in C_{2}$ .

Display full size

6. 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 $(α, β) \in V & V I$ , then the populations of both host and parasitoid will go extinction with long time; (ii) if $(α, β) \in I V$ , 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 $(α, β) \in I I$ , then the population of parasitoid will go extinction while the host will tend to a steady state with long time; (iv) if $(α, β) \in I & I I I$ , 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.

Table 1. Number of fixed points in different parameter regions.

Number of fixed points	Parametric regions	Fixed points
1	$(V I)$ : $0 < α < 1 - \frac{1}{4 β}$ and $β > \frac{1}{4}$	O
2	$C_{1}$ : $α = 1 - \frac{1}{4 β}$ and $β > \frac{1}{4}$	$O, A$
2	$C_{21}$ : $α = 1$ with $β > 1$	$O, B$
2	$(I I)$ : $1 < α < β$	$O, A_{1}$
3	$(I V)$ : $α < β < \frac{1}{4 (1 - α)}$ and $\frac{1}{2} < α < 1$	$O, A_{1}, A_{2}$
3	$(V)$ : $α < β < \frac{1}{4 (1 - α)}$ and $0 < α < \frac{1}{2}$	$O, A_{1}, A_{2}$
3	$(C_{22})$ : $α = 1$ and $0 < β < 1$	$O, B, P$
3	$(I)$ : $α > 1$ and $α > β > 0$	$O, A_{1}, P$
4	$(I I I)$ : $0 < β < α < 1$	$O, A_{1}, A_{2}, P$

Acknowledgements

The authors are grateful for the anonymous referees for their valuable comments and suggestions which led to the improvement of this paper.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

Abdul Qadeer Khan's research is supported by the Higher Education Commission of Pakistan, Jiying Ma's research is supported by the National Natural Science Foundations of China [grant number 11501364] and Dongmei Xiao's research is supported by the National Natural Science Foundations of China [grant numbers 11371248 and 11431008] and the RFDP of Higher Education of China [grant number 20130073110074].

References

R. Asheghi, Bifurcation and dynamics of a discrete predator–prey system, J. Biol. Dyn. 8 (2014), pp. 161–186. doi: 10.1080/17513758.2014.927596 [Taylor & Francis Online], [Web of Science ®], [Google Scholar]
J.R. Beddington, C.A. Free and J.H. Lawton, Dynamic complexity in predator-prey models framed in difference equations, Nature 225 (1975), pp. 58–60. doi: 10.1038/255058a0 [Crossref], [Web of Science ®], [Google Scholar]
A. Bompard, I. Amat, X. Fauvergue and T. Spataro, Host–parasitoid dynamics and the success of biological control when parasitoids are prone to Allee effects, PLoS ONE 8 (2013), pp. 76768. doi: 10.1371/journal.pone.0076768 [Crossref], [PubMed], [Web of Science ®], [Google Scholar]
F. Courchamp, L. Berec and J. Gascoigne, Allee Effects in Ecology and Conservations, Oxford University Press, Oxford, 2008. [Crossref], [Google Scholar]
A. Deredec and F. Courchamp, Combined impacts of Allee effects and parasitism, Oikos 112 (2006), pp. 667–679. doi: 10.1111/j.0030-1299.2006.14243.x [Crossref], [Web of Science ®], [Google Scholar]
S. Elaydi, An Introduction to Difference Equations, 3rd ed., Springer, New York, 2005. [Google Scholar]
S. Elaydi and R.J. Sacker, Population models with Allee effect: A new model, J. Biol. Dyn. 4 (2010), pp. 397–408. doi: 10.1080/17513750903377434 [Taylor & Francis Online], [Google Scholar]
S. Jang, Discrete-time host–parasitoid models with Allee effects: Density dependence versus parasitism, J. Diff. Equ. Appl. 17 (2011), pp. 525–539. doi: 10.1080/10236190903146920 [Taylor & Francis Online], [Web of Science ®], [Google Scholar]
S. Jang and S. Diamond, A host–parasitoid interaction with Allee effects on the host, Comput. Math. Appl. 53 (2007), pp. 89–103. doi: 10.1016/j.camwa.2006.12.013 [Crossref], [Web of Science ®], [Google Scholar]
S. Kapçak, Ü. Ufuktepe and S. Elaydi, Stability and invariant manifolds of a generalized Beddington host–parasitoid model, J. Biol. Dyn. 7 (2013), pp. 233–253. doi: 10.1080/17513758.2013.849764 [Taylor & Francis Online], [Web of Science ®], [Google Scholar]
A. Kent, C. Doncaster and T. Sluckin, Consequences for predators of rescue and Allee effects on prey, Ecol. Mod. 162 (2003), pp. 233–245. doi: 10.1016/S0304-3800(02)00343-5 [Crossref], [Web of Science ®], [Google Scholar]
Y.A. Kuznetson, Elements of Applied Bifurcation Theory, 3rd ed., Springer, New York, 2004. [Crossref], [Google Scholar]
X. Liu and D. Xiao, Complex dynamic behaviors of a discrete-time predator-prey system, Chaos Solitons Fractals 32 (2007), pp. 80–94. doi: 10.1016/j.chaos.2005.10.081 [Crossref], [Web of Science ®], [Google Scholar]
G. Livadiotis, L. Assas, B. Dennis, S. Elaydi and E. Kwessi, A discrete-time host–parasitoid model with an Allee effect, J. Biol. Dyn. 9 (2015), pp. 34–51. doi: 10.1080/17513758.2014.982219 [Taylor & Francis Online], [Google Scholar]
J. Ma and D. Xiao, Nonlinear dynamics of a mathematical model on action potential duration and calcium transient in paced cardiac cells, Discrete Contin. Dyn. Syst. Ser. B 18 (2013), pp. 2377–2396. doi: 10.3934/dcdsb.2013.18.2377 [Crossref], [Web of Science ®], [Google Scholar]
J.D. Murray, Mathematical Biology, Vol. 5, Springer, New York, 2002. [Google Scholar]
M.G. Neubert and M. Kot, The subcritical collapse of predator populations in discrete-time predator–prey models, Math. Biosci. 110 (1992), pp. 45–66. doi: 10.1016/0025-5564(92)90014-N [Crossref], [PubMed], [Web of Science ®], [Google Scholar]
A. Nicholson and V. Bailey, The balance of animal population, Proc. Zool. Soc. Lond. 3 (1935), pp. 551–598. doi: 10.1111/j.1096-3642.1935.tb01680.x [Crossref], [Google Scholar]
G.A. van Voorn, L. Hemerik, M. Boer and B. Kooi, Heteroclinic orbits indicate overexploitation in predator–prey systems with a strong Allee effect, Math. Biosci. 209 (2007), pp. 451–469. doi: 10.1016/j.mbs.2007.02.006 [Crossref], [PubMed], [Web of Science ®], [Google Scholar]

Journal of Biological Dynamics

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

Original Articles

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

ABSTRACT
Formulae display:?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.

1. Introduction

2. Local dynamics of system (3)

Published online:

Published online:

Published online:

Table 1. Number of fixed points in different parameter regions.

Published online:

3. Bifurcations of system (3)

Published online:

4. Boundedness of solutions and global attractors

Published online:

5. Numerical simulations

Published online:

Published online:

Published online:

Published online:

6. Conclusions

Disclosure statement

Funding

References

Related research

Information for

Open access

Opportunities

Help and information

Journal of Biological Dynamics

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

Original Articles

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

ABSTRACTFormulae display:?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.

1. Introduction

2. Local dynamics of system (3)

Published online:

Published online:

Published online:

Table 1. Number of fixed points in different parameter regions.

Published online:

3. Bifurcations of system (3)

Published online:

4. Boundedness of solutions and global attractors

Published online:

5. Numerical simulations

Published online:

Published online:

Published online:

Published online:

6. Conclusions

Acknowledgements

Disclosure statement

Additional information

Funding

References

Related research

Information for

Open access

Opportunities

Help and information

Keep up to date

ABSTRACT
Formulae display:?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.