# Multiple Bifurcations and Chaos in a Discrete Prey-Predator System with Generalized Holling III Functional Response

Authors: Xia Liu; Yanwei Liu; Qiaoping Li
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

A prey-predator system with the strong Allee effect and generalized Holling type III functional response is presented and discretized. It is shown that the combined influences of Allee effect and step size have an important effect on the dynamics of the system. The existences of Flip and Neimark-Sacker bifurcations and strange attractors and chaotic bands are investigated by using the center manifold theorem and bifurcation theory and some numerical methods.


---

## Body

## 1. Introduction

Lotka-Volterra model as a paradigm of the fundamental population models has been developed and widely applied in ecological sciences and other fields. Some linear and nonlinear functional responses are employed to describe the phenomena of predation including the Beddington-DeAngelis type [1, 2], Crowley-Martin type [3], and Holling types [4] that are more realistic and proper for modeling most prey-predator interactions [1–8]. Recently, a modified prey-predator model with a generalized Holling response function of type III was addressed [8]. The system exhibits more richer dynamics including the Hopf bifurcation of codimensions  1 and  2, Bogdanov-Takens bifurcation of codimensions  2 and  3.The researches of the Allee effect on biological population [9–23] imply that the Allee effect probably is a nonignorable factor in ecology, especially for the situation when population density (or size) is low (or small). That fitness of an individual in a small or sparse population decreases as the population density declines leads to Allee effect occurrence. Two most common continuous growth equations to express Allee effect for a single species are given [15–18]
(1)
x
˙
=
r
x
1
-
x
K
(
x
-
m
)
,
x
˙
=
r
1
-
x
K
-
l
x
+
b
x
,
which are called multiplicative Allee effect model [15, 16] and the additive Allee effect model [17, 18], respectively. Parameters m, l, and b
>
0 are the Allee thresholds, where l
x
/
(
x
+
b
) is termed as the strong Allee effect item if l
>
b
r and weak Allee effect item if l
<
b
r. Mathematically, the systems subjected to the Allee effect can be depicted by the bistability switch, which may lead to more complex dynamics. Thus, the consideration of prey-predator model subjected to an Allee effect is more realistic and has attracted renewed interest, which motivates us to make a modification to the system presented in [8] for further study:
(2)
x
˙
=
x
r
1
-
x
K
(
x
-
m
)
-
β
x
y
a
0
x
2
+
b
0
x
+
1
,
y
˙
=
y
-
d
0
+
c
0
x
2
a
0
x
2
+
b
0
x
+
1
,
where ˙
=
d
/
d
t. x and y denote the prey and predator densities at time t, respectively. r, K, m, β, a
0, c
0, and d
0 are positive constants. b
0 is positive or negative. In this paper, we will focus on the case b
0
>
0. In particular, r and K represent the intrinsic growth rate and carrying capacity of prey in the absence of predation, respectively. m
>
0 is so called the strong Allee threshold. The predator consumes the prey with functional response β
x
2
y
/
(
a
0
x
2
+
b
0
x
+
1
), known as Holling type III response and contributes to its growth rate c
0
x
2
y
/
(
a
0
x
2
+
b
0
x
+
1
). d
0 is the maximum death rate in the absence of prey. For simplicity, take nondimensional transformations x
=
K
X, y
=
Y
/
β, s
=
K
t to system (2) and rewrite X, Y, s as x, y, t; we have
(3)
x
˙
=
x
r
1
-
x
(
x
-
A
)
-
x
y
a
x
2
+
b
x
+
1
,
y
˙
=
y
-
d
+
c
x
2
a
x
2
+
b
x
+
1
,
where A
=
m
/
K, a
=
a
0
K
2, b
=
b
0
K, d
=
d
0
/
K, and c
=
c
0
K.Many studies of discrete-time models [19–38] were made and suggested that discrete-time models described by the difference equations are more appropriate than the continuous time models when populations have nonoverlapping generations or the number of populations is small. Moreover, dynamical patterns produced in discrete-time models are much richer than those observed in continuous time models.Applying the forward Euler scheme, we obtain the discrete-time form of system (3) as follows:
(4)
x
k
+
1
=
x
k
+
θ
x
k
r
1
-
x
k
x
k
-
A
-
x
k
y
k
a
x
k
2
+
b
x
k
+
1
≜
f
(
x
k
,
y
k
)
,
y
k
+
1
=
y
k
+
θ
y
k
-
d
+
c
x
k
2
a
x
k
2
+
b
x
k
+
1
≜
g
(
x
k
,
y
k
)
,
where θ
>
0 is the step size.The paper is organized as follows: in Section2 we show the existence and stability of fixed points of (4). In Section 3 we investigate the existence of flip bifurcation and the Neimark-Sacker bifurcation. In Section 4 we show some numerical simulations to illustrate our main results. Finally, a brief discussion of (4) is given to sum up our analysis.


## 2. Existence and Stability of Fixed Points

