# Bifurcation Analysis of a Singular Bioeconomic Model with Allee Effect and Two Time Delays

Authors: Xue Zhang; Qing-ling Zhang; Zhongyi Xiang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A singular prey-predator model with time delays is formulated and analyzed. Allee effect is considered on the growth of the prey population. The singular prey-predator model is transformed into its normal form by using differential-algebraic system theory. We study its dynamics in terms of local analysis and Hopf bifurcation. The existence of periodic solutions via Hopf bifurcation with respect to two delays is established. In particular, we study the direction of Hopf bifurcation and the stability of bifurcated periodic solutions by applying the normal form theory and the center manifold argument. Finally, numerical simulations are included supporting the theoretical analysis and displaying the complex dynamical behavior of the model outside the domain of stability.


---

## Body

## 1. Introduction

The economic theory proposed by Gordon [1] in 1954 was described as follows:net economic revenue (NER) = total revenue (TR) − total cost (TC).This provides theoretical evidence for the establishment of singular bioeconomic model, which is described by differential-algebraic equations. Recently, there has been extensive literature dealing with such systems, describing the interactions between the different species and harvesting effort regarding activity, stabilities of equilibrium, bifurcations, and other dynamics (see, e.g., [2–5] and the references therein).In most of ecosystems, since one species does not respond instantaneously to interactions with other species, some delays due to several reasons, such as gestation, hunting, and maturation, are required. To incorporate this idea in a modelling approach, time delays have been introduced into ecosystems. For a long time, it has been recognized that delays can have very complicated impact on the dynamics of a system (see, e.g., monographs by Hale [6], Hale and Verduyn Lunel [7], Kuang [8], Yuan and Song [9], and Wu [10]). Generally speaking, delays can cause the loss of stability and lead to periodic solutions. Some authors like in [11] described the effects of time delay in a prey-predator model incorporating parasite infection for the prey population. The literature [12] described the dynamics of a stage structured population model with time delay in fluctuating environment. The literature [13] investigated dynamics of delayed prey-predator model with harvesting. Chakraborty et al. [14] introduced a single discrete gestation delay in a differential-algebraic biological economic system and showed Hopf bifurcation in the neighborhood of coexisting equilibrium point through considering the delay as a bifurcation parameter. Zhang et al. [15] studied a ratio-dependent prey-predator singular model and analyzed the direction and stability of periodic solutions. Some other authors have discussed dynamical models with multiple delays. While analyzing, the multiply delayed models are mostly simplified by taking equal magnitude for two delays τ
1
=
τ
2
=
τ [16, 17] or choosing the sum τ of two delays as a bifurcation parameter [18, 19]. However, the delays appearing in different terms of an ecological system are not always equal. Therefore, it is necessary to discuss dynamical system with different delays.In addition, numerous examples demonstrate that the growth of natural populations can exhibit Allee effect [20], which describes a positive relation between population density and the per capita growth rate. When populations get larger, the effect will saturate or disappear. Allee effect may result from some causes, such as mate finding, social dysfunction, inbreeding depression, fool exploitation, and predator avoidance of defense. Therefore, in this paper, we consider a differential-algebraic prey-predator model with two time delays and the Allee effect as follows:
(1)
d
x
d
t
=
x
(
x
A
+
x
-
d
1
-
x
(
t
-
τ
1
)
)
-
x
y
1
+
h
1
x
,
d
y
d
t
=
y
(
b
x
(
t
-
τ
2
)
1
+
h
1
x
(
t
-
τ
2
)
-
d
2
)
-
E
y
,
0
=
E
(
p
y
-
c
)
-
m
,
where A is the Allee effect constant of the prey species, respectively. d
1 and d
2 are the intrinsic mortality rate of the prey and predator population, respectively. h
1 is the handling time. b denotes food utilization efficiency. Prey dynamics is delayed by τ
1 due to crowing and the predator takes time τ
2 to convert the food into its growth; E is harvesting effort for predator, and p
>
0, c
>
0, and m
>
0 are harvesting reward per unit harvesting effort for unit weight of predator, harvesting cost per unit harvesting effort for predator, and the net economic revenue per unit harvesting effort, respectively. All the parameters are positive constants.In this paper, we investigate the effects of two delays on the dynamics of singular bioeconomic model with Allee effect. The existence of periodic solution has been explored through Hopf bifurcation. The domain of stability is defined. Then, the formula for determining the properties of Hopf bifurcation is derived by using the normal form method and center manifold theorem. Finally, numerical simulations show the bifurcation plot with respect to time delay and give the effectiveness of the result mentioned above.


## 2. Stability and Existence of Hopf Bifurcation

