# Dynamical Behavior and Stability Analysis in a Stage-Structured Prey Predator Model with Discrete Delay and Distributed Delay

Authors: Chao Liu; Qingling Zhang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We propose a prey predator model with stage structure for prey. A discrete delay and a distributed delay for predator described by an integral with a strong delay kernel are also considered. Existence of two feasible boundary equilibria and a unique interior equilibrium are analytically investigated. By analyzing associated characteristic equation, local stability analysis of boundary equilibrium and interior equilibrium is discussed, respectively. It reveals that interior equilibrium is locally stable when discrete delay is less than a critical value. According to Hopf bifurcation theorem for functional differential equations, it can be found that model undergoes Hopf bifurcation around the interior equilibrium when local stability switch occurs and corresponding stable limit cycle is observed. Furthermore, directions of Hopf bifurcation and stability of the bifurcating periodic solutions are studied based on normal form theory and center manifold theorem. Numerical simulations are carried out to show consistency with theoretical analysis.


---

## Body

## 1. Introduction

In recent years, much research efforts have been extensively made on interaction and coexistence mechanism of population in prey predator ecosystem by means of Lotka-Volterra dynamical models [1–3]. Generally, in order to reflect the dynamical behavior of mathematical models depending on the past history, time delay is usually incorporated into model, which can be utilized to mathematically describe hunting delay, maturation delay, and gestation delay for population within prey predator ecosystem. With the introduction of time delay, it may cause the loss of stability and other complicated dynamical behavior of model, such as Hopf bifurcation and saddle node bifurcation. It should be noted that there is a well-developed theory of dynamical models which incorporates time delay into model [4]. In particular, the properties of periodic solutions arising from the Hopf bifurcation are of great interest [4, 5].Recently, plenty of dynamical models with discrete and distributed delays have been proposed to discuss the population dynamics of prey predator ecosystem due to variation of maturation and hunting factors [6–13]. Song and Peng [9] proposed a logistic model with discrete and distributed delays
(1)
x
˙
(
t
)
=
r
x
(
t
)
[
1
-
a
1
x
(
t
-
τ
)
-
a
2
∫
-
∞
t
f
(
t
-
s
)
x
(
s
)
d
s
]
,
where r
,
τ
,
a
1, and a
2 are positive constants. Function f is called the delayed kernel, which is the weight given to the population at time t. Under the assumption that f
(
t
)
≥
0 for all t
≥
0 and the normalized condition that ∫
0
∞
f
(
t
)
d
t
=
1, which ensures that the steady state of model (1) is unaffected by the delay, the local stability of model (1) around interior equilibrium and existence of Hopf bifurcation are studied. Furthermore, direction of Hopf bifurcations and the stability of bifurcated periodic solutions are investigated by using the normal form theory and center manifold theorem.By supposing that the predator population at every age stage has the predation ability and the prey population captured by the predator population in the past is all contributing to the predator population at timet, the growth dynamics of the two species can be described by the following delayed Lotka-Volterra two species prey predator model with distributed delays:
(2)
x
˙
(
t
)
=
x
(
t
)
[
r
1
-
a
11
x
(
t
)
-
a
12
∫
-
∞
t
F
(
t
-
s
)
y
(
s
)
d
s
]
,
y
˙
(
t
)
=
y
(
t
)
[
-
r
2
+
a
21
∫
-
∞
t
G
(
t
-
s
)
x
(
s
)
d
s
-
a
22
y
(
s
)
]
,
where constants r
1
>
0 and r
2
>
0 denote intrinsic growth rate of prey population and death rate of predator population, respectively. r
1
/
a
11, a
12, a
21 and a
22 represent the carrying capacity of prey population, predator coefficient of predator population, biomass conversion rate of prey population captured by predator population, and intraspecific competition rate of predator population, respectively. It should be noted that a
i
j
>
0 for i, j
=
1, 2, the delay kernel F
(
s
) and G
(
s
) are bounded nonnegative functions and the following normalized conditions hold:
(3)
∫
0
∞
F
(
s
)
d
s
=
1
,
∫
0
∞
G
(
s
)
d
s
=
1
.Model system (2) with various delay kernels and delayed intraspecific competitions have been extensively investigated in [6–13]. Faria [6] investigated the stability of interior equilibrium of model (2) and Hopf bifurcation of nonconstant periodic solutions around the interior equilibrium. When F
(
s
)
=
δ
(
s
-
τ
) and G
(
s
)
=
δ
(
s
-
η
), τ, η
≥
0, model (2) admits two different discrete time delays; Ruan [7] and Yan and Zhang [11] discussed the stability of the interior equilibrium of model (2) and Hopf bifurcation of nonconstant periodic solutions regarding the sum of two delays τ and η as the bifurcation parameter. Furthermore, dynamic effect of intraspecific competition on population dynamics of model (2) is studied in [8, 9]. By assuming that F
(
s
)
=
δ
(
s
-
τ
), τ
≥
0, the delay kernel function G
(
s
) may take the weak generic kernel function G
(
s
)
=
α
e
-
α
s and strong generic kernel function G
(
s
)
=
α
2
s
e
-
α
s(α
>
0), where the weak generic kernel implies that the importance of events in the past simply decreases exponentially and the further one looks into the past while the strong generic kernel implies that a particular time in the past is more important than any other [1, 2, 5]. Under this assumption, model (2) can be reduced to the following system with a discrete delay and a distributed delay:
(4)
x
˙
(
t
)
=
x
(
t
)
[
r
1
-
a
11
x
(
t
)
-
a
12
y
(
t
-
τ
)
]
,
y
˙
(
t
)
=
y
(
t
)
[
-
r
2
+
a
21
∫
-
∞
t
G
(
t
-
s
)
x
(
s
)
d
s
-
a
22
y
(
t
)
]
.When the delay kernel functionG
(
s
) admits the weak generic kernel, Song and Yuan [10] and Ma et al. [12] discussed the local asymptotical stability of interior equilibrium and Hopf bifurcations of nonconstant periodic solutions based on linearization method and regarding the discrete hunting delay τ as bifurcation parameter. It is shown that interior equilibrium is asymptotically stable when τ is less than a certain critical value and Hopf bifurcation occurs at a critical value of the discrete hunting delay. They also studied the direction of the Hopf bifurcations and the stability of bifurcated periodic solutions occurring through Hopf bifurcations. In the case of strong generic kernel, Zhang et al. [13] studied local stability of all boundary and interior equilibria; conditions for existence of Hopf bifurcation are also investigated. Furthermore, an explicit algorithm determining the direction of Hopf bifurcations and stability of bifurcating periodic solutions occurring through Hopf bifurcations is also given.In the natural world, many species have a life history that takes them through two stages, juvenile stage and adult stage. Individuals in each stage are identical in biological characteristics, and some vital rates (rates of survival, development, and reproduction) of individuals in a population almost always depend on stage structure. Furthermore, many complex biological phenomena arising in prey predator ecosystem always depend on the past history of the system, and it has been recognized that time delay may have complicated impact on the dynamics of prey predator ecosystem [2]. In the past several decades, there has been an increasing interest in prey predator model with stage structure and time delay. In the model proposed by Aiello and Freedman [14], stage structure of single population growth with stage structure and time delay representing for maturation of population is considered. Their model predicts a positive steady state as the global attractor, thereby suggesting that stage structure does not generate sustained oscillations frequently observed in single population in the real world. Subsequent work made by other authors [15–20] suggest that time delay to adulthood should be state dependent. Al-Omari and Gourley [16] suggested that the time delay to adulthood should be state dependent and careful formulation of such state dependent time delays can lead to models that produce periodic solutions. Xu et al. [21] studied the persistence and stability of a delayed prey predator model with stage structure for predator. Gourley and Kuang [22] and Bandyopadhyay and Banerjee [23] formulated a class of general and robust prey predator models with stage structure and constant maturation time delay and performed a systematic mathematical and computational study. They have shown that there is a window in maturation time delay parameter that generates sustainable oscillatory dynamics. Ration dependent prey predator model is proposed in the work done in [12, 24–28] and permanence and stability analysis are also investigated.In model (4) with strong generic kernel proposed in [13], stage structure of prey population is not considered. It is well known that the immature prey population is usually considered as an easier target for predators compared with mature prey population [14]. Hence, the mature prey population predated by predator population can be ignored [1], and it is necessary to investigate dynamic effect of stage structure on prey predator model with discrete delay and distributed delay. By incorporating stage structure of prey population into model (4) with strong generic kernel, work done in [13] is extended in this paper. The organization of rest sections of this paper is as follows: a stage-structured prey predator model with discrete and distributed delay is proposed in the second section. It should be noted that the delay kernel function takes the strong generic kernel function and stage structure for prey population is considered in this paper. In the third section, qualitative analysis of the proposed model is investigated. Conditions for existence of two feasible boundary equilibria and a unique interior equilibrium are analytically discussed. By analyzing corresponding characteristic equation, local stability of two feasible boundary equilibria and an interior equilibrium are discussed, respectively. By taking discrete time delay as a parameter, Hopf bifurcation around the interior equilibrium is studied due to variation of discrete delay. Furthermore, directions of Hopf bifurcation and stability of the bifurcating periodic solutions are studied based on normal form theory and center manifold theorem. In the fourth section, numerical simulations are carried out to show consistency with theoretical results obtained in this paper. Finally, this paper ends with a conclusion.


## 2. Model Formulation

In this paper, a stage-structured prey predator model with discrete and distributed delay is investigated based on the following four hypotheses, which are given as follows.(H1)
Prey population is divided into two-stage groups, that is, immature prey populationx
1
(
t
) and mature prey population x
2
(
t
). The birth of immature prey population is proportional to the existing mature prey population with proportionality constant r
1
>
0 at any time t
>
0, and therefore the term r
1
x
2
(
t
) is in the first equation of model (5). a
11
>
0 and a
12
>
0 represent the intraspecific competition rate of prey population and predation coefficient of predator population, respectively.
(H2)
The rate of transformation of the mature prey population is proportional to the existing immature prey population with proportionality constantβ
>
0, which explains the term -
β
x
1
(
t
) in the first equation and β
x
1
(
t
) in the second equation. 0
<
r
3
<
r
1 denotes death rate of mature prey population.
(H3)
r
2
>
0 denotes death rate of predator population. It is assumed that the immature prey population is usually considered as an easier target for predators; the mature prey population predated by predator population can be ignored. The predator y
(
t
) is considered as a whole group, and only immature prey population is under predation from predator population. a
21
>
0 and a
22
>
0 denote biomass conversion rate of immature prey population captured by predator population and intraspecific competition rate of predator population, respectively.
(H4)
The delay kernel functionG
(
s
) takes strong generic kernel function G
(
s
)
=
α
2
s
e
-
α
s (α
>
0); the strong generic kernel implies that a particular time in the past is more important than any other [1, 2, 5].Based on hypotheses (H1)–(H4) and model (4), a stage-structured prey predator model with discrete hunting delay and distributed maturation delay is established as follows:
(5)
x
˙
1
(
t
)
=
r
1
x
2
(
t
)
-
a
11
x
1
2
(
t
)
-
a
12
x
1
(
t
)
y
(
t
-
τ
)
-
β
x
1
(
t
)
,
x
˙
2
(
t
)
=
β
x
1
(
t
)
-
r
3
x
2
(
t
)
,
y
˙
(
t
)
=
y
(
t
)
[
-
r
2
+
a
21
∫
-
∞
t
G
(
t
-
s
)
x
1
(
s
)
d
s
-
a
22
y
(
t
)
]
,
where x
1
(
t
), x
2
(
t
), and y
(
t
) denote the density of immature prey, mature prey, and predator population, respectively. Other parameters share the same interpretations introduced in hypotheses (H1)–(H4).