To obtain the fixed points of system (4), it is sufficient to solve the following algebraic equations:
(5)
x
r
1
-
x
(
x
-
A
)
-
x
y
a
x
2
+
b
x
+
1
=
0
,
y
-
d
+
c
x
2
a
x
2
+
b
x
+
1
=
0
.
After calculating, we have the following statement.Lemma 1.
System (4) always has the trivial fixed point E
0
=
(
0,0
) and the boundary fixed points E
1
=
(
1,0
), E
2
=
(
A
,
0
). If d
(
a
+
b
+
1
)
<
c
<
d
(
a
+
b
/
A
+
1
/
A
2
), (4) has a unique positive fixed point E
*
=
(
x
*
,
y
*
), where
(6)
x
*
=
b
d
+
b
2
d
2
+
4
d
c
-
a
d
2
c
-
a
d
,
y
*
=
r
1
-
x
*
(
x
*
-
A
)
(
a
x
*
2
+
b
x
*
+
1
)
x
*
.To analyze dynamical properties of (4), we need to compute the Jacobian matrix of (4) at any fixed point E
=
(
x
,
y
)
(7)
M
E
=
f
x
f
y
g
x
g
y
,
where
(8)
f
x
=
1
-
θ
r
3
x
2
-
2
x
A
+
1
+
A
-
θ
x
y
b
x
+
2
a
x
2
+
b
x
+
1
2
;
f
y
=
-
θ
x
2
a
x
2
+
b
x
+
1
;
g
x
=
θ
c
x
y
b
x
+
2
a
x
2
+
b
x
+
1
2
;
g
y
=
1
+
θ
-
d
+
c
x
2
a
x
2
+
b
x
+
1
.
The corresponding characteristic polynomial of M
E can be written as
(9)
F
(
λ
)
=
λ
2
+
p
λ
+
q
,
where p
=
-
f
x
-
g
y, q
=
f
x
g
y
-
f
y
g
x.Letλ
1 and λ
2 be two roots of equation F
(
λ
)
=
0. Clearly, the local stability of the fixed point E
(
x
,
y
) is determined by the modules of eigenvalues |
λ
1
| and |
λ
2
|; that is, if |
λ
1
|
<
1 and |
λ
2
|
<
1, then E
(
x
,
y
) is locally asymptotically stable and is termed a sink. If |
λ
1
|
>
1 and |
λ
2
|
>
1, then E
(
x
,
y
) is unstable and is termed a source. If |
λ
1
|
<
1 and |
λ
2
|
>
1 (or |
λ
1
|
>
1 and |
λ
2
|
<
1), then E
(
x
,
y
) is termed a saddle. If |
λ
1
|
=
1 or |
λ
2
|
=
1, then E
(
x
,
y
) is termed nonhyperbolic. These can be stated by the following lemma.Lemma 2 (see [20]).
SupposeF
(
1
)
>
0; then (i)
|
λ
1
|
<
1 and |
λ
2
|
<
1 if and only if F
(
-
1
)
>
0 and q
<
1;
(ii)
λ
1
=
-
1 and |
λ
2
|
≠
1 if and only if F
(
-
1
)
=
0 and p
≠
0
,
2;
(iii)
|
λ
1
|
>
1 and |
λ
2
|
>
1 if and only if F
(
-
1
)
>
0 and q
>
1;
(iv)
|
λ
1
|
>
1 and |
λ
2
|
<
1 if and only if F
(
-
1
)
<
0.A simple computation show that two eigenvalues ofM
E
0 are λ
1
=
1
-
θ
r
A and λ
2
=
1
-
θ
d. Also note that F
(
1
)
=
θ
2
r
A
d
>
0, F
(
-
1
)
=
(
θ
d
-
2
)
(
θ
r
A
-
2
), p
=
θ
(
A
r
+
d
)
-
2, and q
=
(
θ
r
A
-
1
)
(
θ
d
-
1
). Thus, from Lemma 2, we have the following.Lemma 3.
The fixed pointE
0
=
(
0,0
) of model (4) always exists; moreover, (i)
E
0 is a source if (
θ
d
-
2
)
(
θ
r
A
-
2
)
>
0, θ
r
A
d
>
d
+
r
A;
(ii)
E
0 is a sink if (
θ
d
-
2
)
(
θ
r
A
-
2
)
>
0, θ
r
A
d
<
d
+
r
A;
(iii)
E
0 is a saddle if (
θ
d
-
2
)
(
θ
r
A
-
2
)
<
0;
(iv)
E
0 is nonhyperbolic if (
θ
d
-
2
)
(
θ
r
A
-
2
)
=
0.Remark 4.
Define(10)
S
E
0
1
=
A
,
a
,
b
,
c
,
d
,
θ
:
θ
=
2
d
,
A
≠
d
r
,
S
E
0
2
=
A
,
a
,
b
,
c
,
d
,
θ
:
A
=
2
r
θ
,
θ
≠
2
d
.
From Lemma 3, we see that if (
A
,
a
,
b
,
c
,
d
,
θ
)
∈
S
E
0
1 or S
E
0
2, then one of the two eigenvalues of F
(
λ
)
=
0 is -
1 and the other is neither 1 nor -
1, which implies that system (4) may undergo the flip bifurcation at the fixed point E
0 when parameters (
A
,
a
,
b
,
c
,
d
,
θ
) vary in a small neighborhood of S
E
0
1 or S
E
0
2.
Similarly, the eigenvalues of JacobianM
E
1 are λ
1
=
1
-
r
θ
(
1
-
A
) and λ
2
=
1
+
θ
(
c
-
d
(
a
+
b
+
1
)
)
/
(
a
+
b
+
1
). By Lemma 2, we obtain the following results.Lemma 5.
ForE
1
=
(
1,0
), one has the following:(i)
E
1 is a sink if 0
<
θ
<
2
/
r
(
1
-
A
), (
a
+
b
+
1
)
d
-
2
/
θ
<
c
<
d
(
a
+
b
+
1
);
(ii)
E
1 is a source if one of the following two cases holds:
(1)
θ
>
2
/
r
(
1
-
A
), 0
<
c
<
(
a
+
b
+
1
)
d
-
2
/
θ with d
>
2
/
θ;
(2)
θ
>
2
/
r
(
1
-
A
), c
>
d
(
a
+
b
+
1
);
(iii)
E
1 is a saddle if one of the following three cases holds:
(1)
0
<
θ
<
2
/
r
(
1
-
A
), 0
<
c
<
(
a
+
b
+
1
)
d
-
2
/
θ with d
>
2
/
θ;
(2)
0
<
θ
<
2
/
r
(
1
-
A
), c
>
d
(
a
+
b
+
1
);
(3)
θ
>
2
/
r
(
1
-
A
), (
a
+
b
+
1
)
d
-
2
/
θ
<
c
<
d
(
a
+
b
+
1
);
(iv)
E
1 is nonhyperbolic if one of the following three cases holds:
(1)
θ
=
2
/
r
(
1
-
A
);
(2)
c
=
(
a
+
b
+
1
)
d
-
2
/
θ with θ
>
2
/
d;
(3)
c
=
d
(
a
+
b
+
1
).Remark 6.
Lemma5(iv) implies that there are two parameter surfaces defined, respectively, by
(11)
S
E
1
1
=
(
A
,
a
,
b
,
c
,
d
,
θ
)
:
θ
=
2
r
1
-
A
,
c
≠
(
a
+
b
+
1
)
d
-
2
θ
,
c
≠
d
(
a
+
b
+
1
)
2
r
1
-
A
,
S
E
1
2
=
2
r
1
-
A
(
A
,
a
,
b
,
c
,
d
,
θ
)
:
c
=
(
a
+
b
+
1
)
d
-
2
θ
,
θ
≠
2
r
1
-
A
,
d
>
2
θ
,
such that one of the eigenvalues of M
E
1 is -
1 and the other is neither 1 nor -
1 when (
A
,
a
,
b
,
c
,
d
,
θ
)
∈
S
E
1
1 or S
E
1
2. Thus, system (4) may undergo a flip bifurcation when (
A
,
a
,
b
,
c
,
d
,
θ
) changes in a small neighborhood of S
E
1
1 or S
E
1
2.
After calculating, we also obtain the eigenvalues of the JacobianM
E
2 are λ
1
=
1
+
θ
A
r
(
1
-
A
) and λ
2
=
1
+
θ
-
d
+
c
A
2
/
(
a
A
2
+
b
A
+
1
). Note that 0
<
A
<
1 so λ
1
>
1. Obviously, if |
λ
2
|
>
1, then E
2 is a source. If |
λ
2
|
<
1, then E
2 is a saddle. If |
λ
2
|
=
1, then E
2 is nonhyperbolic. By calculating, one can obtain the following.Lemma 7.
The dynamical behaviors of the fixed pointE
2
=
(
A
,
0
) can be stated as follows. (i)
It is a source ifc
>
d
(
a
A
2
+
b
A
+
1
)
/
A
2 or 0
<
c
<
(
d
θ
-
2
)
(
a
A
2
+
b
A
+
1
)
/
A
2
θ with θ
>
2
/
d.
(ii)
It is a saddle if(
d
θ
-
2
)
(
a
A
2
+
b
A
+
1
)
/
A
2
θ
<
c
<
d
(
a
A
2
+
b
A
+
1
)
/
A
2.
(iii)
It is nonhyperbolic ifc
=
d
(
a
A
2
+
b
A
+
1
)
/
A
2 or c
=
(
d
θ
-
2
)
(
a
A
2
+
b
A
+
1
)
/
A
2
θ with θ
>
2
/
d.Now, we discuss local stability of the fixed pointE
*
(
x
*
,
y
*
). After calculating, we get
(12)
p
=
-
f
x
(
x
*
,
y
*
)
-
g
y
(
x
*
,
y
*
)
=
-
2
-
θ
(
r
C
-
B
)
;
q
=
1
+
θ
(
r
C
-
B
)
+
d
B
θ
2
,
where B
=
x
*
y
*
(
b
x
*
+
2
)
/
(
a
x
*
2
+
b
x
*
+
1
)
2
>
0 and C
=
-
3
x
*
2
+
2
(
1
+
A
)
x
*
-
A. The discriminant of equation F
(
λ
)
=
0 can be calculated Δ
=
θ
2
(
(
r
C
-
B
)
2
-
4
d
B
). Moreover, F
(
1
)
=
d
B
θ
2
>
0 and
(13)
F
(
-
1
)
=
d
B
θ
2
+
2
(
r
C
-
B
)
θ
+
4
.
Regarding q and F
(
-
1
) as two functions of θ, then the discriminants of equations q
(
θ
)
=
0 and F
(
-
1
)
(
θ
)
=
0 are (
r
C
-
B
)
2
-
4
d
B and 4
(
(
r
C
-
B
)
2
-
4
d
B
), respectively.Clearly, the sign ofΔ is determined by (
r
C
-
B
)
2
-
4
d
B. Therefore we have the following.Lemma 8.
Assumed
(
a
+
b
+
1
)
<
c
<
d
(
a
+
b
/
A
+
1
/
A
2
), then (i)
E
* is a source if one of the following three cases holds:
(1)
B
⩽
r
C;
(2)
B
>
r
C, (
r
C
-
B
)
2
⩽
4
d
B, θ
>
(
B
-
r
C
)
/
d
B;
(3)
B
>
r
C, (
r
C
-
B
)
2
>
4
d
B, θ
>
(
B
-
r
C
+
(
B
-
r
C
)
2
-
4
d
B
)
/
d
B;
(ii)
E
* is a sink if one of the following two cases holds:
(1)
B
>
r
C, (
r
C
-
B
)
2
⩽
4
d
B, 0
<
θ
<
(
B
-
r
C
)
/
d
B;
(2)
B
>
r
C, (
r
C
-
B
)
2
>
4
d
B, 0
<
θ
<
(
B
-
r
C
-
(
B
-
r
C
)
2
-
4
d
B
)
/
d
B;
(iii)
E
* is a saddle if (
B
-
r
C
-
(
B
-
r
C
)
2
-
4
d
B
)
/
d
B
<
θ
<
(
B
-
r
C
+
(
B
-
r
C
)
2
-
4
d
B
)
/
d
B and B
>
r
C, (
r
C
-
B
)
2
>
4
d
B;
(iv)
E
* is nonhyperbolic if one of the following three cases holds:
(1)
B
>
r
C, (
r
C
-
B
)
2
⩽
4
d
B, θ
=
(
B
-
r
C
)
/
d
B;
(2)
B
>
r
C, (
r
C
-
B
)
2
>
4
d
B, θ
=
(
B
-
r
C
-
(
B
-
r
C
)
2
-
4
d
B
)
/
d
B;
(3)
B
>
r
C, (
r
C
-
B
)
2
>
4
d
B, θ
=
(
B
-
r
C
+
(
B
-
r
C
)
2
-
4
d
B
)
/
d
B.Proof.
ThatF
(
1
)
>
0 always holds. Set λ
1 and λ
2 are eigenvalues of matrix M
E
*. We will consider two cases to prove Lemma 8.
Case 1. Assume B
⩽
r
C, then F
(
-
1
)
>
0 and q
>
1 which yield |
λ
1
|
>
1 and |
λ
2
|
>
1.
Case 2. Assume B
>
r
C. One can see that if (
r
C
-
B
)
2
⩽
4
d
B, then q
⩾
0, F
(
-
1
)
⩾
0; moreover F
(
-
1
)
=
0 if and only if (
r
C
-
B
)
2
=
4
d
B and θ
=
(
B
-
r
C
)
/
d
B. By comparing q and 1, we obtain that if θ
>
(
B
-
r
C
)
/
d
B, then F
(
-
1
)
>
0 and q
>
1 which lead to the fact that E
* is a source. If 0
<
θ
<
(
B
-
r
C
)
/
d
B the F
(
-
1
)
>
0 and q
<
1 which lead to the fact that E
* is a sink. If θ
=
(
B
-
r
C
)
/
d
B the F
(
-
1
)
=
0 and q
=
1 which lead to |
λ
1
|
=
1 and |
λ
2
|
=
1, that is, E
* is nonhyperbolic.
In Case  2, the fact that the assumption(
r
C
-
B
)
2
>
4
d
B leads to that the equation F
(
-
1
)
=
0 has two distinct positive roots θ
1
=
(
B
-
r
C
-
(
B
-
r
C
)
2
-
4
d
B
)
/
d
B and θ
2
=
(
B
-
r
C
+
(
B
-
r
C
)
2
-
4
d
B
)
/
d
B so that F
(
-
1
)
>
0 when 0
<
θ
<
θ
1 or θ
>
θ
2, F
(
-
1
)
<
0 when θ
1
<
θ
<
θ
2. Also note that q
>
1 if θ
>
(
B
-
r
C
)
/
d
B, q
<
1 if 0
<
θ
<
(
B
-
r
C
)
/
d
B. Thus, we derive that F
(
-
1
)
>
0 and q
>
1 if θ
>
θ
2, F
(
-
1
)
>
0 and q
<
1 if 0
<
θ
<
θ
1, F
(
-
1
)
=
0 if θ
=
θ
1 or θ
2. Following above discussion and Lemma 2, the proof of Lemma 8 is complete.From above analysis and Lemma8, we have two parameter surfaces:
(14)
H
*
=