There exists a positive equilibriumP
0
=
(
x
0
,
y
0
,
E
0
) for model (1), where coordinates x
0, y
0, and E
0 satisfy the following equations:
(2)
x
A
+
x
-
d
1
-
x
-
y
1
+
h
1
x
=
0
,
b
x
(
t
-
τ
2
)
1
+
h
1
x
(
t
-
τ
2
)
-
d
2
-
E
=
0
,
E
(
p
y
-
c
)
-
m
=
0
.For the sake of the simplicity, denote the two differential equations byf
1
(
x
,
y
,
E
) and f
2
(
x
,
y
,
E
) and the algebraic equation by g
(
x
,
y
,
E
).Letu
=
x, v
=
y, and E
¯
=
(
p
E
0
y
/
(
p
y
0
-
c
)
)
+
E. Therefore, the system (1) can be rewritten in the following form:
(3)
d
u
d
t
=
u
(
u
A
+
u
-
d
1
-
u
(
t
-
τ
1
)
)
-
u
v
1
+
h
1
u
,
d
v
d
t
=
v
(
b
u
(
t
-
τ
2
)
1
+
h
1
u
(
t
-
τ
2
)
-
d
2
-
E
¯
+
p
E
0
v
p
y
0
-
c
)
,
0
=
(
E
¯
-
p
E
0
v
p
y
0
-
c
)
(
p
v
-
c
)
-
m
.
Consider the following local parametric ψ of the third equation of the system (3):
(4)
[
u
,
v
,
E
¯
]
T
=
ψ
(
τ
1
,
τ
2
,
Y
)
=
N
0
T
+
U
0
Y
+
V
0
h
(
Y
(
t
)
,
g
(
ψ
(
Y
(
t
)
)
)
)
=
0
,
where U
0
=
(
1
0
0
1
0
0
), V
0
=
(
0
0
1
), Y
=
(
y
1
,
y
2
)
T, N
0
=
(
u
0
,
v
0
,
E
¯
0
), and h
(
y
1
,
y
2
)
:
R
2
→
R is a smooth mapping and has the following expression:
(5)
h
(
y
1
,
y
2
)
=
p
E
0
(
v
0
+
y
2
)
p
v
0
-
c
+
m
p
(
v
0
+
y
2
)
-
c
-
p
E
0
y
0
p
v
0
-
c
-
E
0
.
Then, the linearization part of the system (3) is
(6)
J
(
N
0
)
=
(
D
N
f
1
(
N
0
)
D
N
f
2
(
N
0
)
)
(
D
N
g
(
N
0
)
U
0
T
)
-
1
(
0
0
1
0
0
1
)
=
(
A
u
0
(
A
+
u
0
)
2
+
h
1
u
0
v
0
(
1
+
h
1
u
0
)
2
-
u
0
e
-
λ
τ
1
-
u
0
1
+
h
1
u
0
b
v
0
(
1
+
h
1
u
0
)
2
e
-
λ
τ
2
p
E
0
v
0
p
v
0
-
c
)
.
For the sake of the simplicity, let
(7)
u
1
=
u
11
+
u
12
,
u
11
=
-
u
0
,
u
12
=
A
u
0
(
A
+
u
0
)
2
+
h
u
0
v
0
(
1
+
h
u
0
)
2
,
u
2
=
-
u
0
1
+
h
u
0
,
u
3
=
b
v
0
(
1
+
h
u
0
)
2
,
u
4
=
p
v
0
E
0
p
v
0
-
m
.According to the Jacobian matrixJ, the characteristic equation of (1) at P
0 can be written as follows:
(8)
R
1
(
λ
)
=
λ
2
-
(
u
4
+
u
12
)
λ
+
u
12
u
4
hhhhhh
-
u
11
(
λ
-
u
4
)
e
-
λ
τ
1
-
u
2
u
3
e
-
λ
τ
2
=
0
.When the time delays are not considered, the characteristic polynomial for model (1) can be expressed as the following form:
(9)
λ
2
-
(
u
1
+
u
4
)
λ
+
u
1
u
4
-
u
2
u
3
=
0
.
It is clear that if the following conditions are satisfied:
(10)
u
1
+
u
4
<
0
,
u
1
u
4
-
u
2
u
3
>
0
,
then P
0 is locally asymptotically stable for model (1).Next, we consider the effects of the two time delays on the stability of model (1) in the following cases.Case 1 (τ
1
=
0 and τ
2
>
0).
Due toτ
1
=
0, (8) becomes
(11)
λ
2
-
(
u
1
+
u
4
)
λ
+
u
1
u
4
-
u
2
u
3
=
0
.Assume that a purely imaginary solution of the formλ
=
i
ω
2 exists in (11). Substituting it into (11) and separating the real and imaginary parts, we have
(12)
-
ω
2
2
+
(
u
11
+
u
12
)
u
4
=
u
2
u
3
cos
⁡
ω
2
τ
2
,
-
(
u
11
+
u
12
+
u
4
)
ω
2
=
-
u
2
u
3
sin
ω
2
τ
2
.
Taking square on both sides of (12) and summing them up, we obtain
(13)
ω
2
4
+
(
u
1
2
+
u
4
2
)
ω
2
2
+
u
1
2
u
4
2
-
u
2
2
u
3
2
=
0
.
From condition (10), (13) has unique positive real root ω
20 if
(14)
u
1
u
4
+
u
2
u
3
<
0
.
The critical value of the delay corresponding to ω
20 is given by
(15)
τ
2
k
=
1
ω
20
arccos
-
ω
20
2
+
u
1
u
4
u
2
u
3
+
2
k
π
ω
20
,
k
=
0,1
,
2
,
…
.
If u
1
u
4
+
u
2
u
3
>
0, (13) has no real root and model (1) is asymptotically stable for any time delay τ
2
>
0.Now, differentiating (8) with respect to τ
2 and substituting the eigenvalue i
ω
20, it follows that
(16)
sign
⁡
(
d
d
τ
2
Re
(
λ
)
)
|
τ
2
k
=
sign
⁡
(
ω
4
-
u
1
2
u
4
2
+
u
2
2
u
3
2
)
>
0
.Theorem 1.
Assume that condition (14) holds; then, there exists a τ
20
>
0 given by
(17)
τ
20
=
1
ω
20
arccos
-
ω
20
2
+
u
1
u
4
u
2
u
3
+
2
k
π
ω
20
,
such that the positive equilibrium P
0 of model (1) is locally asymptotically stable for τ
2
<
τ
20 and undergoes a Hopf bifurcation when τ
2
=
τ
20. That is, model (1) has a branch of periodic solutions bifurcating from the positive equilibrium P
0 near τ
2
=
τ
20.Case 2 (τ
1
>
0 and τ
2
=
0).
According to a similar discussion as in Case1, we can obtain the following result.Theorem 2.
There exists aτ
10
>
0 such that the positive equilibrium P
0 of model (1) is locally asymptotically stable for τ
1
<
τ
10 and undergoes a Hopf bifurcation when τ
1
=
τ
10 given by
(18)
τ
10
=
1
ω
10
arccos
-
u
12
ω
10
2
+
u
4
(
-
u
12
u
4
+
u
2
u
3
)
u
11
(
u
4
2
+
ω
10
2
)
,
where i
ω
10 is the root of corresponding characteristic equation.Case 3 (τ
1
>
0 and τ
2
∈
(
0
,
τ
20
)).
This case states thatτ
1 is regarded as a parameter and τ
2 is in its stable interval. Assume that (8) has purely imaginary solution of the form λ
=
i
ω. Substituting it into (8) and separating the real and imaginary parts, we have
(19)
-
ω
2
+
u
12
u
4
-
u
2
u
3
cos
⁡
ω
τ
2
=
-
u
11
u
4
cos
⁡
ω
τ
1
+
u
11
ω
sin
ω
τ
1
,
-
(
u
12
+
u
4
)
ω
+
u
2
u
3
sin
ω
τ
2
=
u
11
ω
cos
⁡
ω
τ
1
+
u
11
u
4
sin
ω
τ
1
.
Eliminating τ
1 leads to
(20)
ω
4
+
[
(
u
4
+
u
12
)
2
-
2
(
u
12
u
4
-
u
2
u
3
cos
⁡
ω
τ
2
)
-
u
11
2
]
ω
2
-
2
u
2
u
3
(
u
4
+
u
12
)
·
(
sin
ω
τ
2
)
ω
+
u
12
2
u
4
2
+
u
2
2
u
3
2
-
u
11
2
u
4
2
-
2
u
12
u
2
u
3
u
4
cos
⁡
ω
τ
2
=
0
.
It can be seen that there exists at least one real positive root ω
10
′ for (20) if the condition (
u
12
u
4
-
u
2
u
3
)
2
<
u
11
2
u
4
2 holds. Equation (19) can also be written as
(21)
-
u
2
u
3
cos
⁡
ω
τ
2
=
-
u
11
u
4
cos
⁡
ω
τ
1
+
u
11
ω
sin
ω
τ
1
hhhhhhhhhhhhhlh
+
ω
2
-
u
12
u
4
,
u
2
u
3
sin
ω
τ
2
=
u
11
ω
cos
⁡
ω
τ
1
+
u
11
u
4
sin
ω
τ
1
hhhhhhhhhhh
+
(
u
12
+
u
4
)
ω
.
Equation (21) can be solved as
(22)
τ
1
k
′
=
1
ω
10
′
×
(
arccos
(
×
(
ω
2
-
u
12
u
4
)
2
+
ω
2
(
u
4
+
u
12
)
2
)
-
1
(
(
ω
2
-
u
12
u
4
)
2
+
(
u
4
+
u
12
)
2
ω
2
h
h
h
h
h
l
h
h
h
h
h
h
+
u
11
2
u
4
2
+
u
11
2
ω
2
-
u
2
2
u
3
2
(
(
ω
2
-
u
12
u
4
)
2
+
(
u
4
+
u
12
)
2
ω
2
)
h
h
h
h
h
h
j
j
j
h
×
(
-
2
u
11
u
4
2
+
ω
2
h
h
h
h
l
h
h
h
h
h
h
l
h
×
(
ω
2
-
u
12
u
4
)
2
+
ω
2
(
u
4
+
u
12
)
2
)
-
1
)
h
h
h
+
ϕ
1
+
ϕ
2
(
arccos
(
×
(
ω
2
-
u
12
u
4
)
2
+
ω
2
(
u
4
+
u
12
)
2
)
-
1
(
(
ω
2
-
u
12
u
4
)
2
+
(
u
4
+
u
12
)
2
ω
2
)
+
2
k
π
ω
10
′
,
k
=
0,1
,
2
,
…
,
where ϕ
1
=
arctan
(
u
4
/
ω
10
′
) and ϕ
2
=
arctan
⁡
(
(
ω
10
′
2
-
u
12
u
4
)
/
ω
10
′
(
u
4
+
u
12
)
).Now, differentiating (8) with respect to τ
1 and substituting the eigenvalue i
ω
10
′ and time delay τ
1
=
τ
10
′, it follows that
(23)
M
(
d
(
Re
λ
)
d
τ
1
)
|
τ
10
′
+
N
(
d
ω
d
τ
1
)
|
τ
10
′
=
P
,
-
N
(
d
(
Re
λ
)
d
τ
1
)
|
τ
10
′
+
M
(
d
ω
d
τ
1
)
|
τ
10
′
=
Q
,
where
(24)
M
=
-
(
u
4
+
u
12
)
-
u
11
cos
⁡
ω
10
′
τ
10
′
+
τ
10
′
u
11
(
-
u
4
cos
⁡
ω
10
′
τ
10
′
+
ω
10
′
sin
ω
10
′
τ
10
′
)
+
τ
2
u
2
u
3
cos
⁡
ω
10
′
τ
2
,
N
=
-
2
ω
10
′
-
u
11
sin
ω
10
′
τ
10
′
-
τ
10
′
u
11
(
ω
10
′
cos
⁡
ω
10
′
τ
10
′
+
u
4
sin
ω
10
′
τ
10
′
)
+
τ
2
u
2
u
3
sin
ω
10
′
τ
2
,
P
=
u
11
ω
10
′
(
ω
10
′
cos
⁡
ω
10
′
τ
10
′
+
u
4
sin
ω
10
′
τ
10
′
)
,
(25)
Q
=
u
11
ω
10
′
(
-
ω
10
′
sin
ω
10
′
τ
10
′
+
u
4
cos
⁡
ω
10
′
τ
10
′
)
.Solving (23), we obtain
(26)
(
d
(
Re
λ
)
d
τ
1
)
|
τ
10
′
=
M
P
-
N
Q
M
2
+
N
2
≠
0
.Theorem 3.
Assume(
u
12
u
4
-
u
2
u
3
)
2
<
u
11
2
u
4
2 and τ
2
∈
[
0
,
τ
20
). The positive equilibrium P
0 of model (1) is asymptotically stable for τ
1
∈
[
0
,
τ
10
′
) and undergoes Hopf bifurcation at τ
1
=
τ
10
′, where
(27)
τ
10
′
=
1
ω
10
′
h
l
×
(
arccos
(
(
ω
10
′
2
-
u
12
u
4
)
2
+
ω
10
′
2
(
u
4
+
u
12
)
2
)
-
1
(
(
ω
10
′
2
-
u
12
u
4
)
2
+
(
u
4
+
u
12
)
2
ω
10
′
2
h
h
h
h
h
h
h
h
h
h
h
×
h
+
u
11
2
u
4
2
+
u
11
2
ω
10
′
2
-
u
2
2
u
3
2
(
(
ω
10
′
2
-
u
12
u
4
)
2
+
(
u
4
+
u
12
)
2
ω
10
′
2
)
h
h
h
h
h
h
h
h
h
h
h
h
×
(
-
2
u
11
u
4
2
+
ω
10
′
2
h
h
h
h
h
h
h
h
h
h
l
h
h
(
-
×
(
ω
10
′
2
-
u
12
u
4
)
2
+
ω
10
′
2
(
u
4
+
u
12
)
2
)
-
1
)
lhlllhh
+
ϕ
1
+
ϕ
2
(
arccos
(
(
ω
10
′
2
-
u
12
u
4
)
2
+
ω
10
′
2
(
u
4
+
u
12
)
2
)
-
1
(
(
ω
10
′
2
-
u
12
u
4
)
2
+
(
u
4
+
u
12
)
2
ω
10
′
2
)
,
where ϕ
1
=
arctan
(
u
4
/
ω
10
′
) and ϕ
2
=
arctan
(
(
ω
10
′
2
-
u
12
u
4
)
/
ω
10
′
(
u
4
+
u
12
)
).Case 4 (τ
1
∈
(
0
,
τ
10
) and τ
2
>
0).
According to a similar discussion as in Case3, we can obtain the following result.Theorem 4.
The positive equilibriumP
0 of model (1) is asymptotically stable for τ
2
∈
[
0
,
τ
20
′
) and undergoes Hopf bifurcation at τ
2
=
τ
20
′, where
(28)
τ
20
′
=
1
ω
20
′
×
(
h
h
h
h
h
h
h
l
×
(
u
4
+
u
12
)
2
ω
20
′
2
+
(
ω
20
′
2
-
u
12
u
4
)
2
)
-
1
arccos
(
h
h
h
h
h
h
h
l
h
h
h
h
h
h
l
×
(
u
4
+
u
12
)
2
ω
20
′
2
+
(
ω
20
′
2
-
u
12
u
4
)
2
h
h
h
h
h
h
h
l
×
(
u
4
+
u
12
)
2
ω
20
′
2
+
(
ω
20
′
2
-
u
12
u
4
)
2
)
-
1
)
-
1
(
(
ω
20
′
2
-
u
12
u
4
)
2
+
(
u
4
+
u
12
)
2
ω
20
′
2
h
h
h
h
h
h
h
h
h
h
h
h
l
h
h
+
u
12
2
u
3
2
-
u
11
2
u
4
2
-
u
11
2
ω
20
′
2
(
(
ω
20
′
2
-
u
12
u
4
)
2
+
(
u
4
+
u
12
)
2
ω
20
′
2
)
h
h
h
h
h
h
h
h
h
h
l
h
×
(
h
h
h
h
h
h
h
l
×
(
u
4
+
u
12
)
2
ω
20
′
2
+
(
ω
20
′
2
-
u
12
u
4
)
2
)
-
1
-
2
u
2
u
3
h
h
h
h
h
h
h
l
h
h
h
h
h
h
l
×
(
u
4
+
u
12
)
2
ω
20
′
2
+
(
ω
20
′
2
-
u
12
u
4
)
2
h
h
h
h
h
h
h
l
×
(
u
4
+
u
12
)
2
ω
20
′
2
+
(
ω
20
′
2
-
u
12
u
4
)
2
)
-
1
)
-
1

)
h
h
h
h
h
h
h
-
ψ
h
h
h
h
h
h
h
l
×
(
u
4
+
u
12
)
2
ω
20
′
2
+
(
ω
20
′
2
-
u
12
u
4
)
2
)
-
1
)
,
where ψ
=
arctan
(
(
u
4
+
u
12
)
ω
20
′
/
(
ω
20
′
2
-
u
12
u
4
)
).