## 3. Qualitative Analysis of Model System

In this section, conditions for existence of two feasible boundary equilibria and a unique interior equilibrium are analytically investigated. By analyzing corresponding characteristic equation, local stability of model around boundary equilibrium and interior equilibrium is discussed, respectively. By taking discrete time delay as parameter, conditions for existence of Hopf bifurcation are discussed based on Hopf bifurcation theorem for functional differential equations. Furthermore, directions of Hopf bifurcation and stability of the bifurcating periodic solutions are studied by using normal form theory and center manifold theorem.


### 3.1. Stability Analysis of Equilibria

For model (5), it follows from normalized condition (3) that there are two feasible boundary equilibria and a unique interior equilibrium, which are given as follows:(i)
two feasible boundary equilibria:E
0
(
0,0
,
0
), E
1
(
β
(
r
1
-
r
3
)
/
a
11
r
3
,
β
2
(
r
1
-
r
3
)
/
a
11
r
3
2
,
0
);
(ii)
a unique interior equilibrium:E
*
(
x
1
*
,
x
2
*
,
y
*
) exists provided that
(6)
a
21
β
(
r
1
-
r
3
)
>
a
11
r
2
r
3
,
wherex
1
*= (
β
a
22
(
r
1
-
r
3
)
+
a
12
r
2
r
3
)
/
r
3
(
a
11
a
22
+
a
12
a
21
), x
2
* = β
[
β
a
22
(
r
1
-
r
3
)
+
a
12
r
2
r
3
]
/
r
3
2
(
a
11
a
22
+
a
12
a
21
), and y
* = (
β
a
21
(
r
1
-
r
3
)
-
a
11
r
2
r
3
)
/
r
3
(
a
11
a
22
+
a
12
a
21
).Theorem 1.
The local stability of model (5) around two-boundary equilibria is as follows:(i)
E
0 is a saddle node of model (5);
(ii)
E
1 is a saddle node of model (5) provided that E
* exists;
(iii)
E
1 is locally stable if 2
a
11
a
12
r
2
r
3
>
β
(
r
1
-
r
3
)
(
a
12
a
21
-
a
11
a
22
) and a
21
β
(
r
1
-
r
3
)
<
a
11
r
2
r
3.Proof.
By simple computation, the characteristic equation of model (5) around E
0 takes the following form:
(7)
(
λ
+
r
2
)
[
λ
2
+
(
β
+
r
3
)
λ
+
β
(
r
3
-
r
1
)
]
=
0
.
By solving (7), three eigenvalues of model (5) around E
0 can be obtained as follows: λ
1
=
-
r
2, λ
2
=
(
-
(
β
+
r
3
)
+
(
β
+
r
3
)
2
+
4
β
(
r
1
-
r
3
)
)
/
2, and λ
3
=
(
-
(
β
+
r
3
)
-
(
β
+
r
3
)
2
+
4
β
(
r
1
-
r
3
)
)
/
2. According to r
1
>
r
3 introduced in (H2), it derives that Re
λ
1
<
0, Re
λ
2
>
0, and Re
λ
3
<
0. Hence, E
0 is a saddle node of model (5).
By simple computation, the characteristic equation of model (5) around E
1 takes the following form:
(8)
(
λ
+
r
2
-
a
21
x
1
*
)
[
λ
2
+
(
r
3
+
β
+
2
a
11
x
1
*
)
λ
k
k
k
k
k
k
k
k
k
k
k
i
k
k
+
r
3
(
β
+
2
a
11
x
1
*
)
-
β
r
1
]
=
0
.
By solving (8), three eigenvalues of model (5) around E
1 are as follows: λ
1
=
a
22
[
a
21
β
(
r
1
-
r
3
)
-
a
11
r
2
r
3
]
/
r
3
(
a
11
a
22
+
a
12
a
21
), which derives that Re
λ
1
>
0 if E
* exists; that is, a
21
β
(
r
1
-
r
3
)
>
a
11
r
2
r
3 holds; Re
λ
1
<
0 if a
21
β
(
r
1
-
r
3
)
<
a
11
r
2
r
3.
Furthermore, the other two eigenvalues of (8) can be determined by the following equation:
(9)
λ
2
+
(
r
3
+
β
+
2
a
11
x
1
*
)
λ
+
r
3
(
β
+
2
a
11
x
1
*
)
-
β
r
1
=
0
.
It follows from Routh-Hurwitz criterion [29] that Re
λ
2
<
0, Re
λ
3
<
0 provided that 2
a
11
a
12
r
2
r
3
>
β
(
r
1
-
r
3
)
(
a
12
a
21
-
a
11
a
22
).
Consequently,E
1 is a saddle node of model (5) provided that E
* exists; E
1 is locally stable if 2
a
11
a
12
r
2
r
3
>
β
(
r
1
-
r
3
)
(
a
12
a
21
-
a
11
a
22
) and a
21
β
(
r
1
-
r
3
)
<
a
11
r
2
r
3.Define new variablesu
(
t
) and v
(
t
) as follows:
(10)
u
(
t
)
=
∫
-
∞
t
α
2
(
t
-
s
)
e
-
α
(
t
-
s
)
x
1
(
s
)
d
s
,
v
(
t
)
=
∫
-
∞
t
α
e
-
α
(
t
-
s
)
x
1
(
s
)
d
s
.According to the law of solving the derivative for an integral with parameterized variables, it can be obtained that(11)
u
˙
(
t
)
=
α
v
(
t
)
-
α
u
(
t
)
,
v
˙
(
t
)
=
α
2
x
1
(
t
)
-
α
2
v
(
t
)
.By virtue of (11), model (5) can be transformed into the following five-dimensional system of FDEs with a discrete delay:
(12)
x
˙
1
(
t
)
=
r
1
x
2
(
t
)
-
a
11
x
1
2
(
t
)
-
a
12
x
1
(
t
)
y
(
t
-
τ
)
-
β
x
1
(
t
)
,
x
˙
2
(
t
)
=
β
x
1
(
t
)
-
r
3
x
2
(
t
)
,
y
˙
(
t
)
=
y
(
t
)
[
-
r
2
+
a
21
u
(
t
)
-
a
22
y
(
t
)
]
,
u
˙
(
t
)
=
α
v
(
t
)
-
α
u
(
t
)
,
v
˙
(
t
)
=
α
2
x
1
(
t
)
-
α
2
v
(
t
)
.It follows from model (12) that u
*
=
v
*
=
x
1
* for the interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
,
u
*
,
v
*
) of model (12), which will be denoted as E
*
(
x
1
*
,
x
2
*
,
y
*
,
x
1
*
,
x
1
*
) in the following part of this paper. Furthermore, the interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
) of model (5) is transformed into the interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
,
x
1
*
,
x
1
*
) of model (12). Hence, the qualitative local stability analysis around E
*
(
x
1
*
,
x
2
*
,
y
*
,
x
1
*
,
x
1
*
) of model (12) is equivalent to the qualitative local stability analysis around interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
) of model (5).According to the leading matrix of model (12) around the interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
,
x
1
*
,
x
1
*
), the characteristic equation takes the following form:
(13)
|
λ
+
2
a
11
x
1
*
+
a
12
y
*
+
β
-
r
1
a
12
x
1
*
e
-
λ
τ
0
0
-
β
λ
+
r
3
0
0
0
0
0
λ
+
a
22
y
*
-
a
21
y
*
0
0
0
0
λ
+
α
-
α
-
α
2
0
0
0
λ
+
α
2
|
=
0
,
which derives that
(14)
M
(
λ
)
+
N
(
λ
)
e
-
λ
τ
=
0
,
where
(15)
M
(
λ
)
=
λ
5
+
m
1
λ
4
+
m
2
λ
3
+
m
3
λ
2
+
m
4
λ
+
m
5
,
N
(
λ
)
=
n
4
λ
+
n
5
,
m
1
=
2
a
11
x
1
*
+
(
a
12
+
a
22
)
y
*
+
β
+
α
+
α
2
+
r
3
,
m
2
=
r
3
(
2
a
11
x
1
*
+
a
12
y
*
)
+
β
(
r
3
-
r
1
)
+
α
3
+
a
22
y
*
(
α
+
α
2
)
+
(
a
22
y
*
+
α
+
α
2
)
×
(
2
a
11
x
1
*
+
a
12
y
*
+
β
+
r
3
)
,
m
3
=
(
a
22
y
*
+
α
+
α
2
)
×
[
r
3
(
2
a
11
x
1
*
+
a
12
y
*
)
+
β
(
r
3
-
r
1
)
]
+
(
2
a
11
x
1
*
+
a
12
y
*
+
β
+
r
3
)
×
[
α
3
+
a
22
y
*
(
α
+
α
2
)
]
+
a
22
α
3
y
*
,
m
4
=
[
r
3
(
2
a
11
x
1
*
+
a
12
y
*
)
+
β
(
r
3
-
r
1
)
]
[
α
3
+
a
22
y
*
(
α
+
α
2
)
]
+
a
22
α
3
y
*
(
2
a
11
x
1
*
+
a
12
y
*
+
β
+
r
3
)
,
m
5
=
a
22
α
3
y
*
[
r
3
(
2
a
11
x
1
*
+
a
12
y
*
)
+
β
(
r
3
-
r
1
)
]
,
n
4
=
a
12
a
21
α
3
x
1
*
y
*
,
n
5
=
a
12
a
21
r
3
α
3
x
1
*
y
*
.Whenτ
=
0, (14) takes the following form:
(16)
λ
5
+
m
1
λ
4
+
m
2
λ
3
+
m
3
λ
2
+
(
m
4
+
n
4
)
λ
+
m
5
+
n
5
=
0
.Based on the values ofm
k, n
j, k
=
1, 2, 3, 4, 5, j
=
4, 5 defined in (14), it is easy to show that m
1
>
0 and m
5
+
n
5
>
0.Define quantitiesδ
k, k
=
1,2
,
3,4
,
5 as follows:
(17)
δ
k
=
|
m
1
1
0
0
0
m
3
m
2
m
1
1
0
m
5
+
n
5
m
4
+
n
4
m
3
m
2
m
1
0
0
m
5
+
n
5
m
4
+
n
4
m
3
0
0
0
0
m
5
+
n
5
|
=
a
k
δ
k
-
1
.Further computations show thatδ
k
>
0 for k
=
2,3
,
4. According to the above analysis and Routh-Hurwitz criteria [2], it follows from m
1
>
0, m
5
+
n
5
>
0, and δ
k
>
0 for k
=
2,3
,
4 that model (12) is locally stable around E
*
(
x
1
*
,
x
2
*
,
y
*
,
x
1
*
,
x
1
*
) in the case of τ
=
0, which implies that model (5) is locally stable around E
*
(
x
1
*
,
x
2
*
,
y
*
).Subsequently, the dynamic effect of discrete delay on model dynamics is investigated. In the case of positive discrete delayτ
>
0, let i
ω (ω
>
0) be a root of (14). Separating the real and imaginary parts of (14) with λ
=
i
ω yields the following equations:
(18)
ω
5
+
m
2
ω
3
+
m
4
ω
=
n
5
sin
⁡
ω
τ
-
n
4
ω
cos
⁡
⁡
ω
τ
,
m
3
ω
2
-
m
1
ω
4
-
m
5
=
n
5
cos
⁡
⁡
ω
τ
+
n
4
ω
sin
⁡
ω
τ
.By adding up the squares of the corresponding sides of (18), it yields the following algebraic equation with respect to ω:
(19)
ω
10
+
A
1
ω
8
+
A
2
ω
6
+
A
3
ω
4
+
A
4
ω
2
+
A
5
=
0
,
where A
1
=
2
m
2
+
m
1
2, A
2
=
m
2
2
+
2
m
4
-
2
m
1
m
3, A
3
=
m
3
2
+
2
(
m
2
m
4
+
m
1
m
5
), A
4
=
m
4
2
-
2
m
3
m
5
-
n
4
2, and A
5
=
m
5
2
-
n
5
2.Letξ
=
ω
2; (19) can be rewritten as follows:
(20)
g
(
ξ
)
∶
=
ξ
5
+
A
1
ξ
4
+
A
2
ξ
3
+
A
3
ξ
2
+
A
4
ξ
+
A
5
=
0
.By virtue of derivative ofg
(
ξ
) with respect to ξ, it gives that
(21)
g
′
(
ξ
)
=
5
ξ
4
+
4
A
1
ξ
3
+
3
A
2
ξ
2
+
2
A
3
ξ
+
A
4
.By simple computation, it is easy to show thatA
k
>
0, k
=
1, 2, 3, 4, which derives that g
′
(
ξ
)
>
0 on (
0
,
+
∞
) and g
(
ξ
) is strictly monotonically increasing on (
0
,
+
∞
). Consequently, g
(
ξ
) has a unique positive root if A
5
<
0 and g
(
ξ
) does not have a positive root if A
5
≥
0.SupposeA
5
<
0, which derives that
(22)
(
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
>
0
,
and then there exists a unique root of g
(
ξ
)
=
0, which is denoted by ξ
0, and the unique positive root of (14) is ω
0
=
ξ
0.SupposeA
5
≥
0, which derives that
(23)
(
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
≤
0
,
and then (14) does not have a positive root.By eliminatingsin
⁡
ω
τ, it follows from (18) that
(24)
τ
k
=
1
ω
0
[
+
2
k
π
arccos
(
(
n
4
2
ω
0
2
+
n
5
2
)
-
1
(
-
n
4
ω
0
6
-
(
m
1
n
5
+
m
2
n
4
)
ω
0
4
kkkkkkkkkkkkk
+
(
m
3
n
5
-
m
4
n
4
)
ω
0
2
-
m
5
n
5
)
h
h
h
h
h
h
h
h
k
i
h
h
×
(
n
4
2
ω
0
2
+
n
5
2
)
-
1
)
kkkkk
+
2
k
π
(
(
n
4
2
ω
0
2
+
n
5
2
)
-
1
(
-
n
4
ω
0
6
-
(
m
1
n
5
+
m
2
n
4
)
ω
0
4
]
,
and then (14) has a pair of purely imaginary roots ±
i
ω
0 when τ
=
τ
k for k
=
0,1
,
….By using Butler’s lemma [29], the interior equilibrium E
* remains stable for 0
≤
τ
<
τ
0. Based on the above analysis, local stability analysis of model (5) around the interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
) can be concluded in the following Theorem 2.Theorem 2.
If condition (6) holds, then the following results hold.(i)
If(
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
≤
0, then model (5) is locally stable around E
*
(
x
1
*
,
x
2
*
,
y
*
) for any values τ
≥
0.
(ii)
If(
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
>
0, then model (5) is locally stable around E
*
(
x
1
*
,
x
2
*
,
y
*
) for τ
∈
[
0
,
τ
0
).
Subsequently, the conditions for Hopf bifurcation in [29] are utilized to investigate whether there is a phenomenon of Hopf bifurcation as τ increases through τ
k.Theorem 3.
If(
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
>
0, then model (5) undergoes Hopf bifurcation at the interior equilibrium E
* when τ
=
τ
k, k
=
0,1
,
2
,
…. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium E
* when τ
>
τ
0 and ∥
τ
-
τ
0
∥
≪
1.Proof.
As mentioned above, letλ
=
i
ω
0 represent the purely imaginary root of (14). It follows from (14) that |
M
(
i
ω
0
)
|
=
|
N
(
i
ω
0
)
|, which determines a set of possible values of ω
0.
In the following part, we determine the direction of motion ofλ
=
i
ω
0 as τ is varied; namely, we determine
(25)
Θ
=
sign
⁡
[
d
(
Re
λ
)
d
τ
]
λ
=
i
ω
0
=
sign
⁡
[
Re
(
d
λ
d
τ
)
-
1
]
λ
=
i
ω
0
.
By differentiating (14) with respect to τ, it can be obtained that
(26)
[
5
λ
4
+
4
m
1
λ
3
+
3
m
2
λ
2
+
2
m
3
λ
+
m
4
+
n
4
e
-
λ
τ
k
k
-
τ
(
n
4
λ
+
n
5
)
e
-
λ
τ
]
d
λ
d
τ
=
λ
(
n
4
λ
+
n
5
)
e
-
λ
τ
,
⟹
(
d
λ
d
τ
)
-
1
=
5
λ
4
+
4
m
1
λ
3
+
3
m
2
λ
2
+
2
m
3
λ
+
m
4
λ
(
n
4
λ
+
n
5
)
e
-
λ
τ
+
n
4
λ
(
n
4
λ
+
n
5
)
-
τ
λ
=
5
λ
4
+
4
m
1
λ
3
+
3
m
2
λ
2
+
2
m
3
λ
+
m
4
-
λ
(
λ
5
+
m
1
λ
4
+
m
2
λ
3
+
m
3
λ
2
+
m
4
λ
+
m
5
)
+
n
4
λ
(
n
4
λ
+
n
5
)
-
τ
λ
=
4
λ
5
+
3
m
1
λ
4
+
2
m
2
λ
3
+
m
3
λ
2
-
m
5
-
λ
2

(
λ
5
+
m
1
λ
4
+
m
2
λ
3
+
m
3
λ
2
+
m
4
λ
+
m
5
)
-
n
5
λ
2
(
n
4
λ
+
n
5
)
-
τ
λ
,
Θ
=
sign
⁡
[
Re
(
d
λ
d
τ
)
-
1
]
λ
=
i
ω
0
=
1
ω
0
2
sign
⁡
[
Re
(
n
5
n
4
ω
0
i
+
n
5
(
3
m
1
ω
0
4
-
m
3
ω
0
2
-
m
5
k
k
k
k
k
k
k
k
k
k
k
k
h
+
2
ω
0
3
(
2
ω
0
2
-
m
2
)
i
)
k
k
k
k
k
k
k
k
k
k
k
k
k
×
(
m
1
ω
0
4
-
m
3
ω
0
2
+
m
5
k
k
k
k
k
h
k
k
k
k
k
k
k
+
(
ω
0
5
-
m
2
ω
0
3
+
m
4
ω
0
)
i
)
-
1
k
k
k
k
k
k
k
k
k
k
k
k
+
n
5
n
4
ω
0
i
+
n
5
+
ω
0
τ
i
)
]
=
1
ω
0
2
sign
⁡
[
+
(
ω
0
5
-
m
2
ω
0
3
+
m
4
ω
0
)
2
)
-
1
(
3
m
1
2
ω
0
8
+
4
m
1
m
3
ω
0
6
k
k
k
k
k
k
k
k
k
k
+
(
2
m
1
m
5
+
m
3
2
)
ω
0
4
+
n
5
2
-
m
5
2
)
k
k
k
k
k
k
k
k
k
×
(
+
(
ω
0
5
-
m
2
ω
0
3
+
m
4
ω
0
)
2
(
m
1
ω
0
4
-
m
3
ω
0
2
+
m
5
)
2
k
k
k
k
k
k
k
k
k
k
k
+
(
ω
0
5
-
m
2
ω
0
3
+
m
4
ω
0
)
2
)
-
1
]
.
If(
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
>
0, then n
5
-
m
5
>
0 which signifies that
(27)
Θ
=
sign
⁡
[
Re
(
d
λ
d
τ
)
-
1
]
λ
=
i
ω
0
>
0
.
Hence, the transversality condition [29] holds and model (5) undergoes Hopf bifurcation at the interior equilibrium E
* when τ
=
τ
k, k
=
0,1
,
2
,
…. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium E
* when τ
>
τ
0 and ∥
τ
-
τ
0
∥
≪
1.


### 3.2. Direction and Stability of Hopf Bifurcation

By using normal theory and center manifold theorem [29, 30], directions of Hopf bifurcation and stability of the bifurcating periodic solutions are discussed in this section. As analyzed in Section 3.1, some transformations associated with the interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
,
x
1
*
,
x
1
*
) of model (12) are given as follows:
(28)
z
1
=
x
1
-
x
1
*
,
z
2
=
x
2
-
x
2
*
,
z
3
=
y
-
y
*
,
z
4
=
u
-
x
1
*
,
z
5
=
v
-
x
1
*
,
z
¯
i
(
t
)
=
z
i
(
τ
t
)
,
τ
=
μ
+
τ
k
,
μ
∈
ℝ
=
(
-
∞
,
+
∞
)
,
and then μ
=
0 is the Hopf bifurcation value of model (12). Bars of variables are dropped for simplicity of notations; model (12) is transformed to a functional differential equation in C
=
C
(
[
-
1,0
]
,
ℝ
3
),
(29)
z
˙
(
t
)
=
L
μ
(
z
t
)
+
f
(
μ
,
z
t
)
,
where C
=
C
(
[
-
1,0
]
,
ℝ
5
) is the Banach space of continuous functions mapping [
-
τ
,
0
] into ℝ
5, z
(
t
)
=
(
z
1
(
t
)
,
z
2
(
t
)
,
z
3
(
t
)
,
z
4
(
t
)
,
z
5
(
t
)
)
T
∈
ℝ
5, z
t
(
θ
)
=
z
(
t
+
θ
) for θ
∈
[
-
τ
,
0
], and L
μ
:
C
→
ℝ
5, f
:
ℝ
×
C
→
ℝ
5 are defined as follows, respectively:
(30)
L
μ
(
ϕ
)
=
(
τ
k
+
μ
)
×
(
-
a
11
x
1
*
-
r
1
x
2
*
x
1
*
r
1
0
0
0
β
-
r
3
0
0
0
0
0
-
a
22
y
*
a
21
y
*
0
0
0
0
-
α
α
α
2
0
0
0
-
α
2
)
×
(
ϕ
1
(
0
)
ϕ
2
(
0
)
ϕ
3
(
0
)
ϕ
4
(
0
)
ϕ
5
(
0
)
)
+
(
τ
k
+
μ
)
(
0
0
-
a
12
x
1
*
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
)
(
ϕ
1
(
-
1
)
ϕ
2
(
-
1
)
ϕ
3
(
-
1
)
ϕ
4
(
-
1
)
ϕ
5
(
-
1
)
)
,
(31)
f
(
μ
,
ϕ
)
=
(
τ
k
+
μ
)
(
-
a
12
ϕ
1
(
0
)
ϕ
3
(
-
1
)
0
0
0
0
)
.It is easy to show thatL
μ is a continuous linear function mapping C into ℝ
5. According to Riesz representation theorem [4], there exists a 5
×
5 matrix function η
(
θ
,
μ
) of bounded variation for θ
∈
[
-
1,0
] such that
(32)
L
μ
(
ϕ
)
=
∫
-
1
0
d
η
(
θ
,
μ
)
ϕ
(
θ
)
,
where ϕ
∈
C
(
[
-
1,0
]
,
ℝ
5
).In fact, we can choose(33)
η
(
θ
,
μ
)
=
(
τ
k
+
μ
)
×
(
-
a
11
x
1
*
-
r
1
x
2
*
x
1
*
r
1
0
0
0
β
-
r
3
0
0
0
0
0
-
a
22
y
*
a
21
y
*
0
0
0
0
-
α
α
α
2
0
0
0
-
α
2
)
×
δ
(
θ
)
-
(
τ
k
+
μ
)
(
0
0
-
a
12
x
1
*
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
)
×
δ
(
θ
+
1
)
,
where δ denotes the Dirac delta function.Ifϕ is any given function in C
(
[
-
1,0
]
,
ℝ
5
) and z
(
ϕ
) is the unique solution of the linearized equation z
˙
(
t
)
=
L
μ
(
z
t
) of (29) with initial function ϕ at zero, then the solution operator T
~
(
t
)
:
C
→
C is defined by
(34)
T
~
(
t
)
ϕ
=
z
t
(
ϕ
)
,
t
≥
0
.It follows from Lemma  7.1.1 in [29] that T
~
(
t
), t
≥
0 is a strongly continuous semigroup of linear transformation on [
0
,
+
∞
) and the infinitesimal generator A
μ of T
~
(
t
), t
≥
0 is as follows:
(35)
A
μ
(
ϕ
)
=
{
d
ϕ
(
θ
)
d
θ
,
θ
∈
[
-
1,0
)
∫
-
1
0
d
η
(
μ
,
s
)
ϕ
(
s
)
,
θ
=
0
,
for ϕ
∈
C
1
(
[
-
1,0
]
,
ℝ
5
), the space of functions mapping the interval [
-
1,0
] into ℝ
5 which have continuous first derivative, and also define
(36)
R
(
μ
)
(
ϕ
)
=
{
0
,
θ
∈
[
-
1,0
)
f
(
μ
,
ϕ
)
,
θ
=
0
,
and then model (29) is equivalent to
(37)
z
˙
t
=
A
(
μ
)
z
t
+
R
(
μ
)
z
t
.Forψ
∈
C
1
(
[
0,1
]
,
(
ℝ
5
)
*
), the space of functions mapping interval [
0,1
] into the five-dimensional row vectors which have continuous first derivative, define
(38)
A
*
ϕ
(
s
)
=
{
-
d
ψ
(
s
)
d
s
,
s
∈
(
0,1
]
∫
-
1
0
d
η
T
(
t
,
0
)
ψ
(
-
t
)
,
s
=
0
,
and a bilinear inner product
(39)
〈
ψ
(
s
)
,
ϕ
(
θ
)
〉
=
ψ
¯
(
0
)
ϕ
(
0
)
-
∫
-
1
0
∫
ρ
=
0
θ
ψ
¯
(
ρ
-
θ
)
d
η
(
θ
)
ϕ
(
ρ
)
d
ρ
,
where η
(
θ
)
=
η
(
θ
,
0
). It follows from the above analysis that A
(
0
) and A
* are adjoint operators.By virtue of discussion in Section3.1, ±
i
ω
0
τ
k are eigenvalues of A
(
0
). Hence, they are also eigenvalues of A
*. In the following, eigenvectors of A
(
0
) and A
* correspond to i
ω
0
τ
k and -
i
ω
0
τ
k, respectively.Supposeq
(
θ
)
=
(
1
,
a
,
b
,
c
,
d
)
T
e
i
ω
0
τ
k
θ is the eigenvectors of A
(
0
) corresponding to i
ω
0
τ
k, which derives that A
(
0
)
q
(
θ
)
=
i
ω
0
τ
k
q
(
θ
). By using the definition of A
(
0
), (30), and (31), it gives that
(40)
τ
k
(
-
a
11
x
1
*
-
r
1
x
2
*
x
1
*
r
1
0
0
0
β
-
r
3
0
0
0
0
0
-
a
22
y
*
a
21
y
*
0
0
0
0
-
α
α
α
2
0
0
0
-
α
2
)
q
(
0
)
+
τ
k
(
0
0
-
a
12
x
1
*
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
)
q
(
-
1
)
=
i
ω
0
τ
k
q
(
0
)
.Forq
(
-
1
)
=
q
(
0
)
e
-
i
ω
0
τ
k, then it can be obtained that
(41)
a
=
β
r
3
+
i
ω
0
,
b
=
x
1
*
[
β
r
1
-
(
a
11
x
1
*
+
i
ω
0
)
(
r
3
+
i
ω
0
)
]
-
r
1
x
2
*
(
r
3
+
i
ω
0
)
a
12
x
1
*
2
e
-
i
ω
0
τ
k
(
r
3
+
i
ω
0
)
,
c
=
α
3
(
α
+
i
ω
0
)
(
α
2
+
i
ω
0
)
,
d
=
α
2
α
2
+
i
ω
0
.Similarly, it follows from simple computation that eigenvector ofA
* corresponds to -
i
ω
0
τ
k; that is, q
*
(
s
)
=
J
(
1
,
a
*
,
b
*
,
c
*
,
d
*
)
T
e
i
ω
0
τ
k
s
(42)
a
*
=
r
1
r
3
+
i
ω
0
,
b
*
=
-
a
12
x
1
*
e
-
i
ω
0
τ
k
a
22
y
*
+
i
ω
0
,
c
*
=
-
a
12
a
21
x
1
*
e
-
i
ω
0
τ
k
(
α
+
i
ω
0
)
(
a
22
y
*
+
i
ω
0
)
,
d
*
=
-
α
a
12
a
21
x
1
*
e
-
i
ω
0
τ
k
(
α
+
i
ω
0
)
(
α
2
+
i
ω
0
)
(
a
22
y
*
+
i
ω
0
)
.In order to assume〈
q
*
(
s
)
,
q
(
θ
)
〉
=
1, we need to determine the value of J in the following part.By virtue of(A.9), it derives that
(43)
〈
q
*
(
s
)
,
q
(
θ
)
〉
=
J
¯
(
1
,
a
¯
*
,
b
¯
*
,
c
¯
*
,
d
¯
*
)
(
1
,
a
,
b
,
c
,
d
)
T
-
∫
-
1
0
∫
ρ
=
0
θ
J
¯
(
1
,
a
*
,
b
*
,
c
*
,
d
*
)
e
-
i
ω
0
τ
k
(
ρ
-
θ
)
d
η
(
θ
)
k
k
k
k
k
k
k
k
h
k
k
×
(
1
,
a
,
b
,
c
,
d
)
T
e
i
ω
0
τ
k
ρ
d
ρ
=
J
¯
(
1
+
a
a
¯
*
+
b
b
¯
*
+
c
c
¯
*
+
d
d
¯
*
)
-
J
¯
∫
-
1
0
(
1
,
a
¯
*
,
b
¯
*
,
c
¯
*
,
d
¯
*
)
θ
e
i
ω
0
τ
k
θ
d
η
(
θ
)
(
1
,
a
,
b
,
c
,
d
)
T
=
J
¯
(
1
+
a
a
¯
*
+
b
b
¯
*
+
c
c
¯
*
+
d
d
¯
*
+
B
)
,
where
(44)
B
=
[
r
1
x
2
*
x
1
*
α
2
d
¯
*
+
β
a
¯
*
+
a
(
r
1
-
a
¯
*
(
r
3
+
i
ω
0
)
)
k
-
a
11
x
1
*
-
i
ω
0
-
r
1
x
2
*
x
1
*
]
-
b
[
a
12
x
1
*
e
-
i
ω
0
τ
k
+
b
¯
*
(
a
22
y
*
+
i
ω
0
)
]
+
c
[
a
21
b
¯
*

y
*
-
c
¯
*
(
α
+
i
ω
0
)
]
+
d
[
α
c
¯
*
-
(
α
2
+
i
ω
0
)
d
¯
*
]
.Hence, we can chooseJ as follows:
(45)
J
=
1
1
+
a
a
¯
*
+
b
b
¯
*
+
c
c
¯
*
+
d
d
¯
*
+
B
.Nextly, we will compute the coordinate to describe the centre manifoldC
0 at μ
=
0. Let z
t be the solution of (37) when μ
=
0.Define(46)
χ
(
t
)
=
〈
q
*
,
z
t
〉
,
W
(
t
,
θ
)
=
z
t
(
θ
)
-
2
Re
{
χ
(
t
)
q
(
θ
)
}
.On the center manifoldC
0, it derives that
(47)
W
(
t
,
θ
)
=
W
(
χ
(
t
)
,
χ
¯
(
t
)
,
θ
)
,
where
(48)
W
(
χ
(
t
)
,
χ
¯
(
t
)
,
θ
)
=
W
20
(
θ
)
χ
2
2
+
W
11
(
θ
)
χ
χ
¯
+
W
02
χ
¯
2
2
+
⋯
,
χ and χ
¯ are local coordinates for center manifold C
0 in the direction of q
* and q
¯
*.It is noted thatW is real if z
t is real, and we only consider real solutions. For solution z
t
∈
C
0 of (37), since μ
=
0, it derives that
(49)
χ
˙
(
t
)
=
i
ω
0
τ
k
χ
+
q
¯
*
(
0
)
f
(
0
,
W
(
χ
,
χ
¯
,
0
)
+
2
Re
{
χ
q
(
θ
)
}
)
≜
i
ω
0
τ
k
χ
+
q
¯
*
(
0
)
f
0
(
χ
,
χ
¯
)
.The above equation can be rewritten as follows:(50)
χ
˙
(
t
)
=
i
ω
0
τ
k
χ
(
t
)
+
g
(
χ
,
χ
¯
)
,
where
(51)
g
(
χ
,
χ
¯
)
=
q
¯
*
(
0
)
f
0
(
χ
,
χ
¯
)
=
g
20
χ
2
2
+
g
11
χ
χ
¯
+
g
02
χ
¯
2
2
+
g
21
χ
2
χ
¯
2
+
⋯
.It follows from (46) and (48) that
(52)
z
t
(
θ
)
=
W
(
t
,
θ
)
+
2
Re
{
χ
(
t
)
q
(
θ
)
}
=
W
20
(
θ
)
χ
2
2
+
W
11
(
θ
)
χ
χ
¯
+
W
02
(
θ
)
χ
¯
2
2
+
(
1
,
a
,
b
,
c
,
d
)
T
e
i
ω
0
τ
k
θ
χ
+
(
1
,
a
¯
,
b
¯
,
c
¯
,
d
¯
)
T
e
-
i
ω
0
τ
k
θ
χ
¯
+
⋯
.By virtue of (31), (32), and (33), it derives that
(53)
g
(
χ
,
χ
¯
)
=
q
¯
*
(
0
)
f
0
(
χ
,
χ
¯
)
=
q
¯
*
(
0
)
f
(
0
,
z
t
)
=
τ
k
J
¯
(
1
,
a
¯
*
,
b
¯
*
,
c
¯
*
,
d
¯
*
)
(
-
a
12
z
1
t
(
0
)
z
3
t
(
-
1
)
0
0
0
0
)
=
-
τ
k
J
¯
a
12
z
1
t
(
0
)
z
3
t
(
-
1
)
=
-
τ
k
J
¯
a
12
[
χ
+
χ
¯
+
W
20
(
1
)
(
0
)
χ
2
2
+
W
11
(
1
)
(
0
)
χ
χ
¯
kkkkkkkkk
+
W
02
(
1
)
(
0
)
χ
¯
2
2
+
o
(
|
(
χ
,
χ
¯
)
|
3
)
]
×
[
χ
+
χ
¯
+
W
20
(
3
)
(
-
1
)
χ
2
2
+
W
11
(
3
)
(
-
1
)
χ
χ
¯
k
k
k
k
k
+
W
02
(
3
)
(
-
1
)
χ
¯
2
2
+
o
(
|
(
χ
,
χ
¯
)
|
3
)
]
=
-
τ
k
J
¯
a
12
[
+
W
20
(
1
)
(
0
)
+
W
20
(
3
)
(
-
1
)
2
)
χ
2
+
2
χ
χ
¯
+
χ
¯
2
k
k
k
k
h
k
k
k
k
+
(
W
20
(
1
)
(
0
)
+
W
20
(
3
)
(
-
1
)
2
W
11
(
3
)
(
-
1
)
+
W
11
(
1
)
(
0
)
k
k
k
k
k
k
k
k
k
k
k
+
W
20
(
1
)
(
0
)
+
W
20
(
3
)
(
-
1
)
2
)
χ
2
χ
¯
+
⋯
]
.By comparing the coefficients with (51), it gives that
(54)
g
20
=
-
2
a
12
τ
k
J
¯
,
g
11
=
-
2
a
12
τ
k
J
¯
,
g
02
=
-
2
a
12
τ
k
J
¯
,
g
21
=
-
2
a
12
τ
k
J
¯
×
(
W
11
(
3
)
(
-
1
)
+
W
11
(
1
)
(
0
)
+
W
20
(
1
)
(
0
)
+
W
20
(
3
)
(
-
1
)
2
)
.Sinceg
21 is associated with W
20
(
θ
) and W
11
(
θ
), further attempts should be carried out to compute W
20
(
θ
) and W
11
(
θ
), which can be found in the appendix.Furthermore, we can computeg
21 based on (54). Hence, the following values can be computed as follows:
(55)
d
1
(
0
)
=
i
2
ω
0
τ
k
(
g
20
g
11
-
2
|
g
11
|
2
-
|
g
02
|
2
3
)
+
g
21
2
,
δ
2
=
-
Re
{
d
1
(
0
)
}
Re
{
λ
′
(
τ
k
)
}
,
γ
2
=
2
Re
(
d
1
(
0
)
)
,
T
2
=
Im
⁡
{
d
1
(
0
)
}
+
δ
2
Im
⁡
{
λ
′
(
τ
k
)
}
ω
0
τ
k
.Some properties of bifurcating periodic solution in the center manifold at the critical valuesτ
k follow from [29]. Based on the analysis in Section 3.1 of this paper, the following theorem can be concluded.Theorem 4.
The properties of Hopf bifurcation are determined by values computed in (55).(i)
μ
2 determines directions of Hopf bifurcation: if μ
2
>
0  (
μ
2
<
0
), then Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for τ
>
τ
k
(
τ
<
τ
k
).
(ii)
γ
2 determines the stability of bifurcating periodic solutions: bifurcating periodic solutions are stable (unstable) if γ
2
<
0
(
γ
2
>
0
).
(iii)
T
2 determines the period of bifurcating periodic solutions: period increases (decreases) if T
2
>
0
(
T
2
<
0
).


## 3.1. Stability Analysis of Equilibria

For model (5), it follows from normalized condition (3) that there are two feasible boundary equilibria and a unique interior equilibrium, which are given as follows:(i)
two feasible boundary equilibria:E
0
(
0,0
,
0
), E
1
(
β
(
r
1
-
r
3
)
/
a
11
r
3
,
β
2
(
r
1
-
r
3
)
/
a
11
r
3
2
,
0
);
(ii)
a unique interior equilibrium:E
*
(
x
1
*
,
x
2
*
,
y
*
) exists provided that
(6)
a
21
β
(
r
1
-
r
3
)
>
a
11
r
2
r
3
,
wherex
1
*= (
β
a
22
(
r
1
-
r
3
)
+
a
12
r
2
r
3
)
/
r
3
(
a
11
a
22
+
a
12
a
21
), x
2
* = β
[
β
a
22
(
r
1
-
r
3
)
+
a
12
r
2
r
3
]
/
r
3
2
(
a
11
a
22
+
a
12
a
21
), and y
* = (
β
a
21
(
r
1
-
r
3
)
-
a
11
r
2
r
3
)
/
r
3
(
a
11
a
22
+
a
12
a
21
).Theorem 1.
The local stability of model (5) around two-boundary equilibria is as follows:(i)
E
0 is a saddle node of model (5);
(ii)
E
1 is a saddle node of model (5) provided that E
* exists;
(iii)
E
1 is locally stable if 2
a
11
a
12
r
2
r
3
>
β
(
r
1
-
r
3
)
(
a
12
a
21
-
a
11
a
22
) and a
21
β
(
r
1
-
r
3
)
<
a
11
r
2
r
3.Proof.
By simple computation, the characteristic equation of model (5) around E
0 takes the following form:
(7)
(
λ
+
r
2
)
[
λ
2
+
(
β
+
r
3
)
λ
+
β
(
r
3
-
r
1
)
]
=
0
.
By solving (7), three eigenvalues of model (5) around E
0 can be obtained as follows: λ
1
=
-
r
2, λ
2
=
(
-
(
β
+
r
3
)
+
(
β
+
r
3
)
2
+
4
β
(
r
1
-
r
3
)
)
/
2, and λ
3
=
(
-
(
β
+
r
3
)
-
(
β
+
r
3
)
2
+
4
β
(
r
1
-
r
3
)
)
/
2. According to r
1
>
r
3 introduced in (H2), it derives that Re
λ
1
<
0, Re
λ
2
>
0, and Re
λ
3
<
0. Hence, E
0 is a saddle node of model (5).
By simple computation, the characteristic equation of model (5) around E
1 takes the following form:
(8)
(
λ
+
r
2
-
a
21
x
1
*
)
[
λ
2
+
(
r
3
+
β
+
2
a
11
x
1
*
)
λ
k
k
k
k
k
k
k
k
k
k
k
i
k
k
+
r
3
(
β
+
2
a
11
x
1
*
)
-
β
r
1
]
=
0
.
By solving (8), three eigenvalues of model (5) around E
1 are as follows: λ
1
=
a
22
[
a
21
β
(
r
1
-
r
3
)
-
a
11
r
2
r
3
]
/
r
3
(
a
11
a
22
+
a
12
a
21
), which derives that Re
λ
1
>
0 if E
* exists; that is, a
21
β
(
r
1
-
r
3
)
>
a
11
r
2
r
3 holds; Re
λ
1
<
0 if a
21
β
(
r
1
-
r
3
)
<
a
11
r
2
r
3.
Furthermore, the other two eigenvalues of (8) can be determined by the following equation:
(9)
λ
2
+
(
r
3
+
β
+
2
a
11
x
1
*
)
λ
+
r
3
(
β
+
2
a
11
x
1
*
)
-
β
r
1
=
0
.
It follows from Routh-Hurwitz criterion [29] that Re
λ
2
<
0, Re
λ
3
<
0 provided that 2
a
11
a
12
r
2
r
3
>
β
(
r
1
-
r
3
)
(
a
12
a
21
-
a
11
a
22
).
Consequently,E
1 is a saddle node of model (5) provided that E
* exists; E
1 is locally stable if 2
a
11
a
12
r
2
r
3
>
β
(
r
1
-
r
3
)
(
a
12
a
21
-
a
11
a
22
) and a
21
β
(
r
1
-
r
3
)
<
a
11
r
2
r
3.Define new variablesu
(
t
) and v
(
t
) as follows:
(10)
u
(
t
)
=
∫
-
∞
t
α
2
(
t
-
s
)
e
-
α
(
t
-
s
)
x
1
(
s
)
d
s
,
v
(
t
)
=
∫
-
∞
t
α
e
-
α
(
t
-
s
)
x
1
(
s
)
d
s
.According to the law of solving the derivative for an integral with parameterized variables, it can be obtained that(11)
u
˙
(
t
)
=
α
v
(
t
)
-
α
u
(
t
)
,
v
˙
(
t
)
=
α
2
x
1
(
t
)
-
α
2
v
(
t
)
.By virtue of (11), model (5) can be transformed into the following five-dimensional system of FDEs with a discrete delay:
(12)
x
˙
1
(
t
)
=
r
1
x
2
(
t
)
-
a
11
x
1
2
(
t
)
-
a
12
x
1
(
t
)
y
(
t
-
τ
)
-
β
x
1
(
t
)
,
x
˙
2
(
t
)
=
β
x
1
(
t
)
-
r
3
x
2
(
t
)
,
y
˙
(
t
)
=
y
(
t
)
[
-
r
2
+
a
21
u
(
t
)
-
a
22
y
(
t
)
]
,
u
˙
(
t
)
=
α
v
(
t
)
-
α
u
(
t
)
,
v
˙
(
t
)
=
α
2
x
1
(
t
)
-
α
2
v
(
t
)
.It follows from model (12) that u
*
=
v
*
=
x
1
* for the interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
,
u
*
,
v
*
) of model (12), which will be denoted as E
*
(
x
1
*
,
x
2
*
,
y
*
,
x
1
*
,
x
1
*
) in the following part of this paper. Furthermore, the interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
) of model (5) is transformed into the interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
,
x
1
*
,
x
1
*
) of model (12). Hence, the qualitative local stability analysis around E
*
(
x
1
*
,
x
2
*
,
y
*
,
x
1
*
,
x
1
*
) of model (12) is equivalent to the qualitative local stability analysis around interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
) of model (5).According to the leading matrix of model (12) around the interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
,
x
1
*
,
x
1
*
), the characteristic equation takes the following form:
(13)
|
λ
+
2
a
11
x
1
*
+
a
12
y
*
+
β
-
r
1
a
12
x
1
*
e
-
λ
τ
0
0
-
β
λ
+
r
3
0
0
0
0
0
λ
+
a
22
y
*
-
a
21
y
*
0
0
0
0
λ
+
α
-
α
-
α
2
0
0
0
λ
+
α
2
|
=
0
,
which derives that
(14)
M
(
λ
)
+
N
(
λ
)
e
-
λ
τ
=
0
,
where
(15)
M
(
λ
)
=
λ
5
+
m
1
λ
4
+
m
2
λ
3
+
m
3
λ
2
+
m
4
λ
+
m
5
,
N
(
λ
)
=
n
4
λ
+
n
5
,
m
1
=
2
a
11
x
1
*
+
(
a
12
+
a
22
)
y
*
+
β
+
α
+
α
2
+
r
3
,
m
2
=
r
3
(
2
a
11
x
1
*
+
a
12
y
*
)
+
β
(
r
3
-
r
1
)
+
α
3
+
a
22
y
*
(
α
+
α
2
)
+
(
a
22
y
*
+
α
+
α
2
)
×
(
2
a
11
x
1
*
+
a
12
y
*
+
β
+
r
3
)
,
m
3
=
(
a
22
y
*
+
α
+
α
2
)
×
[
r
3
(
2
a
11
x
1
*
+
a
12
y
*
)
+
β
(
r
3
-
r
1
)
]
+
(
2
a
11
x
1
*
+
a
12
y
*
+
β
+
r
3
)
×
[
α
3
+
a
22
y
*
(
α
+
α
2
)
]
+
a
22
α
3
y
*
,
m
4
=
[
r
3
(
2
a
11
x
1
*
+
a
12
y
*
)
+
β
(
r
3
-
r
1
)
]
[
α
3
+
a
22
y
*
(
α
+
α
2
)
]
+
a
22
α
3
y
*
(
2
a
11
x
1
*
+
a
12
y
*
+
β
+
r
3
)
,
m
5
=
a
22
α
3
y
*
[
r
3
(
2
a
11
x
1
*
+
a
12
y
*
)
+
β
(
r
3
-
r
1
)
]
,
n
4
=
a
12
a
21
α
3
x
1
*
y
*
,
n
5
=
a
12
a
21
r
3
α
3
x
1
*
y
*
.Whenτ
=
0, (14) takes the following form:
(16)
λ
5
+
m
1
λ
4
+
m
2
λ
3
+
m
3
λ
2
+
(
m
4
+
n
4
)
λ
+
m
5
+
n
5
=
0
.Based on the values ofm
k, n
j, k
=
1, 2, 3, 4, 5, j
=
4, 5 defined in (14), it is easy to show that m
1
>
0 and m
5
+
n
5
>
0.Define quantitiesδ
k, k
=
1,2
,
3,4
,
5 as follows:
(17)
δ
k
=
|
m
1
1
0
0
0
m
3
m
2
m
1
1
0
m
5
+
n
5
m
4
+
n
4
m
3
m
2
m
1
0
0
m
5
+
n
5
m
4
+
n
4
m
3
0
0
0
0
m
5
+
n
5
|
=
a
k
δ
k
-
1
.Further computations show thatδ
k
>
0 for k
=
2,3
,
4. According to the above analysis and Routh-Hurwitz criteria [2], it follows from m
1
>
0, m
5
+
n
5
>
0, and δ
k
>
0 for k
=
2,3
,
4 that model (12) is locally stable around E
*
(
x
1
*
,
x
2
*
,
y
*
,
x
1
*
,
x
1
*
) in the case of τ
=
0, which implies that model (5) is locally stable around E
*
(
x
1
*
,
x
2
*
,
y
*
).Subsequently, the dynamic effect of discrete delay on model dynamics is investigated. In the case of positive discrete delayτ
>
0, let i
ω (ω
>
0) be a root of (14). Separating the real and imaginary parts of (14) with λ
=
i
ω yields the following equations:
(18)
ω
5
+
m
2
ω
3
+
m
4
ω
=
n
5
sin
⁡
ω
τ
-
n
4
ω
cos
⁡
⁡
ω
τ
,
m
3
ω
2
-
m
1
ω
4
-
m
5
=
n
5
cos
⁡
⁡
ω
τ
+
n
4
ω
sin
⁡
ω
τ
.By adding up the squares of the corresponding sides of (18), it yields the following algebraic equation with respect to ω:
(19)
ω
10
+
A
1
ω
8
+
A
2
ω
6
+
A
3
ω
4
+
A
4
ω
2
+
A
5
=
0
,
where A
1
=
2
m
2
+
m
1
2, A
2
=
m
2
2
+
2
m
4
-
2
m
1
m
3, A
3
=
m
3
2
+
2
(
m
2
m
4
+
m
1
m
5
), A
4
=
m
4
2
-
2
m
3
m
5
-
n
4
2, and A
5
=
m
5
2
-
n
5
2.Letξ
=
ω
2; (19) can be rewritten as follows:
(20)
g
(
ξ
)
∶
=
ξ
5
+
A
1
ξ
4
+
A
2
ξ
3
+
A
3
ξ
2
+
A
4
ξ
+
A
5
=
0
.By virtue of derivative ofg
(
ξ
) with respect to ξ, it gives that
(21)
g
′
(
ξ
)
=
5
ξ
4
+
4
A
1
ξ
3
+
3
A
2
ξ
2
+
2
A
3
ξ
+
A
4
.By simple computation, it is easy to show thatA
k
>
0, k
=
1, 2, 3, 4, which derives that g
′
(
ξ
)
>
0 on (
0
,
+
∞
) and g
(
ξ
) is strictly monotonically increasing on (
0
,
+
∞
). Consequently, g
(
ξ
) has a unique positive root if A
5
<
0 and g
(
ξ
) does not have a positive root if A
5
≥
0.SupposeA
5
<
0, which derives that
(22)
(
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
>
0
,
and then there exists a unique root of g
(
ξ
)
=
0, which is denoted by ξ
0, and the unique positive root of (14) is ω
0
=
ξ
0.SupposeA
5
≥
0, which derives that
(23)
(
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
≤
0
,
and then (14) does not have a positive root.By eliminatingsin
⁡
ω
τ, it follows from (18) that
(24)
τ
k
=
1
ω
0
[
+
2
k
π
arccos
(
(
n
4
2
ω
0
2
+
n
5
2
)
-
1
(
-
n
4
ω
0
6
-
(
m
1
n
5
+
m
2
n
4
)
ω
0
4
kkkkkkkkkkkkk
+
(
m
3
n
5
-
m
4
n
4
)
ω
0
2
-
m
5
n
5
)
h
h
h
h
h
h
h
h
k
i
h
h
×
(
n
4
2
ω
0
2
+
n
5
2
)
-
1
)
kkkkk
+
2
k
π
(
(
n
4
2
ω
0
2
+
n
5
2
)
-
1
(
-
n
4
ω
0
6
-
(
m
1
n
5
+
m
2
n
4
)
ω
0
4
]
,
and then (14) has a pair of purely imaginary roots ±
i
ω
0 when τ
=
τ
k for k
=
0,1
,
….By using Butler’s lemma [29], the interior equilibrium E
* remains stable for 0
≤
τ
<
τ
0. Based on the above analysis, local stability analysis of model (5) around the interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
) can be concluded in the following Theorem 2.Theorem 2.
If condition (6) holds, then the following results hold.(i)
If(
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
≤
0, then model (5) is locally stable around E
*
(
x
1
*
,
x
2
*
,
y
*
) for any values τ
≥
0.
(ii)
If(
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
>
0, then model (5) is locally stable around E
*
(
x
1
*
,
x
2
*
,
y
*
) for τ
∈
[
0
,
τ
0
).
Subsequently, the conditions for Hopf bifurcation in [29] are utilized to investigate whether there is a phenomenon of Hopf bifurcation as τ increases through τ
k.Theorem 3.
If(
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
>
0, then model (5) undergoes Hopf bifurcation at the interior equilibrium E
* when τ
=
τ
k, k
=
0,1
,
2
,
…. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium E
* when τ
>
τ
0 and ∥
τ
-
τ
0
∥
≪
1.Proof.
As mentioned above, letλ
=
i
ω
0 represent the purely imaginary root of (14). It follows from (14) that |
M
(
i
ω
0
)
|
=
|
N
(
i
ω
0
)
|, which determines a set of possible values of ω
0.
In the following part, we determine the direction of motion ofλ
=
i
ω
0 as τ is varied; namely, we determine
(25)
Θ
=
sign
⁡
[
d
(
Re
λ
)
d
τ
]
λ
=
i
ω
0
=
sign
⁡
[
Re
(
d
λ
d
τ
)
-
1
]
λ
=
i
ω
0
.
By differentiating (14) with respect to τ, it can be obtained that
(26)
[
5
λ
4
+
4
m
1
λ
3
+
3
m
2
λ
2
+
2
m
3
λ
+
m
4
+
n
4
e
-
λ
τ
k
k
-
τ
(
n
4
λ
+
n
5
)
e
-
λ
τ
]
d
λ
d
τ
=
λ
(
n
4
λ
+
n
5
)
e
-
λ
τ
,
⟹
(
d
λ
d
τ
)
-
1
=
5
λ
4
+
4
m
1
λ
3
+
3
m
2
λ
2
+
2
m
3
λ
+
m
4
λ
(
n
4
λ
+
n
5
)
e
-
λ
τ
+
n
4
λ
(
n
4
λ
+
n
5
)
-
τ
λ
=
5
λ
4
+
4
m
1
λ
3
+
3
m
2
λ
2
+
2
m
3
λ
+
m
4
-
λ
(
λ
5
+
m
1
λ
4
+
m
2
λ
3
+
m
3
λ
2
+
m
4
λ
+
m
5
)
+
n
4
λ
(
n
4
λ
+
n
5
)
-
τ
λ
=
4
λ
5
+
3
m
1
λ
4
+
2
m
2
λ
3
+
m
3
λ
2
-
m
5
-
λ
2

(
λ
5
+
m
1
λ
4
+
m
2
λ
3
+
m
3
λ
2
+
m
4
λ
+
m
5
)
-
n
5
λ
2
(
n
4
λ
+
n
5
)
-
τ
λ
,
Θ
=
sign
⁡
[
Re
(
d
λ
d
τ
)
-
1
]
λ
=
i
ω
0
=
1
ω
0
2
sign
⁡
[
Re
(
n
5
n
4
ω
0
i
+
n
5
(
3
m
1
ω
0
4
-
m
3
ω
0
2
-
m
5
k
k
k
k
k
k
k
k
k
k
k
k
h
+
2
ω
0
3
(
2
ω
0
2
-
m
2
)
i
)
k
k
k
k
k
k
k
k
k
k
k
k
k
×
(
m
1
ω
0
4
-
m
3
ω
0
2
+
m
5
k
k
k
k
k
h
k
k
k
k
k
k
k
+
(
ω
0
5
-
m
2
ω
0
3
+
m
4
ω
0
)
i
)
-
1
k
k
k
k
k
k
k
k
k
k
k
k
+
n
5
n
4
ω
0
i
+
n
5
+
ω
0
τ
i
)
]
=
1
ω
0
2
sign
⁡
[
+
(
ω
0
5
-
m
2
ω
0
3
+
m
4
ω
0
)
2
)
-
1
(
3
m
1
2
ω
0
8
+
4
m
1
m
3
ω
0
6
k
k
k
k
k
k
k
k
k
k
+
(
2
m
1
m
5
+
m
3
2
)
ω
0
4
+
n
5
2
-
m
5
2
)
k
k
k
k
k
k
k
k
k
×
(
+
(
ω
0
5
-
m
2
ω
0
3
+
m
4
ω
0
)
2
(
m
1
ω
0
4
-
m
3
ω
0
2
+
m
5
)
2
k
k
k
k
k
k
k
k
k
k
k
+
(
ω
0
5
-
m
2
ω
0
3
+
m
4
ω
0
)
2
)
-
1
]
.
If(
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
>
0, then n
5
-
m
5
>
0 which signifies that
(27)
Θ
=
sign
⁡
[
Re
(
d
λ
d
τ
)
-
1
]
λ
=
i
ω
0
>
0
.
Hence, the transversality condition [29] holds and model (5) undergoes Hopf bifurcation at the interior equilibrium E
* when τ
=
τ
k, k
=
0,1
,
2
,
…. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium E
* when τ
>
τ
0 and ∥
τ
-
τ
0
∥
≪
1.


## 3.2. Direction and Stability of Hopf Bifurcation

By using normal theory and center manifold theorem [29, 30], directions of Hopf bifurcation and stability of the bifurcating periodic solutions are discussed in this section. As analyzed in Section 3.1, some transformations associated with the interior equilibrium E
*
(
x
1
*
,
x
2
*
,
y
*
,
x
1
*
,
x
1
*
) of model (12) are given as follows:
(28)
z
1
=
x
1
-
x
1
*
,
z
2
=
x
2
-
x
2
*
,
z
3
=
y
-
y
*
,
z
4
=
u
-
x
1
*
,
z
5
=
v
-
x
1
*
,
z
¯
i
(
t
)
=
z
i
(
τ
t
)
,
τ
=
μ
+
τ
k
,
μ
∈
ℝ
=
(
-
∞
,
+
∞
)
,
and then μ
=
0 is the Hopf bifurcation value of model (12). Bars of variables are dropped for simplicity of notations; model (12) is transformed to a functional differential equation in C
=
C
(
[
-
1,0
]
,
ℝ
3
),
(29)
z
˙
(
t
)
=
L
μ
(
z
t
)
+
f
(
μ
,
z
t
)
,
where C
=
C
(
[
-
1,0
]
,
ℝ
5
) is the Banach space of continuous functions mapping [
-
τ
,
0
] into ℝ
5, z
(
t
)
=
(
z
1
(
t
)
,
z
2
(
t
)
,
z
3
(
t
)
,
z
4
(
t
)
,
z
5
(
t
)
)
T
∈
ℝ
5, z
t
(
θ
)
=
z
(
t
+
θ
) for θ
∈
[
-
τ
,
0
], and L
μ
:
C
→
ℝ
5, f
:
ℝ
×
C
→
ℝ
5 are defined as follows, respectively:
(30)
L
μ
(
ϕ
)
=
(
τ
k
+
μ
)
×
(
-
a
11
x
1
*
-
r
1
x
2
*
x
1
*
r
1
0
0
0
β
-
r
3
0
0
0
0
0
-
a
22
y
*
a
21
y
*
0
0
0
0
-
α
α
α
2
0
0
0
-
α
2
)
×
(
ϕ
1
(
0
)
ϕ
2
(
0
)
ϕ
3
(
0
)
ϕ
4
(
0
)
ϕ
5
(
0
)
)
+
(
τ
k
+
μ
)
(
0
0
-
a
12
x
1
*
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
)
(
ϕ
1
(
-
1
)
ϕ
2
(
-
1
)
ϕ
3
(
-
1
)
ϕ
4
(
-
1
)
ϕ
5
(
-
1
)
)
,
(31)
f
(
μ
,
ϕ
)
=
(
τ
k
+
μ
)
(
-
a
12
ϕ
1
(
0
)
ϕ
3
(
-
1
)
0
0
0
0
)
.It is easy to show thatL
μ is a continuous linear function mapping C into ℝ
5. According to Riesz representation theorem [4], there exists a 5
×
5 matrix function η
(
θ
,
μ
) of bounded variation for θ
∈
[
-
1,0
] such that
(32)
L
μ
(
ϕ
)
=
∫
-
1
0
d
η
(
θ
,
μ
)
ϕ
(
θ
)
,
where ϕ
∈
C
(
[
-
1,0
]
,
ℝ
5
).In fact, we can choose(33)
η
(
θ
,
μ
)
=
(
τ
k
+
μ
)
×
(
-
a
11
x
1
*
-
r
1
x
2
*
x
1
*
r
1
0
0
0
β
-
r
3
0
0
0
0
0
-
a
22
y
*
a
21
y
*
0
0
0
0
-
α
α
α
2
0
0
0
-
α
2
)
×
δ
(
θ
)
-
(
τ
k
+
μ
)
(
0
0
-
a
12
x
1
*
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
)
×
δ
(
θ
+
1
)
,
where δ denotes the Dirac delta function.Ifϕ is any given function in C
(
[
-
1,0
]
,
ℝ
5
) and z
(
ϕ
) is the unique solution of the linearized equation z
˙
(
t
)
=
L
μ
(
z
t
) of (29) with initial function ϕ at zero, then the solution operator T
~
(
t
)
:
C
→
C is defined by
(34)
T
~
(
t
)
ϕ
=
z
t
(
ϕ
)
,
t
≥
0
.It follows from Lemma  7.1.1 in [29] that T
~
(
t
), t
≥
0 is a strongly continuous semigroup of linear transformation on [
0
,
+
∞
) and the infinitesimal generator A
μ of T
~
(
t
), t
≥
0 is as follows:
(35)
A
μ
(
ϕ
)
=
{
d
ϕ
(
θ
)
d
θ
,
θ
∈
[
-
1,0
)
∫
-
1
0
d
η
(
μ
,
s
)
ϕ
(
s
)
,
θ
=
0
,
for ϕ
∈
C
1
(
[
-
1,0
]
,
ℝ
5
), the space of functions mapping the interval [
-
1,0
] into ℝ
5 which have continuous first derivative, and also define
(36)
R
(
μ
)
(
ϕ
)
=
{
0
,
θ
∈
[
-
1,0
)
f
(
μ
,
ϕ
)
,
θ
=
0
,
and then model (29) is equivalent to
(37)
z
˙
t
=
A
(
μ
)
z
t
+
R
(
μ
)
z
t
.Forψ
∈
C
1
(
[
0,1
]
,
(
ℝ
5
)
*
), the space of functions mapping interval [
0,1
] into the five-dimensional row vectors which have continuous first derivative, define
(38)
A
*
ϕ
(
s
)
=
{
-
d
ψ
(
s
)
d
s
,
s
∈
(
0,1
]
∫
-
1
0
d
η
T
(
t
,
0
)
ψ
(
-
t
)
,
s
=
0
,
and a bilinear inner product
(39)
〈
ψ
(
s
)
,
ϕ
(
θ
)
〉
=
ψ
¯
(
0
)
ϕ
(
0
)
-
∫
-
1
0
∫
ρ
=
0
θ
ψ
¯
(
ρ
-
θ
)
d
η
(
θ
)
ϕ
(
ρ
)
d
ρ
,
where η
(
θ
)
=
η
(
θ
,
0
). It follows from the above analysis that A
(
0
) and A
* are adjoint operators.By virtue of discussion in Section3.1, ±
i
ω
0
τ
k are eigenvalues of A
(
0
). Hence, they are also eigenvalues of A
*. In the following, eigenvectors of A
(
0
) and A
* correspond to i
ω
0
τ
k and -
i
ω
0
τ
k, respectively.Supposeq
(
θ
)
=
(
1
,
a
,
b
,
c
,
d
)
T
e
i
ω
0
τ
k
θ is the eigenvectors of A
(
0
) corresponding to i
ω
0
τ
k, which derives that A
(
0
)
q
(
θ
)
=
i
ω
0
τ
k
q
(
θ
). By using the definition of A
(
0
), (30), and (31), it gives that
(40)
τ
k
(
-
a
11
x
1
*
-
r
1
x
2
*
x
1
*
r
1
0
0
0
β
-
r
3
0
0
0
0
0
-
a
22
y
*
a
21
y
*
0
0
0
0
-
α
α
α
2
0
0
0
-
α
2
)
q
(
0
)
+
τ
k
(
0
0
-
a
12
x
1
*
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
)
q
(
-
1
)
=
i
ω
0
τ
k
q
(
0
)
.Forq
(
-
1
)
=
q
(
0
)
e
-
i
ω
0
τ
k, then it can be obtained that
(41)
a
=
β
r
3
+
i
ω
0
,
b
=
x
1
*
[
β
r
1
-
(
a
11
x
1
*
+
i
ω
0
)
(
r
3
+
i
ω
0
)
]
-
r
1
x
2
*
(
r
3
+
i
ω
0
)
a
12
x
1
*
2
e
-
i
ω
0
τ
k
(
r
3
+
i
ω
0
)
,
c
=
α
3
(
α
+
i
ω
0
)
(
α
2
+
i
ω
0
)
,
d
=
α
2
α
2
+
i
ω
0
.Similarly, it follows from simple computation that eigenvector ofA
* corresponds to -
i
ω
0
τ
k; that is, q
*
(
s
)
=
J
(
1
,
a
*
,
b
*
,
c
*
,
d
*
)
T
e
i
ω
0
τ
k
s
(42)
a
*
=
r
1
r
3
+
i
ω
0
,
b
*
=
-
a
12
x
1
*
e
-
i
ω
0
τ
k
a
22
y
*
+
i
ω
0
,
c
*
=
-
a
12
a
21
x
1
*
e
-
i
ω
0
τ
k
(
α
+
i
ω
0
)
(
a
22
y
*
+
i
ω
0
)
,
d
*
=
-
α
a
12
a
21
x
1
*
e
-
i
ω
0
τ
k
(
α
+
i
ω
0
)
(
α
2
+
i
ω
0
)
(
a
22
y
*
+
i
ω
0
)
.In order to assume〈
q
*
(
s
)
,
q
(
θ
)
〉
=
1, we need to determine the value of J in the following part.By virtue of(A.9), it derives that
(43)
〈
q
*
(
s
)
,
q
(
θ
)
〉
=
J
¯
(
1
,
a
¯
*
,
b
¯
*
,
c
¯
*
,
d
¯
*
)
(
1
,
a
,
b
,
c
,
d
)
T
-
∫
-
1
0
∫
ρ
=
0
θ
J
¯
(
1
,
a
*
,
b
*
,
c
*
,
d
*
)
e
-
i
ω
0
τ
k
(
ρ
-
θ
)
d
η
(
θ
)
k
k
k
k
k
k
k
k
h
k
k
×
(
1
,
a
,
b
,
c
,
d
)
T
e
i
ω
0
τ
k
ρ
d
ρ
=
J
¯
(
1
+
a
a
¯
*
+
b
b
¯
*
+
c
c
¯
*
+
d
d
¯
*
)
-
J
¯
∫
-
1
0
(
1
,
a
¯
*
,
b
¯
*
,
c
¯
*
,
d
¯
*
)
θ
e
i
ω
0
τ
k
θ
d
η
(
θ
)
(
1
,
a
,
b
,
c
,
d
)
T
=
J
¯
(
1
+
a
a
¯
*
+
b
b
¯
*
+
c
c
¯
*
+
d
d
¯
*
+
B
)
,
where
(44)
B
=
[
r
1
x
2
*
x
1
*
α
2
d
¯
*
+
β
a
¯
*
+
a
(
r
1
-
a
¯
*
(
r
3
+
i
ω
0
)
)
k
-
a
11
x
1
*
-
i
ω
0
-
r
1
x
2
*
x
1
*
]
-
b
[
a
12
x
1
*
e
-
i
ω
0
τ
k
+
b
¯
*
(
a
22
y
*
+
i
ω
0
)
]
+
c
[
a
21
b
¯
*

y
*
-
c
¯
*
(
α
+
i
ω
0
)
]
+
d
[
α
c
¯
*
-
(
α
2
+
i
ω
0
)
d
¯
*
]
.Hence, we can chooseJ as follows:
(45)
J
=
1
1
+
a
a
¯
*
+
b
b
¯
*
+
c
c
¯
*
+
d
d
¯
*
+
B
.Nextly, we will compute the coordinate to describe the centre manifoldC
0 at μ
=
0. Let z
t be the solution of (37) when μ
=
0.Define(46)
χ
(
t
)
=
〈
q
*
,
z
t
〉
,
W
(
t
,
θ
)
=
z
t
(
θ
)
-
2
Re
{
χ
(
t
)
q
(
θ
)
}
.On the center manifoldC
0, it derives that
(47)
W
(
t
,
θ
)
=
W
(
χ
(
t
)
,
χ
¯
(
t
)
,
θ
)
,
where
(48)
W
(
χ
(
t
)
,
χ
¯
(
t
)
,
θ
)
=
W
20
(
θ
)
χ
2
2
+
W
11
(
θ
)
χ
χ
¯
+
W
02
χ
¯
2
2
+
⋯
,
χ and χ
¯ are local coordinates for center manifold C
0 in the direction of q
* and q
¯
*.It is noted thatW is real if z
t is real, and we only consider real solutions. For solution z
t
∈
C
0 of (37), since μ
=
0, it derives that
(49)
χ
˙
(
t
)
=
i
ω
0
τ
k
χ
+
q
¯
*
(
0
)
f
(
0
,
W
(
χ
,
χ
¯
,
0
)
+
2
Re
{
χ
q
(
θ
)
}
)
≜
i
ω
0
τ
k
χ
+
q
¯
*
(
0
)
f
0
(
χ
,
χ
¯
)
.The above equation can be rewritten as follows:(50)
χ
˙
(
t
)
=
i
ω
0
τ
k
χ
(
t
)
+
g
(
χ
,
χ
¯
)
,
where
(51)
g
(
χ
,
χ
¯
)
=
q
¯
*
(
0
)
f
0
(
χ
,
χ
¯
)
=
g
20
χ
2
2
+
g
11
χ
χ
¯
+
g
02
χ
¯
2
2
+
g
21
χ
2
χ
¯
2
+
⋯
.It follows from (46) and (48) that
(52)
z
t
(
θ
)
=
W
(
t
,
θ
)
+
2
Re
{
χ
(
t
)
q
(
θ
)
}
=
W
20
(
θ
)
χ
2
2
+
W
11
(
θ
)
χ
χ
¯
+
W
02
(
θ
)
χ
¯
2
2
+
(
1
,
a
,
b
,
c
,
d
)
T
e
i
ω
0
τ
k
θ
χ
+
(
1
,
a
¯
,
b
¯
,
c
¯
,
d
¯
)
T
e
-
i
ω
0
τ
k
θ
χ
¯
+
⋯
.By virtue of (31), (32), and (33), it derives that
(53)
g
(
χ
,
χ
¯
)
=
q
¯
*
(
0
)
f
0
(
χ
,
χ
¯
)
=
q
¯
*
(
0
)
f
(
0
,
z
t
)
=
τ
k
J
¯
(
1
,
a
¯
*
,
b
¯
*
,
c
¯
*
,
d
¯
*
)
(
-
a
12
z
1
t
(
0
)
z
3
t
(
-
1
)
0
0
0
0
)
=
-
τ
k
J
¯
a
12
z
1
t
(
0
)
z
3
t
(
-
1
)
=
-
τ
k
J
¯
a
12
[
χ
+
χ
¯
+
W
20
(
1
)
(
0
)
χ
2
2
+
W
11
(
1
)
(
0
)
χ
χ
¯
kkkkkkkkk
+
W
02
(
1
)
(
0
)
χ
¯
2
2
+
o
(
|
(
χ
,
χ
¯
)
|
3
)
]
×
[
χ
+
χ
¯
+
W
20
(
3
)
(
-
1
)
χ
2
2
+
W
11
(
3
)
(
-
1
)
χ
χ
¯
k
k
k
k
k
+
W
02
(
3
)
(
-
1
)
χ
¯
2
2
+
o
(
|
(
χ
,
χ
¯
)
|
3
)
]
=
-
τ
k
J
¯
a
12
[
+
W
20
(
1
)
(
0
)
+
W
20
(
3
)
(
-
1
)
2
)
χ
2
+
2
χ
χ
¯
+
χ
¯
2
k
k
k
k
h
k
k
k
k
+
(
W
20
(
1
)
(
0
)
+
W
20
(
3
)
(
-
1
)
2
W
11
(
3
)
(
-
1
)
+
W
11
(
1
)
(
0
)
k
k
k
k
k
k
k
k
k
k
k
+
W
20
(
1
)
(
0
)
+
W
20
(
3
)
(
-
1
)
2
)
χ
2
χ
¯
+
⋯
]
.By comparing the coefficients with (51), it gives that
(54)
g
20
=
-
2
a
12
τ
k
J
¯
,
g
11
=
-
2
a
12
τ
k
J
¯
,
g
02
=
-
2
a
12
τ
k
J
¯
,
g
21
=
-
2
a
12
τ
k
J
¯
×
(
W
11
(
3
)
(
-
1
)
+
W
11
(
1
)
(
0
)
+
W
20
(
1
)
(
0
)
+
W
20
(
3
)
(
-
1
)
2
)
.Sinceg
21 is associated with W
20
(
θ
) and W
11
(
θ
), further attempts should be carried out to compute W
20
(
θ
) and W
11
(
θ
), which can be found in the appendix.Furthermore, we can computeg
21 based on (54). Hence, the following values can be computed as follows:
(55)
d
1
(
0
)
=
i
2
ω
0
τ
k
(
g
20
g
11
-
2
|
g
11
|
2
-
|
g
02
|
2
3
)
+
g
21
2
,
δ
2
=
-
Re
{
d
1
(
0
)
}
Re
{
λ
′
(
τ
k
)
}
,
γ
2
=
2
Re
(
d
1
(
0
)
)
,
T
2
=
Im
⁡
{
d
1
(
0
)
}
+
δ
2
Im
⁡
{
λ
′
(
τ
k
)
}
ω
0
τ
k
.Some properties of bifurcating periodic solution in the center manifold at the critical valuesτ
k follow from [29]. Based on the analysis in Section 3.1 of this paper, the following theorem can be concluded.Theorem 4.
The properties of Hopf bifurcation are determined by values computed in (55).(i)
μ
2 determines directions of Hopf bifurcation: if μ
2
>
0  (
μ
2
<
0
), then Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for τ
>
τ
k
(
τ
<
τ
k
).
(ii)
γ
2 determines the stability of bifurcating periodic solutions: bifurcating periodic solutions are stable (unstable) if γ
2
<
0
(
γ
2
>
0
).
(iii)
T
2 determines the period of bifurcating periodic solutions: period increases (decreases) if T
2
>
0
(
T
2
<
0
).


## 4. Numerical Simulation

In this section, values of parameters are partially taken from numerical simulation of [13] and set in appropriate units, which can be found as follows: r
1
=
7.889, a
11
=
0.7, a
12
=
0.7, β
=
2.657, r
3
=
0.9, a
21
=
2, r
2
=
1, a
22
=
0.8, and α
=
2. Numerical simulations are provided to support the theoretical findings obtained in Section 3 of this paper.By virtue of given values of parameters, it can be computed that(56)
a
21
β
(
r
1
-
r
3
)
-
a
11
r
2
r
3
=
36.5095
>
0
,
which implies that (6) holds and there exists a unique interior equilibrium of model (5), that is, E
*
(
5.058,14.932,11.395
). Furthermore, it follows from given values of parameters that
(57)
(
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
=
13.0081
>
0
,
which shows that (22) holds and (19) has a unique positive root. Based on (24), the corresponding τ
0
=
2.4576 can be computed. According to Theorem 2, model is locally stable around the interior equilibrium E
*
(
5.058,14.932,11.395
) when τ
∈
[
0,2.4576
). Dynamical responses of model (5) with τ
=
1.26 are plotted in Figure 1 and phase portrait of model (5) with τ
=
1.26 is plotted in Figure 2. It should be noted that τ
=
1.26 in Figures 1 and 2 is randomly selected in the interval [
0,2.4576
), which is enough to merit the above mathematical study. It follows from (A.18) that μ
2
=
3.9152
>
0, γ
2
=
-
0.7831
<
0, and T
2
=
1.9724
>
0. As τ increases through τ
0, a periodic solution caused by the phenomenon of Hopf bifurcation occurs, that is, a family of periodic solutions bifurcate from the interior equilibrium E
* based on Theorem 3. Since μ
2
>
0 and γ
2
<
0, the Hopf bifurcation is supercritical, the directions of the Hopf bifurcation are τ
>
τ
0, and these bifurcating periodic solutions from the interior equilibrium E
* at τ
0 are stable based on Theorem 4. Dynamical responses of model (5) with τ
=
2.5 are plotted in Figure 3, since τ
=
2.5
>
τ
0, which shows that model (5) is unstable around the interior equilibrium and stable periodic solutions bifurcate from the interior equilibrium E
*
(
5.058,14.932,11.395
). Furthermore, Figure 4 indicates a limit cycle corresponding to the stable periodic solutions plotted in Figure 3.Figure 1
Dynamical responses of model (5) with τ
=
1.26, which shows that model (5) is stable around the interior equilibrium (
5.058,14.932,11.395
).Figure 2
Phase portrait of model (5) with τ
=
1.26 corresponding to stable dynamical responses in Figure 1, which shows that model (5) is asymptotically stable around the interior equilibrium (
5.058,14.932,11.395
).Figure 3
Dynamical responses of model (5) with τ
=
2.5, which shows that model (5) is unstable around the interior equilibrium and stable periodic solutions bifurcate from the interior equilibrium (
5.058,14.932,11.395
).Figure 4
A limit cycle corresponding to the stable periodic solutions plotted in Figure3.


## 5. Conclusion

In this paper, a dynamical prey predator model with discrete hunting delay and distributed maturation delay is proposed to investigate the dynamic effect of time delay and stage structure on population dynamics. Many species in the natural world have a life history that takes them through two stages, juvenile stage and adult stage. Individuals in each stage are identical in biological characteristics, and some vital rates (rates of survival, development, and reproduction) of individuals in a population almost always depend on stage structure. In model (4) with strong generic kernel [13], stage structure of prey population is not considered. It is well known that the immature prey population is usually considered as an easier target for predators compared with mature prey population. Hence, the mature prey population predated by predator population can be ignored [1], and it is necessary to investigate the dynamic effect of stage structure on population dynamics of model (4).By incorporating stage structure for prey population into the model, the work done in [13] is extended in this paper, where discrete delay and distributed delay for predator population described by an integral with a strong delay kernel are also studied in the proposed model. In Section 3 of this paper, conditions for existence of two boundary equilibria and a unique interior equilibrium are discussed. Local stability around two feasible boundary equilibria is analyzed in Theorem 1. By taking the discrete delay as a parameter, local stability analysis around the unique interior equilibrium is investigated due to variation of discrete delay in Theorem 2. It reveals that discrete delay has dynamic effect on population dynamics, which shows that model (5) is locally stable around interior equilibrium for any nonnegative discrete delay provided that (
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
≤
0 and model (5) remains locally stable around interior equilibrium within τ
∈
[
0
,
τ
0
) under the assumption (
a
12
a
21
-
a
11
a
22
)
[
a
12
r
2
r
3
+
a
22
β
(
r
1
-
r
3
)
]
>
0. According to Hopf bifurcation theorem for functional differential equations, it can be found that Hopf bifurcation occurs when the discrete delay crosses through a sequence of critical values and corresponding stable limit cycle is also observed, which can be found in Theorem 3. Further attempts are made to study the directions of Hopf bifurcation and stability of the bifurcating periodic solutions in Theorem 4.Compared with the previous related work [13], theoretical results obtained in this paper concentrate on stage-structured prey predator model with discrete hunting delay and distributed maturation delay, which are potentially beneficial for administrative agency to have a deep insight on the coexistence and interaction mechanism of prey predator ecosystem in the real world. It makes work done in this paper has some positive and new feature.



---

*Source: 184174-2014-04-27.xml*