a
+
b
A
+
1
A
2
(
A
,
a
,
b
,
c
,
d
,
θ
)
:
θ
=
B
-
r
C
d
B
,
d
(
a
+
b
+
1
)
<
c
<
d
a
+
b
A
+
1
A
2
,
0
<
B
-
r
C
<
2
d
B
a
+
b
A
+
1
A
2
;
S
E
*
=

a
+
b
A
+
1
A
2
(
A
,
a
,
b
,
c
,
d
,
θ
)
:
θ
=
θ
1
,
d
(
a
+
b
+
1
)
<
c
<
d
a
+
b
A
+
1
A
2
,
B
-
r
C
>
2
d
B
,
θ
(
B
-
r
C
)
≠
2,4
a
+
b
A
+
1
A
2
.From Lemma8, we can see that if (
A
,
a
,
b
,
c
,
d
,
θ
)
∈
H
*, then two eigenvalues of M
E
* are a pair of conjugate complex numbers and satisfy |
λ
1
|
=
|
λ
2
|
=
1. Thus, when (
A
,
a
,
b
,
c
,
d
,
θ
) vary in a small neighborhood of H
* model (4) undergoes a Neimark-Sacker bifurcation from the fixed point E
*. Similarly, if (
A
,
a
,
b
,
c
,
d
,
θ
)
∈
S
E
*, then one of two eigenvalues of M
E
* is -
1 and the other is neither 1 nor -
1; that is, when (
A
,
a
,
b
,
c
,
d
,
θ
) vary in a small neighborhood of S
E
*, a flip bifurcation from the fixed point E
* can occur.