## 3. Direction and the Stability of Hopf Bifurcation

In the previous section, we obtain the conditions under which a family of periodic solutions bifurcate from the positive equilibrium at the critical values of time delaysτ
1 and τ
2. As pointed out by Hassard et al., by employing the normal form and center manifold theory, the formulae for determining the directions, stability, and period of Hopf bifurcation can be presented. Following the ideas of Hassard et al., this section discusses the directions, stability, and period of Hopf bifurcation. These properties are studied with respect to τ
1 for fixed τ
2
∈
(
0
,
τ
20
). Throughout this section, it is considered that the system (1) undergoes Hopf bifurcation at τ
1
=
τ
10
′, τ
2
*
∈
(
0
,
τ
20
) at the equilibrium P
0. Without loss of generality, this section assumes that τ
2
*
<
τ
10
′.Letx
1
=
y
1
-
y
1
*, x
2
=
y
2
-
y
2
*, t
=
t
/
τ
1, x
¯
i
(
t
)
=
x
i
(
τ
1
t
), and τ
1
=
τ
10
′
+
μ and dropping the bars for simplification of notations, the system (1) can be written as a functional differential equation in C
=
C
(
[
-
1,0
]
,
R
2
)
(29)
x
˙
(
t
)
=
L
μ
(
x
t
)
+
F
(
μ
,
x
t
)
,
where x
(
t
)
=
(
x
1
(
t
)
,
x
2
(
t
)
)
T
∈
R
2, and L
μ
:
C
→
R, F
:
R
×
C
→
R are given, respectively, by
(30)
L
μ
(
ϕ
)
=
(
τ
10
′
+
μ
)
h
i
×
(
(
A
u
0
(
A
+
u
0
)
2
+
h
1
u
0
v
0
(
1
+
h
1
u
0
)
2
-
u
0
1
+
h
1
u
0
0
p
E
0
v
0
p
v
0
-
c
)
h
h
h
h
h
h
i
×
ϕ
(
0
)
+
(
0
0
b
v
0
(
1
+
h
1
u
0
)
2
0
)
ϕ
(
-
τ
2
*
τ
10
′
)
h
h
h
h
h
i
h
+
(
-
u
0
0
0
0
)
ϕ
(
-
1
)
(
(
A
u
0
(
A
+
u
0
)
2
+
h
1
u
0
v
0
(
1
+
h
1
u
0
)
2
-
u
0
1
+
h
1
u
0
0
p
E
0
v
0
p
v
0
-
c
)
)
,
F
(
μ
,
ϕ
)
=
(
τ
10
′
+
μ
)
×
(
(
A
2
(
A
+
u
0
)
3
+
h
1
v
0
(
1
+
h
1
u
0
)
3
)
ϕ
1
2
(
0
)
hhh
-
ϕ
1
(
0
)
ϕ
1
(
-
1
)
-
1
(
1
+
h
1
u
0
)
2
ϕ
1
(
0
)
ϕ
2
(
0
)
hhh
-
(
A
2
(
A
+
u
0
)
4
+
h
1
2
v
0
(
1
+
h
1
u
0
)
4
)
·
ϕ
1
3
(
0
)
hhh
+
h
1
(
1
+
h
1
u
0
)
3
ϕ
1
2
(
0
)
ϕ
2
(
0
)
+
⋯
hhh
-
b
h
1
v
0
(
1
+
h
1
u
0
)
3
ϕ
1
2
(
-
τ
2
*
τ
10
′
)
hhh
+
b
(
1
+
h
1
u
0
)
2
ϕ
1
(
-
τ
2
*
τ
10
′
)
ϕ
2
(
0
)
hhh
-
p
E
0
c
(
p
v
0
-
c
)
2
ϕ
2
2
(
0
)
+
b
v
0
h
1
2
(
1
+
h
1
u
0
)
4
ϕ
1
3
(
-
τ
2
*
τ
10
′
)
hhh
-
b
h
1
(
1
+
h
1
u
0
)
3
ϕ
1
2
(
-
τ
2
*
τ
10
′
)
ϕ
2
(
0
)
h
h
h
+
m
p
2
c
p
v
0
-
c
ϕ
2
3
(
0
)
+
⋯
)
,
where ϕ
(
θ
)
=
(
ϕ
1
(
θ
)
,
ϕ
2
(
θ
)
)
T
∈
C. By the Riesz representation theorem, there exists a 2
×
2 matrix η
(
θ
,
μ
) whose components are bounded variation for θ
∈
[
-
1,0
], such that
(31)
L
μ
ϕ
=
∫
-
1
0
‍
d
η
(
θ
,
μ
)
ϕ
(
θ
)
,
for
ϕ
∈
C
,
where we choose
(32)
η
(
θ
,
μ
)
=
{
(
τ
10
′
+
μ
)
×
(
A
u
0
(
A
+
u
0
)
2
+
h
1
u
0
v
0
(
1
+
h
1
u
0
)
2
-
u
0
-
u
0
1
+
h
1
u
0
b
v
0
(
1
+
h
1
u
0
)
2
p
E
0
v
0
p
v
0
-
c
)
,
ccccccccccccccccccccccccccccccchcccccccccc
θ
=
0
,
(
τ
10
′
+
μ
)
(
-
u
0
0
b
v
0
(
1
+
h
1
u
0
)
2
0
)
,
θ
∈
[
-
τ
2
*
τ
10
′
,
0
)
,
(
τ
10
′
+
μ
)
(
-
u
0
0
0
0
)
,
θ
∈
(
-
1
,
-
τ
2
*
τ
10
′
)
,
0
2
×
2
,
θ
=
-
1
.
For ϕ
∈
C
1
(
[
-
1,0
]
,
R
2
), define
(33)
A
(
μ
)
ϕ
=
{
d
ϕ
(
θ
)
d
θ
,
θ
∈
[
-
1,0
)
,
∫
-
1
0
‍
d
η
(
s
,
μ
)
ϕ
(
s
)
,
θ
=
0
,
R
(
μ
)
ϕ
=
{
0
,
θ
∈
[
-
1,0
)
,
F
(
μ
,
ϕ
)
,
θ
=
0
.
Then, system (29) is equivalent to
(34)
x
˙
t
=
A
(
μ
)
x
t
+
R
(
μ
)
x
t
,
where x
t
(
θ
)
=
x
(
t
+
θ
) for θ
∈
[
-
1,0
].Forψ
∈
C
1
(
[
0,1
]
,
(
R
2
)
*
), define
(35)
A
*
ψ
(
s
)
=
{
-
d
ψ
(
s
)
d
s
,
s
∈
(
0,1
]
,
∫
-
1
0
‍
d
η
T
(
t
,
0
)
ψ
(
-
t
)
,
s
=
0
,
and bilinear form
(36)
〈
ψ
(
s
)
,
ϕ
(
θ
)
〉
=
ψ
¯
(
0
)
ϕ
(
0
)
-
∫
-
1
0
‍
∫
ξ
=
0
θ
‍
ψ
¯
(
ξ
-
θ
)
d
η
(
θ
)
ϕ
(
ξ
)
d
ξ
,
where η
(
θ
)
=
η
(
θ
,
0
), A
=
A
(
0
), and A
* are adjoint operators. Since ±
i
ω
10
′
τ
10
′ are eigenvalues of A
(
0
), they are also eigenvalues of A
*. Next, we need to compute the eigenvector of A
(
0
) and A
* corresponding to i
ω
10
′
τ
10
′ and -
i
ω
10
′
τ
10
′, respectively.In terms of the discussion mentioned above, we see that±
i
ω
10
′
τ
10
′ are eigenvalues of A
(
0
) and A
*. Next, we calculate the eigenvector q
(
θ
) of A
(
0
) corresponding to i
ω
10
′
τ
10
′ and the eigenvector q
*
(
s
) of A
* belonging to -
i
ω
10
′
τ
10
′, respectively. By the definition of A
(
0
) and A
*, it is not difficult to obtain that
(37)
α
=
(
i
ω
10
′
+
u
0
e
-
i
ω
10
′
τ
10
′
-
A
u
0
(
A
+
u
0
)
2
-
h
1
u
0
v
0
(
1
+
h
1
u
0
)
2
)
×
(
-
u
0
(
1
+
h
1
u
0
)
)
-
1
.
Similarly, let q
*
(
s
)
=
D
(
1
,
α
*
)
e
i
ω
10
′
τ
10
s
′ be the eigenvector of A
* corresponding to -
i
ω
10
′
τ
10
′, and we can obtain
(38)
α
*
=
(
-
i
ω
10
′
+
u
0
e
i
ω
10
′
τ
10
′
-
A
u
0
(
A
+
u
0
)
2
-
h
1
u
0
v
0
(
1
+
h
1
u
0
)
2
)
×
(
b
v
0
(
1
+
h
1
u
0
)
2
e
i
ω
10
′
τ
2
*
)
-
1
.
It follows from (36) that
(39)
〈
q
*
(
s
)
,
q
(
θ
)
〉
=
L
¯
(
1
,
α
¯
*
)
(
1
,
α
)
T
-
∫
-
1
0
‍
∫
ξ
=
0
θ
‍
L
¯
(
1
,
α
¯
*
)
e
-
i
ω
10
′
τ
10
′
(
ξ
-
θ
)
d
η
(
θ
)
(
1
,
α
)
T
e
i
ω
10
′
τ
10
′
ξ
d
ξ
=
L
¯
(
1
+
α
¯
*
α
-
u
0
e
-
i
ω
10
′
τ
10
′
+
b
v
0
τ
2
*
α
¯
*
τ
10
′
(
1
+
h
1
u
0
)
2
e
-
i
ω
10
′
τ
2
*
)
.
Then, we obtain
(40)
L
¯
=
(
1
+
α
¯
*
α
-
u
0
e
-
i
ω
10
′
τ
10
′
+
b
v
0
τ
2
*
α
¯
*
τ
10
′
(
1
+
h
1
u
0
)
2
e
-
i
ω
10
′
τ
2
*
)
-
1
,
such that 〈
q
*
(
s
)
,
q
(
θ
)
〉
=
1, 〈
q
*
(
s
)
,
q
¯
(
θ
)
〉
=
0.Based on the ideas in Hassard et al. [21], we can obtain g
(
z
,
z
¯
)
=
q
¯
*
(
0
)
F
0
(
z
,
z
¯
) and the following expressions:
(41)
x
1
t
(
0
)
=
W
(
1
)
(
t
,
0
)
+
z
+
z
¯
,
x
2
t
(
0
)
=
W
(
2
)
(
t
,
0
)
+
α
z
+
α
¯
z
¯
,
x
1
t
(
-
τ
2
*
τ
10
′
)
=
W
(
1
)
(
t
,
-
τ
2
*
τ
10
′
)
+
z
e
-
i
ω
10
′
τ
2
*
+
z
¯
e
i
ω
10
′
τ
2
*
,
x
2
t
(
-
τ
2
*
τ
10
′
)
=
W
(
2
)
(
t
,
-
τ
2
*
τ
10
′
)
+
α
z
e
-
i
ω
10
′
τ
2
*
+
α
¯
z
¯
e
i
ω
10
′
τ
2
*
,
x
1
t
(
-
1
)
=
W
(
1
)
(
t
,
-
1
)
+
z
e
-
i
ω
10
′
τ
10
′
+
z
¯
e
i
ω
10
′
τ
10
′
,
x
2
t
(
-
1
)
=
W
(
2
)
(
t
,
-
1
)
+
α
z
e
-
i
ω
10
′
τ
10
′
+
α
¯
z
¯
e
i
ω
10
′
τ
10
′
,
where those important parameters can be computed as follows:
(42)
g
20
=
2
D
¯
τ
10
′
[
(
A
2
(
A
+
u
0
)
3
+
h
1
v
0
(
1
+
h
1
u
0
)
3
)
hhhhhhhh
-
e
-
i
ω
10
′
τ
10
′
-
α
(
1
+
h
1
u
0
)
2
hhhhhhhh
+
α
¯
*
(
b
α
(
1
+
h
1
u
0
)
2
e
-
i
ω
10
′
τ
2
*
hhhhhhhhhhhhhh
-
b
h
1
v
0
(
1
+
h
1
u
0
)
3
e
-
2
i
ω
10
′
τ
2
*
hhhhhhhhhhhhhh
-
p
E
0
c
α
2
(
p
v
0
-
c
)
2
)
]
,
g
11
=
2
D
¯
τ
10
′
[
(
A
2
(
A
+
u
0
)
3
+
h
1
v
0
(
1
+
h
1
u
0
)
3
)
hhhhhhhh
-
Re
{
e
i
ω
10
′
τ
10
′
}
-
Re
{
α
}
(
1
+
h
1
u
0
)
2
hhhhhhhh
+
α
¯
*
(
b
(
1
+
h
1
u
0
)
2
Re
{
α
e
i
ω
10
′
τ
2
*
}
hhhhhhhhhhhhhh
-
b
h
1
v
0
(
1
+
h
1
u
0
)
3
-
p
E
0
c
α
α
¯
(
p
v
0
-
c
)
2
)
]
,
g
02
=
2
D
¯
τ
10
′
[
(
A
2
(
A
+
u
0
)
3
+
h
1
v
0
(
1
+
h
1
u
0
)
3
)
hhhhhhhh
-
e
i
ω
10
′
τ
10
′
-
α
¯
(
1
+
h
1
u
0
)
2
hhhhhhhh
+
α
¯
*
(
b
α
¯
(
1
+
h
1
u
0
)
2
e
i
ω
10
′
τ
2
*
hhhhhhhhhhhhhh
-
b
h
1
v
0
(
1
+
h
1
u
0
)
3
e
2
i
ω
10
′
τ
2
*
hlhhhhhhhhhhhlh
-
p
E
0
c
α
¯
2
(
p
v
0
-
c
)
2
)
]
,
g
21
=
2
D
¯
τ
10
′
{
(
A
2
(
A
+
u
0
)
3
+
h
1
v
0
(
1
+
h
1
u
0
)
3
)
h
h
h
h
h
h
h
h
×
(
2
W
11
(
1
)
(
0
)
+
W
20
(
1
)
(
0
)
)
h
h
h
h
h
h
h
h
-
(
W
11
(
1
)
(
-
1
)
+
W
20
(
1
)
(
-
1
)
2
+
W
20
1
(
0
)
2
e
i
ω
10
′
τ
10
′
h
h
h
h
h
h
h
h
h
h
h
h
+
W
11
(
1
)
(
0
)
e
-
i
ω
10
′
τ
10
′
W
20
(
1
)
(
-
1
)
2
)
-
1
(
1
+
h
1
u
0
)
2
hhhhhhjjhl
×
(
W
11
2
(
0
)
+
W
20
(
2
)
(
0
)
2
h
h
h
h
h
h
h
l
h
h
h
h
h
+
W
20
(
1
)
(
0
)
2
α
¯
+
W
11
(
1
)
(
0
)
α
)
h
h
h
h
h
h
h
h
-
3
(
A
2
(
A
+
u
0
)
4
+
v
0
h
1
2
(
1
+
h
1
u
0
)
4
)
h
h
h
h
h
h
h
h
+
h
1
(
1
+
h
1
u
0
)
3
(
α
¯
+
2
α
)
jjjjj
h
h
h
h
h
h
+
α
¯
*
(
b
(
1
+
h
1
u
0
)
2
h
h
h
h
h
h
h
j
j
j
j
h
h
h
h
×
(
(
-
τ
2
*
τ
10
′
)
α
W
11
(
2
)
(
0
)
e
-
i
ω
10
′
τ
2
*
hhhhhhhhjjhhhhhhlh
+
W
20
(
2
)
(
0
)
2
e
i
ω
10
′
τ
2
*
hhhhhhhhjjhhhhhhlh
+
W
20
(
1
)
(
-
τ
2
*
/
τ
10
′
)
2
α
¯
hhhhhhhjjhhhhhhlllh
+
W
11
(
1
)
(
-
τ
2
*
τ
10
′
)
α
)
-
b
h
1
v
0
(
1
+
h
1
u
0
)
3
hhhhhhhhhhhlhlh
×
(
2
W
11
(
1
)
(
-
τ
2
*
τ
10
′
)
e
-
i
ω
10
′
τ
2
*
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
+
W
20
(
1
)
(
-
τ
2
*
τ
10
′
)
e
i
ω
10
′
τ
2
*
)
hhhhhhhhhhlhhlh
-
p
E
0
c
(
p
v
0
-
c
)
2
(
2
W
11
(
2
)
(
0
)
α
+
W
20
(
2
)
(
0
)
α
¯
)
hhhhhhhhhhhhlh
+
3
b
v
0
h
1
2
(
1
+
h
1
u
0
)
4
e
-
i
ω
10
′
τ
2
*
-
b
h
1
(
1
+
h
1
u
0
)
3
hhhhhhhhhhhhlh
×
(
α
¯
e
-
2
i
ω
10
′
τ
2
*
+
2
α
)
+
3
m
p
2
c
α
2
α
¯
(
p
v
0
-
c
)
4
)
}
.
According to the computation process similar to Hassard’s, parameters of g
21 can be written as follows:
(43)
S
1
=
(
A
2
(
A
+
u
0
)
3
+
h
1
v
0
(
1
+
h
1
u
0
)
3
-
e
-
i
ω
10
′
τ
10
′
-
α
(
1
+
h
1
u
0
)
2
b
α
(
1
+
h
1
u
0
)
2
e
-
i
ω
10
′
τ
2
*
-
b
h
1
v
0
(
1
+
h
1
u
0
)
3
e
-
2
i
ω
10
′
τ
2
*
-
p
E
0
c
α
2
(
p
v
0
-
c
)
2
)
,
S
2
=
(
A
2
(
A
+
u
0
)
3
+
h
1
v
0
(
1
+
h
1
u
0
)
3
-
Re
(
e
i
ω
10
′
τ
10
′
)
-
Re
(
α
)
(
1
+
h
1
u
0
)
2
b
(
1
+
h
1
u
0
)
2
Re
(
α
e
i
ω
10
τ
2
*
)
-
b
h
1
v
0
(
1
+
h
1
u
0
)
3
-
p
E
0
c
α
α
¯
(
p
v
0
-
c
)
2
)
,
R
1
=
2
(
2
i
ω
10
′
-
A
u
0
(
A
+
u
0
)
2
-
h
1
u
0
v
0
(
1
+
h
1
u
0
)
2
+
u
0
e
-
2
i
ω
10
′
τ
10
′
u
0
1
+
h
1
u
0
-
b
v
0
1
+
h
1
u
0
e
-
2
i
ω
10
′
τ
2
*
2
i
ω
10
′
-
p
E
0
v
0
p
v
0
-
c
)
-
1


S
1
,
R
2
=
2
(
u
0
-
A
u
0
(
A
+
u
0
)
2
-
h
1
u
0
v
0
(
1
+
h
1
u
0
)
2
u
0
1
+
h
1
u
0
-
b
v
0
(
1
+
h
1
u
0
)
2
-
p
E
0
v
0
p
v
0
-
c
)
-
1
S
2
,
W
20
(
θ
)
=
i
g
20
ω
10
′
τ
10
′
q
(
0
)
e
i
ω
10
′
τ
10
′
θ
+
i
g
¯
02
3
ω
10
′
τ
10
′
q
¯
(
0
)
e
-
i
ω
10
′
τ
10
′
θ
+
R
1
e
2
i
ω
10
′
τ
10
′
θ
,
W
11
(
θ
)
=
-
i
g
11
ω
10
′
τ
10
′
q
(
0
)
e
i
ω
10
′
τ
10
′
θ
+
i
g
¯
11
ω
10
′
τ
10
′
q
¯
(
0
)
e
-
i
ω
10
′
τ
10
′
θ
+
R
2
.Therefore, we can computeg
21 by the parameters and delays. Then, the following results can be expressed as
(44)
c
1
(
0
)
=
i
2
ω
10
′
τ
10
′
(
g
20
g
11
-
2
|
g
11
|
2
-
1
3
|
g
02
|
2
)
+
g
21
2
,
μ
2
=
-
Re
{
c
1
(
0
)
}
Re
{
λ
′
(
τ
10
′
)
}
,
β
2
=
2
Re
{
c
1
(
0
)
}
,
T
2
=
-
Im
⁡
{
c
1
(
0
)
}
+
μ
2
Im
⁡
{
λ
′
(
τ
10
′
)
}
ω
10
′
τ
10
′
.
These parameters give a description of the bifurcating periodic solutions in the center manifold of system (1) at critical values τ
1
=
τ
10
′, which can be expressed as follows.(i)
μ
2 determines the direction of the Hopf bifurcation: if μ
2
>
0
(
μ
2
<
0
), the Hopf bifurcation is supercritical (subcritical).
(ii)
β
2 gives the stability of periodic solution: if β
2
<
0
(
β
2
>
0
), the periodic solution is stable (unstable).
(iii)
T
2 expresses the period of the bifurcating periodic solutions: if T
2
>
0
(
T
2
<
0
), the period of the bifurcating periodic solution increases (decreases).


## 4. Numerical Simulation

With the help of MATLAB, we present some numerical simulations to verify and extend our theoretical analysis proved in Sections2 and 3. The parameters of model (1) are chosen as A
=
0.02, d
1
=
0.01, h
1
=
1.25, b
=
7.5, d
2
=
2, p
=
35, c
=
1, and m
=
0.5.According to the parameters given above, we can obtain the positive equilibrium pointP
0
(
0.4055,0.8099,0.0183
) of model (1). It is easy to obtain that P
0 is asymptotically stable in the absence of delay (see Figure 1). Moreover, since u
1
u
4
+
u
2
u
3
=
-
0.7233
<
0, there is a positive root ω
20
=
0.8389 for (13). Hence, when τ
1
=
0, the critical value of the delay corresponding to ω
20 is τ
20
=
0.2247, which is shown in Figure 2. When τ
2
=
0, by a similar computation, we can obtain that the critical delay for τ
1 is τ
10
=
1.272.Without any time delay, the positive equilibrium pointP
0 is asymptotically stable.
(a)
(b)
(c)
(d)Figure 2
Whenτ
2
=
0.288
>
τ
20, bifurcating periodic solutions from the equilibrium point P
0 occur.Based on the results of Cases3 and 4 in Section 2, we draw a critical curve τ
2
=
f
(
τ
1
) with respect to two parameters τ
1 and τ
2 (see Figure 3). The main algorithm for drawing the critical curve is searching the values τ
10
′(τ
20
′) for fixed positive values of τ
2
∈
(
0
,
τ
20
)  (
τ
1
∈
(
0
,
τ
10
)
). From Figure 3, it is clear that the domain surrounded by the critical curve f
(
τ
1
) and the two axes, that is, the gray region, is stable for model (1).Figure 3
The critical curve of stable domain with respect toτ
1 and τ
2.In addition, according to the algorithm derived in Section3, we can obtain the following values: c
1
(
0
)
=
-
1.9243
-
3.0887
i, ω
10
′
=
0.7538, τ
10
′
=
0.6546, λ
′
(
τ
10
′
)
=
0.0879
+
0.1826
i, μ
2
=
21.8851, β
2
=
-
3.8487, and T
2
=
-
1.8402. By the discussion in Section 3, we know that model (1) can undergo a supercritical Hopf bifurcation at the positive equilibrium P
0 and the bifurcating periodic solution occurs when τ
1 crosses τ
10
′ to the right with fixed τ
2
*
∈
(
0
,
τ
20
) and the bifurcating periodic solution is stable. Figure 4 plots the bifurcating periodic solution of the model (1) at τ
1
=
τ
10
′
=
0.6546 and τ
2
*
=
0.15.When atτ
1
=
τ
10
′
=
0.6546 and τ
2
*
=
0.15, bifurcating periodic solutions of model (1) from P
0 occur.
(a)
(b)
(c)
(d)To explore the possibility of occurrence of chaos, bifurcation diagrams are plotted for the key parameterτ
1. Bifurcation diagram in Figure 5 of model (1) in (τ
1
-
x) plane is given for the fixed τ
2
=
0.15. We can see that model (1) experiences the processes of periodic, periodic doubling cascade and chaos in the region 0.98
<
τ
1
<
1.018. When 0.98
<
τ
1
<
0.996, model (1) exhibits a stable limit cycle. When τ
1 is increasing, the limit cycle is unstable and there is a cascade of periodic doubling bifurcation (see Figure 6) leading to chaos.Figure 5
Bifurcation diagram of model (1) showing the effect of τ
1 on the dynamical behavior.The projection of the attracting set of model (1) in the x
-
y plane for (a) τ
1
=
1.01, (b) τ
1
=
1.0125, and (c) τ
1
=
1.017, respectively.
(a)
(b)
(c)


## 5. Conclusions

In this paper, dynamical complexity of a multiple-delayed predator-prey model with Allee effect is analyzed by using differential-algebraic system theory. The prey-predator model with single delay has been investigated by many researchers (see, e.g., [22–25]). However, to our best knowledge, there are few papers on the bifurcations of differential-algebraic population dynamics with two or multiple delays. On the other hand, it should be noted that almost the existing bioeconomic models (see [2–5, 14, 15]) only investigate the simplest case of a logistic prey growth function. Compared with these works, the introduction of Allee effect makes the work studied in this paper novel.This paper provides a new and efficient method for the qualitative analysis of the Hopf bifurcation of differential-algebraic biological economic system with multiple delays. From the analysis of the proposed model, we have obtained some interesting and useful results. The main contribution of this paper is that the stability domain with respect to two delays is defined and plotted. By choosing the time delayτ
1 as a bifurcation parameter, the direction of Hopf bifurcation and stability of the bifurcating periodic orbit are discussed by using normal form and center manifold. Moreover, numerical simulations are carried out which substantiate the theoretical results. At the same time, bifurcation diagrams and attracting sets are plotted, which validate the existence of chaos due to periodic doubling bifurcation. However, those theoretical analyses, such as the existence of chaos and conditions of periodic doubling bifurcation, will be the topics of the future research.



---

*Source: 745296-2014-03-10.xml*