## 3. Analysis of Bifurcation

On the basis of the analysis in Section2, by choosing parameter θ as a bifurcation parameter, we mainly investigate the flip bifurcation and Neimark-Sacker bifurcation of positive fixed point E
(
x
*
,
y
*
) by using the center manifold theorem and bifurcation theory of [39–41] in this section.Giving a perturbationθ
* (|
θ
*
|
≪
1) to parameter θ we consider a perturbation of model (4) as follows:
(15)
x
k
+
1
=
x
k
+
(
θ
+
θ
*
)
x
k
×
r
1
-
x
k
(
x
k
-
A
)
-
x
k
y
k
a
x
k
2
+
b
x
k
+
1
,
y
k
+
1
=
y
k
+
θ
+
θ
*
y
k
-
d
+
c
x
k
2
a
x
k
2
+
b
x
k
+
1
.LetE
(
x
,
y
) be a feasible fixed point of model (4). We translate E
(
x
,
y
) to the origin by using transformations X
k
=
x
k
-
x, Y
k
=
y
k
-
y. Then we have (rewrite X
k, Y
k as x
k, y
k)
(16)
x
k
+
1
=
x
k
+
(
θ
+
θ
*
)
(
x
k
+
x
)
×

(
x
k
+
x
)
(
y
k
+
y
)
a
x
k
+
x
2
+
r
1
-
x
k
-
x
(
x
k
+
x
-
A
)
-
(
x
k
+
x
)
(
y
k
+
y
)
a
x
k
+
x
2
+
b
(
x
k
+
x
)
+
1
,
y
k
+
1
=
y
k
+
θ
+
θ
*
y
k
+
y
×
-
d
+
c
x
k
+
x
2
a
x
k
+
x
2
+
b
(
x
k
+
x
)
+
1
.
Expanding model (16) as Taylor series at (
x
k
,
y
k
,
θ
*
)
=
(
0,0
,
0
) to the second order, then we have
(17)
x
k
+
1
=
a
100
x
k
+
a
010
y
k
+
a
200
x
k
2
+
a
110
x
k
y
k
+
a
101
x
k
θ
*
+
a
011
y
k
θ
*
+
O
(
x
k
,
y
k
,
θ
*
4
)
;
y
k
+
1
=
b
100
x
k
+
b
010
y
k
+
b
200
x
k
2
+
b
110
x
k
y
k
+
b
101
x
k
θ
*
+
b
011
y
k
θ
*
+
O
(
x
k
,
y
k
,
θ
*
4
)
,
where
(18)
a
100
=
1
+
θ
r
(
2
x
-
A
-
3
x
2
+
2
x
A
)
-
θ
x
y
(
b
x
+
2
)
a
x
2
+
b
x
+
1
2
;
a
010
=
-
x
2
θ
a
x
2
+
b
x
+
1
;
a
200
=
θ
r
(
1
+
A
-
3
x
)
+
θ
y
(
3
a
x
2
+
a
x
3
b
-
1
)
a
x
2
+
b
x
+
1
3
;
a
110
=
-
θ
x
(
b
x
+
2
)
a
x
2
+
b
x
+
1
2
;
a
101
=
a
100
-
1
θ
;
a
011
=
-
x
2
a
x
2
+
b
x
+
1
;
b
100
=
θ
c
x
y
(
b
x
+
2
)
a
x
2
+
b
x
+
1
2
;
b
010
=
1
+
θ
-
d
+
c
x
2
a
x
2
+
b
x
+
1
;
b
200
=
-
θ
y
c
(
a
b
x
3
+
3
a
x
2
-
1
)
a
x
2
+
b
x
+
1
3
;
b
110
=
θ
c
x
(
b
x
+
2
)
a
x
2
+
b
x
+
1
2
;
b
101
=
c
x
y
(
b
x
+
2
)
a
x
2
+
b
x
+
1
2
;
b
011
=
-
d
+
c
x
2
a
x
2
+
b
x
+
1
.The flip bifurcation of model (4) at E
* will be investigated firstly when (
A
,
a
,
b
,
c
,
d
,
θ
) varies in the small neighborhood of S
E
*. Similar arguments can be applied to the other cases S
E
0
1, S
E
0
2, S
E
1
1, S
E
1
2, and S
E
2.Let(
A
,
a
,
b
,
c
,
d
,
θ
)
∈
S
E
* arbitrarily. From the discussion in Section 2, we can see model (4) has a unique positive fixed point E
*; its eigenvalues are λ
1
=
-
1 and -
1
<
λ
2
<
0. In this case, x
=
x
*, y
=
y
*, and θ
=
θ
1. Also note that p
=
-
a
100
-
b
010
=
-
a
100
-
1 and λ
1
+
λ
2
=
-
p so λ
2
=
2
+
a
100. Moreover, there is a transformation
(19)
x
k
y
k
=
a
010
a
010
-
1
-
a
100
λ
2
-
a
100
u
k
v
k
so that model (17) becomes
(20)
u
k
+
1
=
-
u
k
+
f
1
(
x
k
,
y
k
,
θ
*
)
+
O
(
x
k
,
y
k
,
θ
*
4
)
,
v
k
+
1
=
λ
2
v
k
+
f
2
(
x
k
,
y
k
,
θ
*
)
+
O
(
x
k
,
y
k
,
θ
*
4
)
,
where
(21)
f
1
x
k
,
y
k
,
θ
*
=
2
a
200
-
a
010
b
200
a
010
1
+
λ
2
x
k
2
+
2
a
110
-
a
010
b
110
a
010
1
+
λ
2
x
k
y
k
+
2
a
101
-
a
010
b
101
a
010
(
1
+
λ
2
)
x
k
θ
*
+
2
a
011
a
010
(
1
+
λ
2
)
y
k
θ
*
,
f
2
x
k
,
y
k
,
θ
*
=
a
010
b
200
-
p
a
200
a
010
1
+
λ
2
x
k
2
+
a
010
b
110
-
p
a
110
a
010
1
+
λ
2
x
k
y
k
+
a
010
b
101
-
p
a
101
a
010
(
1
+
λ
2
)
x
k
θ
*
-
p
a
011
a
010
(
1
+
λ
2
)
y
k
θ
*
.By the center manifold theorem [41], we can get the approximate representation of the center manifold W
c
(
0,0
,
0
) of model (20) at the fixed point (
0,0
) in a small neighborhood of θ
1 as follows:
(22)
W
c
0,0
,
0
=

u
k
,
θ
*
3
(
u
k
,
v
k
,
θ
*
)
∈
R
3
:
v
k
=
η
0
u
k
2
+
η
1
u
k
θ
*
+
O
u
k
,
θ
*
3
,
where
(23)
η
0
=
a
010
(
a
010
b
200
-
p
a
200
)
1
-
λ
2
2
+
p
(
a
010
b
110
-
p
a
110
)
1
-
λ
2
2
;
η
1
=
-
a
010
b
101
-
p
a
101
1
+
λ
2
2
+
p
2
a
011
a
010
1
+
λ
2
2
.
Furthermore, on the center manifold W
c
(
0,0
,
0
), we have
(24)
x
k
=
a
010
(
u
k
+
η
0
u
k
2
+
η
1
u
k
θ
*
)
+
O
(
u
k
,
θ
*
3
)
;
y
k
=
p
u
k
+
2
(
η
0
u
k
2
+
η
1
u
k
θ
*
)
+
O
(
u
k
,
θ
*
3
)
;
x
k
2
=
a
010
2
(
3
u
k
2
+
2
η
0
u
k
3
+
2
η
1
u
k
2
θ
*
)
+
O
(
u
k
,
θ
*
4
)
;
x
k
y
k
=
-
a
010
u
k
2
(
-
2
η
0
u
k
-
u
k
η
0
p
-
p
-
2
η
1
θ
*
-
η
1
θ
*
p
)
+
O
(
u
k
,
θ
*
4
)
;
x
k
θ
*
=
a
010
(
u
k
+
η
0
u
k
2
+
η
1
u
k
θ
*
)
θ
*
+
O
(
u
k
,
θ
*
4
)
;
y
k
θ
*
=
(
p
u
k
+
2
η
0
u
k
2
+
2
η
1
u
k
θ
*
)
θ
*
+
O
(
u
k
,
θ
*
4
)
.
Therefore, the map F which is model (4) restricted to the center manifold W
c
(
0,0
,
0
) takes the form
(25)
F
u
k
=
-
u
k
+
h
1
u
k
2
+
h
2
u
k
θ
*
+
h
3
u
k
3
+
h
4
u
k
2
θ
*
+
h
5
u
k
θ
*
2
+
O
(
u
k
,
θ
*
4
)
,
where
(26)
h
1
=
p
(
2
a
110
-
a
010
b
110
)
1
+
λ
2
+
a
010
(
2
a
200
-
a
010
b
200
)
1
+
λ
2
;
h
2
=
2
p
a
011
a
010
(
1
+
λ
2
)
+
2
a
101
-
a
010
b
101
1
+
λ
2
;
h
3
=
2
η
0
(
1
+
p
)
(
2
a
110
-
a
010
b
110
)
1
+
λ
2
+
2
η
0
a
010
(
2
a
200
-
a
010
b
200
)
1
+
λ
2
;
h
4
=
η
1
2
+
p
2
a
110
-
a
010
b
110
1
+
λ
2
+
4
η
0
a
011
a
010
1
+
λ
2
+
2
η
1
a
010
(
2
a
200
-
a
010
b
200
)
1
+
λ
2
+
η
0
(
2
a
101
-
a
010
b
101
)
1
+
λ
2
;
h
5
=
4
η
1
a
011
a
010
(
1
+
λ
2
)
+
η
1
(
2
a
101
-
a
010
b
101
)
1
+
λ
2
.
In order to undergo a flip bifurcation for map (25), we require that two discriminatory quantities α
1 and α
2 are not zero, where
(27)
α
1
=
∂
2
F
∂
u
k
∂
θ
*
+
1
2
∂
F
∂
θ
*
∂
2
F
∂
u
k
2
(
0,0
)
=
h
2
≠
0
;
α
2
=
1
6
∂
3
F
∂
u
k
3
+
1
2
∂
2
F
∂
u
k
2
2
(
0,0
)
=
h
3
+
h
1
2
.
From above analysis and the theorem in [40], we have the following.Theorem 9.
Ifα
2
≠
0, then model (4) undergoes a flip bifurcation at equilibrium E
*
(
x
*
,
y
*
) when the parameter θ varies in the small neighborhood of θ
1. Moreover, if α
2
>
0 (resp. α
2
<
0), then the period 2 points that bifurcate from E
* are stable (resp., unstable).Next, we discuss the Neimark-Sacker bifurcation ofE
* when parameters (
A
,
a
,
b
,
c
,
d
,
θ
) vary in the small neighborhood of H
*. Taking θ
=
(
B
-
r
C
)
/
d
B (for simplicity, denoting (
B
-
r
C
)
/
d
B by θ
0) and choosing θ
* as a perturbation of bifurcation parameter, a perturbed form of model (4) is obtained as follows:
(28)
x
k
+
1
=
x
k
+
θ
0
+
θ
*
x
k
×
r
1
-
x
k
(
x
k
-
A
)
-
x
k
y
k
a
x
k
2
+
b
x
k
+
1
,
y
k
+
1
=
y
k
+
(
θ
0
+
θ
*
)
y
k
-
d
+
c
x
k
2
a
x
k
2
+
b
x
k
+
1
,
where |
θ
*
|
≪
1 which is a small perturbation.We transform the unique positive fixed pointE
* into the origin by using X
k
=
x
k
-
x
*, Y
k
=
y
k
-
y
*:
(29)
X
k
+
1
=
a
10
X
k
+
a
01
Y
k
+
a
11
X
k
Y
k
+
a
20
X
k
2
+
a
02
Y
k
2
+
a
30
X
k
3
+
a
21
X
k
2
Y
k
+
X
k
,
Y
k
3
,
Y
k
+
1
=
b
10
X
k
+
b
01
Y
k
+
b
11
X
k
Y
k
+
b
20
X
k
2
+
b
02
Y
k
2
+
b
30
X
k
3
+
b
21
X
k
2
Y
k
+
X
k
,
Y
k
3
,
where
(30)
a
10
=
1
+
θ
0
+
θ
*
r
2
x
*
-
A
-
3
x
*
2
+
2
x
*
A
-
θ
0
+
θ
*
x
*
y
*
b
x
*
+
2
a
x
*
2
+
b
x
*
+
1
2
;
a
01
=
-
d
θ
0
+
θ
*
c
;
a
11
=
-
θ
0
+
θ
*
x
*
b
x
*
+
2
a
x
*
2
+
b
x
*
+
1
2
;
a
20
=
θ
0
+
θ
*
r
1
+
A
-
3
x
*
+
θ
0
+
θ
*
y
*
3
a
x
*
2
+
a
x
*
3
b
-
1
a
x
*
2
+
b
x
*
+
1
3
;
a
30
=
-
θ
0
+
θ
*
y
*
a
x
*
2
-
1
a
x
*
2
b
+
4
a
x
*
+
b
a
x
*
2
+
b
x
*
+
1
4
-
θ
0
+
θ
*
r
;
a
21
=
θ
0
+
θ
*
a
b
x
*
3
+
3
a
x
*
2
-
1
a
x
*
2
+
b
x
*
+
1
3
;
b
10
=
θ
0
+
θ
*
c
x
*
y
*
b
x
*
+
2
a
x
*
2
+
b
x
*
+
1
2
;
b
01
=
1
;
b
20
=
-
θ
0
+
θ
*
y
*
c
a
b
x
*
3
+
3
a
x
*
2
-
1
a
x
*
2
+
b
x
*
+
1
3
;
b
11
=
θ
0
+
θ
*
c
x
*
b
x
*
+
2
a
x
*
2
+
b
x
*
+
1
2
;
b
30
=
-
c
a
30
;
b
21
=
-
c
a
21
.
The characteristic equation associated with linearization system of model (29) at (
0,0
) is
(31)
λ
2
+
p
θ
*
λ
+
q
θ
*
=
0
,
where
(32)
p
θ
*
=
-
2
-
θ
0
+
θ
*
r
C
-
B
,
q
(
θ
*
)
=
1
+
θ
0
+
θ
*
r
C
-
B
+
d
B
θ
0
+
θ
*
2
.
When parameters (
A
,
a
,
b
,
c
,
d
,
θ
)
∈
H
* and θ
* vary in a small neighborhood of θ
*
=
0, the roots of the characteristic equation (31) are a pair of complex conjugate numbers λ
1 and λ
2 denoted by
(33)
λ
1,2
=
1
+
θ
0
+
θ
*
2
r
C
-
B
±
i
4
d
B
-
r
C
-
B
2
,
(34)
λ
1,2
=
q
θ
*
=
1
+
θ
0
+
θ
*
r
C
-
B
+
d
B
θ
0
+
θ
*
2
.
One can see that
(35)
d
λ
1,2
d
θ
*
θ
*
=
0
=
B
-
r
C
2
>
0
.On the other hand, it is required thatλ
1,2
k
≠
1, m
=
1,2
,
3,4 when θ
*
=
0, which is equivalent to p
(
0
)
≠
±
2,0
,
-
1. Also note that p
(
0
)
2
-
4
q
(
0
)
<
0 and q
(
0
)
=
1 when (
A
,
a
,
b
,
c
,
d
,
θ
)
∈
H
*. Thus, p
(
0
)
2
<
4, which indicates p
(
0
)
≠
±
2. We only need to require that p
(
0
)
≠
0,1; that is,
(36)
r
C
-
B
2
≠
2
d
B
,
3
d
B
,
which implies that λ
1,2 defined by (33) do not lie in the intersection of the unit circle with the coordinate axes when θ
*
=
0.Next we study the normal form of model (29) when θ
*
=
0.Letα
=
Re
(
λ
1
), β
=
Im
⁡
(
λ
1
), and
(37)
T
=
a
01
0
α
-
a
10
-
β
=
a
01
0
-
θ
0
(
r
C
-
B
)
2
-
θ
0
4
d
B
-
r
C
-
B
2
2
and perform the transformations X
k
=
a
01
u
k, Y
k
=
(
α
-
a
10
)
u
k
-
β
v
k to (29); we can obtain
(38)
u
k
+
1
=
α
u
k
-
β
v
k
+
G
1
u
k
,
v
k
,
v
k
+
1
=
β
u
k
+
α
v
k
+
G
2
u
k
,
v
k
,
where
(39)
G
1
u
k
,
v
k
=
a
20
a
01
-
a
11
a
10
-
b
01
2
+
a
02
a
10
-
b
01
2
4
a
01
u
k
2
+
a
02
β
2
a
01
v
k
2
-
β
a
11
a
01
+
a
02
b
01
-
a
10
2
a
01
u
k
v
k
+
a
01
2
a
30
a
01
-
a
21
a
10
-
b
01
2
u
k
3
-
a
01
a
21
β
u
k
2
v
k
+
O
(
u
k
,
v
k
4
)
;
G
2
(
u
k
,
v
k
)
=
-
1
8
a
01
β
×
b
20
a
01
3
(
a
10
-
b
01
)
mmmmi
×
a
02
a
10
-
b
01
2
-
2
a
01
(
a
11
-
b
02
)
(
a
10
-
b
01
)
mmmmmm
+
4
a
01
2
(
a
20
-
b
11
)
+
8
b
20
a
01
3
u
k
2
+
1
2
a
01

a
10
-
b
01
2
a
01
(
a
11
-
2
b
02
)
(
a
10
-
b
01
)
mmmimimm
-
a
02
a
10
-
b
01
2
+
2
b
11
a
01
2
u
k
v
k
-
β
2
a
01
b
02
+
a
02
a
10
-
b
01
2
a
01
v
k
2
-
a
01
4
β
a
10
-
b
01
2
2
a
01
(
a
30
-
b
21
)
(
a
10
-
b
01
)
-
a
21
a
10
-
b
01
2
+
4
b
30
a
01
2
u
k
3
+
[
a
21
(
a
10
-
b
01
)
+
2
b
21
a
01
]
a
01
2
u
k
2
v
k
.Let(40)
ξ
11
=
1
4
∂
2
G
1
∂
u
k
2
+
∂
2
G
1
∂
v
k
2
+
i
∂
2
G
2
∂
u
k
2
-
∂
2
G
2
∂
v
k
2
(
u
k
,
v
k
)
=
(
0,0
)
,
ξ
20
=
1
8
∂
2
G
1
∂
u
k
2
-
∂
2
G
1
∂
v
k
2
+
2
∂
2
G
2
∂
u
k
∂
v
k
+
i
∂
2
G
2
∂
u
k
2
-
∂
2
G
2
∂
v
k
2
-
2
∂
2
G
1
∂
u
k
∂
v
k
(
u
k
,
v
k
)
=
(
0,0
)
,
ξ
02
=
1
8
∂
2
G
1
∂
u
k
2
-
∂
2
G
1
∂
v
k
2
-
2
∂
2
G
2
∂
u
k
∂
v
k
+
i
∂
2
G
2
∂
u
k
2
-
∂
2
G
2
∂
v
k
2
+
2
∂
2
G
1
∂
u
k
∂
v
k
(
u
k
,
v
k
)
=
(
0,0
)
,
ξ
21
=
1
16
∂
3
G
1
∂
u
k
3
+
∂
3
G
1
∂
u
k
∂
v
k
2
+
∂
3
G
2
∂
u
k
2
∂
v
k
+
∂
3
G
2
∂
v
k
3
m
m
+
i
∂
3
G
2
∂
u
k
3
+
∂
3
G
2
∂
u
k
∂
v
k
2
-
∂
3
G
1
∂
u
k
2
∂
v
k
-
∂
3
G
1
∂
v
k
3
(
u
k
,
v
k
)
=
(
0,0
)
.
By [39], we know that if the following discriminatory quantity holds, then model (38) undergoes a Neimark-Sacker bifurcation.
(41)
μ
*
=
-
Re
1
-
2
λ
λ
-
2
1
-
λ
ξ
20
ξ
11
-
ξ
11
2
2
-
ξ
02
2
+
Re
(
λ
-
ξ
21
)
≠
0
,
where
(42)
λ
=
α
+
i
β
,
ξ
11
2
=
1
16

∂
2
G
1
∂
u
k
2
+
∂
2
G
1
∂
v
k
2

2

+

∂
2
G
2
∂
u
k
2
+
∂
2
G
2
∂
v
k
2

2
u
k
,
v
k
=
0,0
,
ξ
02
2
=
1
64
∂
2
G
1
∂
u
k
2
-
∂
2
G
1
∂
v
k
2
-
2
∂
2
G
2
∂
u
k
∂
v
k
2
+
∂
2
G
2
∂
u
k
2
-
∂
2
G
2
∂
v
k
2
+
2
∂
2
G
1
∂
u
k
∂
v
k
2
(
u
k
,
v
k
)
=
(
0,0
)
.
Form above analysis and the theorem in [39], we have the following.Theorem 10.
Let(
r
C
-
B
)
2
≠
2
d
B
,
3
d
B and (41) hold. Then system (4) undergoes Neimark-Sacker bifurcation at the positive fixed point E
*
(
x
*
,
y
*
) when the parameter θ varies in the small neighborhood of θ
0. Moreover, if μ
*
<
0 (resp., μ
*
>
0), then an attracting (resp., repelling) invariant closed curve bifurcates from the fixed point E
* for θ
>
θ
0 (resp., θ
<
θ
0).


## 4. Numerical Simulations

To confirm the above theoretical analysis, in this section we present the bifurcation diagrams and phase portraits for system (4). θ is chosen as the bifurcation parameter to verify the existence of flip bifurcation and Neimark-Sacker bifurcation.Firstly, letr
=
0.5, A
=
0.2, a
=
0.4, b
=
5, d
=
0.1, c
=
1.168, and θ
=
1.290556794. One can verify that all conditions in H
* are satisfied. That is, if θ varies in a small neighborhood of |
θ
-
1.290556794
|
<
ϵ (0
<
ϵ
≪
1) and the other parameters are kept fixed, then system (4) undergoes a Neimark-Sacker bifurcation, which is illustrated in Figure 1.(a) The bifurcation diagram of system (4) with θ varying in (
0.2,2
). The system undergoes a Neimark-Sacker bifurcation when θ
≈
1.290556794 (denoted by red point). (b) A stable fixed point E
*
(
0.5928,0.5537
) exists when θ
=
1.2 which occurs before the bifurcation. (c) An invariant closed curve around the fixed point E
* created after the bifurcation which exists for θ
=
1.3. (d) A chaotic attractor occurs when θ
=
6.6.
(a)
(b)
(c)
(d)Letr
=
10, A
=
0.5, a
=
3, b
=
9, d
=
0.1, and c
=
6.62. By θ
1
=
(
B
-
r
C
-
(
B
-
r
C
)
2
-
4
d
B
)
/
d
B, we obtain θ
1
≈
0.456; that is, (
r
,
A
,
a
,
b
,
d
,
c
,
θ
)
∈
S
E
* when θ
=
θ
1, which leads to α
1
≠
0 and α
2
>
0. Moreover, E
* is stable when θ
<
θ
1 and loses its stability when θ reaches and passes θ
1, which leads to a cascade of the period doubling with θ increasing (see Figure 2).(a) The bifurcation diagram of system (4) with θ varying in (
0.3,0.6
). The system undergoes a flip bifurcation when θ
≈
0.456. (b) With θ increasing, system (4) undergoes period-1 point (red line corresponding to θ
=
0.36), period-2 points (blue dot lines corresponding to θ
=
0.46), period-4 points (green dot lines corresponding to θ
=
0.5), and period-8 points (black dot lines corresponding to θ
=
0.518).
(a)
(b)


## 5. Conclusion

In this paper, a discrete-time prey-predator model with generalized Holling type III functional response and Allee effect has been presented. The stability of fixed points and the existence of flip bifurcation and Neimark-Sacker bifurcation are investigated. Our analysis shows that the strong Allee effect can change topological structure and dynamics of the system. The combined influence of Allee effect and the other parameters induce system (4) to yield more rich dynamical behaviors including chaotic bands and chaotic attractors which do not exist in the continuous case. These results reported in this paper could be very useful for the biologists who devote to study the discrete-time prey-predator models.



---

*Source: 245421-2015-02-25.xml*