# Stability and Hopf Bifurcation Analysis of a Vector-Borne Disease with Time Delay

Authors: Yuanyuan Chen; Ya-Qing Bi
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

A delay-differential modelling of vector-borne is investigated. Its dynamics are studied in terms of local analysis and Hopf bifurcation theory, and its linear stability and Hopf bifurcation are demonstrated by studying the characteristic equation. The stability and direction of Hopf bifurcation are determined by applying the normal form theory and the center manifold argument.


---

## Body

## 1. Introduction

Vector-borne diseases are an important public health problem. Vector-borne diseases are infectious diseases caused by virus, bacteria, and so on which are primarily transmitted by disease biological agents, called vector carrying the disease.Malaria is the most prevalent vector-borne disease, which is transmitted to the human host through a bite by an infected mosquito. It can lead to serous affecting the brain, lungs, kidneys, and other organs, and it caused the greatest number of deaths. Approximately, 40 percent of the world's population is at risk, and 2 million deaths per year can be attributed to malaria, half of those in children under 5 years old. Especially in Africa, more than one million children mostly under 5 years die each year. No effective vaccines are available for the disease. In many years, the effective way to prevent the malaria and other mosquito-borne disease is to control mosquito.Several theoretical studies have proposed vector-borne models. Reference [1] used a mathematical model to show that bringing a mosquito population below a certain threshold was sufficient to eliminate malaria. Reference [2] studied both a baseline ODE version of the model and a model with a discrete time delay and gave the conditions under which equilibrium is globally stable and the disease dies out. Reference [3] showed that reducing the number of mosquitoes is an inefficient control strategy that would have little effect on the epidemiology of malaria in areas of intense transmission. Reference [4] used a mathematical model to evaluate the impact from the programs of selective mass drug administration and vector control through mosquito nets. References [5, 6] models took into account the acquired immunity to malaria depends on exposure (i.e., that immunity is boosted by additional infections).For a long time, it has been recognized that delay may have very complicated impact on the dynamics of a system. Delay can cause the loss of stability and can bifurcate various periodic solutions. Recently, there has been extensive work dealing with time delay systems (see, e.g., [7–11]). As far as we know, there are few works on the delayed vector-borne system, let alone the existence of Hopf bifurcation, and the stability and direction of bifurcating periodic solutions. In this paper, we focus on investigating these problems.This paper is organized as follows. In Section2, we provide a vector-borne model and analyze the property of the nonnegative equilibria. In Section 3, we get the existence of the Hopf bifurcation. In Section 4, the stability and direction of periodic solutions bifurcating from the Hopf bifurcation are determined by using the normal form theory and center manifold argument introduced by Hassard et al. [12].


## 2. Property of the Nonnegative Equilibria

We can describe the dynamics of the disease in the host population as follows:(1)
S
˙
(
t
)
=
b
1
-
λ
1
S
(
t
)
I
(
t
)
-
λ
2
S
(
t
)
V
(
t
)
-
μ
1
S
(
t
)
,
I
˙
(
t
)
=
λ
1
S
(
t
)
I
(
t
)
+
λ
2
S
(
t
)
V
(
t
)
-
γ
I
(
t
)
-
μ
1
I
(
t
)
,
R
˙
(
t
)
=
γ
I
(
t
)
-
μ
1
R
(
t
)
.
Here S
(
t
), I
(
t
), and R
(
t
) represent the population density of susceptible, infectious, and recovered at time t, respectively. The total host population size at time t is given by N
1
(
t
). The host population dies at a natural rate μ
1, and the host grows with intrinsic growth rate b
1. λ
1 is the rate of direct transmission, while λ
2 is the biting rate of a pathogen-carrier vector. The host recovers at the rate γ. The recovered individuals are assumed to acquire permanent immunity.The system that describes the dynamics of the vector is given by(2)
M
˙
(
t
)
=
b
2
-
λ
3
M
(
t
)
I
(
t
)
-
μ
2
M
(
t
)
,
V
˙
(
t
)
=
λ
3
M
(
t
)
I
(
t
)
-
μ
2
V
(
t
)
.
Here, V
(
t
) is the number of vectors at time t carrying the pathogen at time t, and M
(
t
) represents the population density of pathogen-free vector at time t. The total vectors population size at time t is given by N
2
(
t
). b
2 and μ
2 are the birth rate and death rate of vector population, respectively. Suspectable vectors start carrying the pathogen after getting into contact with an infective host at a rate λ
3. We assume that the vectors carry the microparasite for life once they became carrier of it.The time delayτ
>
0 is introduced in the system (2) to describe the dynamics of the vector. At time t, the susceptible vectors bite the host τ time ago, and the vector became infectious. The delay model of the system takes the following form:
(3)
M
˙
(
t
)
=
b
2
-
λ
3
M
(
t
)
I
(
t
-
τ
)
-
μ
2
M
(
t
)
,
V
˙
(
t
)
=
λ
3
M
(
t
)
I
(
t
-
τ
)
-
μ
2
V
(
t
)
.
The systems (1) and (3) satisfy the initial conditions: S
(
θ
)
=
S
0, I
(
θ
)
=
I
0, R
(
θ
)
=
R
0, M
(
θ
)
=
M
0, V
(
θ
)
=
V
0, and θ
∈
[
-
τ
,
0
). The total host population size N
1
(
t
) can be determined by N
1
(
t
)
=
S
(
t
)
+
I
(
t
)
+
R
(
t
) or
(4)
N
˙
1
(
t
)
=
b
1
-
μ
1
N
1
(
t
)
.
The total number of vectors N
2
(
t
) can be determined by N
2
(
t
)
=
M
(
t
)
+
V
(
t
) or
(5)
N
2
(
t
)
=
b
2
-
μ
2
N
2
(
t
)
.
The total population size of both host and vector populations are asymptotically constant; that is, lim
⁡
t
→
∞
N
1
(
t
)
=
b
1
/
μ
1 and lim
⁡
t
→
∞
N
2
(
t
)
=
b
2
/
μ
2. Without loss of generality, we assume that N
1
(
t
)
=
b
1
/
μ
1 and N
2
(
t
)
=
b
2
/
μ
2 for all t
≥
0 provided that S
0
+
I
0
+
R
0
=
b
1
/
μ
1 and M
0
+
V
0
=
b
2
/
μ
2.The systems (1) and (3) are equivalent to the dynamics of the following system:
(6)
S
˙
(
t
)
=
b
1
-
λ
1
S
(
t
)
I
(
t
)
-
λ
2
S
(
t
)
V
(
t
)
-
μ
1
S
(
t
)
,
I
˙
(
t
)
=
λ
1
S
(
t
)
I
(
t
)
+
λ
2
S
(
t
)
V
(
t
)
-
γ
I
(
t
)
-
μ
1
I
(
t
)
,
V
˙
(
t
)
=
λ
3
(
b
2
μ
2
-
V
(
t
)
)
I
(
t
-
τ
)
-
μ
2
V
(
t
)
.The initial condition of system (6) is
(7)
(
S
(
θ
)
,
I
(
θ
)
,
V
(
θ
)
∈
C
+
=
C
(
(
-
τ
,
0
]
,
R
+
3
)
)
,
S
0
,
I
0
,
V
0
>
0
.
System (6) has two equilibria E
10
=
(
b
1
/
μ
1
,
0,0
) and E
*
=
(
S
*
,
I
*
,
V
*
), where
(8)
S
*
=
b
1
-
(
γ
+
μ
1
)
I
*
μ
1
,
V
*
=
λ
3
b
2
I
*
μ
2
(
μ
2
+
λ
3
I
*
)
.
I
* is determined by the following equation:
(9)
b
1
-
(
γ
+
μ
1
)
I
μ
1
(
λ
1
+
λ
2
λ
3
b
2
μ
2
(
λ
3
I
+
μ
2
)
)
=
γ
+
μ
1
.Equation (9) has a unique positive root, when (
γ
+
μ
1
)
μ
1
μ
2
2
<
b
1
(
λ
1
μ
2
2
+
λ
2
λ
3
b
2
).For the equilibriumE
10, the characteristic equation is
(10)
(
μ
1
+
λ
)
(
γ
+
μ
1
-
λ
1
b
1
μ
1
+
λ
)
(
μ
2
+
λ
)
=
0
.
We can easily get the following theorem by some calculation.Theorem 1.
E
10 is asymptotically stable if λ
1
b
1
/
μ
1
-
γ
-
μ
1
<
0; it is unstable if λ
1
b
1
/
μ
1
-
γ
-
μ
1
>
0.


## 3. Existence of Hopf Bifurcation

We studyE
* under the condition (
γ
+
μ
1
)
μ
1
μ
2
2
<
b
1
(
λ
1
μ
2
2
+
λ
2
λ
3
b
2
). The characteristic equation of E
* is
(11)
λ
3
+
(
2
μ
1
+
μ
2
+
γ
+
λ
1
I
*
+
λ
3
I
*
+
λ
2
V
*
-
λ
1
S
*
)
λ
2
+
[
(
μ
2
+
λ
3
I
*
)
(
2
μ
1
+
γ
+
λ
1
I
*
+
λ
2
V
*
-
λ
1
S
*
)
+
μ
1
(
μ
1
+
γ
-
λ
1
S
*
)
+
(
λ
1
I
*
+
λ
2
V
*
)
(
μ
1
+
γ
)
]
λ
+
(
μ
2
+
λ
3
I
*
)
[
μ
1
(
μ
1
+
γ
-
λ
1
S
*
)
+
(
λ
1
I
*
+
λ
2
V
*
)
(
μ
1
+
γ
)
]
+
λ
2
λ
3
S
*
V
*
(
λ
+
μ
1
)
e
-
λ
τ
=
0
.
Let
(12)
A
=
2
μ
1
+
μ
2
+
γ
+
λ
1
I
*
+
λ
3
I
*
+
λ
2
V
*
-
λ
1
S
*
,
(13)
B
=
(
μ
2
+
λ
3
I
*
)
(
2
μ
1
+
γ
+
λ
1
I
*
+
λ
2
V
*
-
λ
1
S
*
)
+
μ
1
(
μ
1
+
γ
-
λ
1
S
*
)
+
(
λ
1
I
*
+
λ
2
V
*
)
(
μ
1
+
γ
)
,
(14)
C
=
μ
1
(
μ
2
-
λ
3
I
*
)
(
μ
1
+
γ
-
λ
1
S
*
)
+
(
λ
1
I
*
+
λ
2
V
*
)
(
μ
2
-
λ
3
I
*
)
,
(15)
D
=
λ
2
λ
3
S
*
V
*
,
E
=
μ
1
λ
2
λ
3
S
*
V
*
.
Then (11) can be rewritten as
(16)
λ
3
+
A
λ
2
+
B
λ
+
C
+
[
D
λ
+
E
]
e
-
λ
τ
=
0
.Lemma 2.
Equation (11) has a unique pair of purely imaginary roots if A
2
-
2
B
>
0 and C
2
-
E
2
<
0.Proof.
Ifλ
=
i
w, w
>
0 is a root of (16), separating real and imaginary parts, we have the following:
(17)
A
w
2
-
C
=
E
cos
⁡
w
τ
+
D
w
sin
w
τ
,
w
3
-
B
w
=
D
w
cos
⁡
w
τ
-
E
sin
w
τ
.
Squaring and adding both equations, we have
(18)
w
6
+
Q
1
w
4
+
Q
2
w
2
+
Q
3
=
0
,
where
(19)
Q
1
=
A
2
-
2
B
,
Q
2
=
B
2
-
2
A
C
-
D
2
,
Q
3

=

C
2
-
E
2
.
We know that
(20)
Q
1
=
A
2
-
2
B
>
0
,
(21)
Q
3
=
C
2
-
E
2
<
0
.
Then this lemma implies that there is a unique positive root w
0 satisfying (16). That is, (11) has a unique pair of purely imaginary roots ±
i
w
0.
From (17), τ
n can be obtained:
(22)
τ
n
=
1
w
0
cos
⁡
-
1
D
w
0
4
+
(
A
E
-
B
D
)
w
0
2
-
E
C
E
2
+
D
2
w
0
2
+
2
n
π
w
0
,
n
=
0,1
,
2
,
…
.Theorem 3.
If the following conditions(23)
A
2
-
2
B
>
0
,
C
2
-
E
2
<
0
,
A
(
B
+
D
)
>
C
+
E
,
(24)
(
B
2
-
2
A
C
)
E
2
>
C
2
D
2
are satisfied, system (3) undergoes Hopf bifurcation at E
* when τ
=
τ
n, n
=
0,1
,
2
,
…; furthermore, E
* is locally asymptotically stable if τ
∈
[
0
,
τ
0
) and unstable if τ
>
τ
0.Proof.
Differentiating (16) with respect to τ, we get
(25)
[
3
λ
2
+
2
A
λ
+
B
+
D
e
-
λ
τ
-
τ
(
D
λ
+
E
)
e
-
λ
τ
]
d
λ
d
τ
=
λ
(
D
λ
+
E
)
e
-
λ
τ
.
That is
(26)
(
d
λ
d
τ
)
-
1
=
-
3
λ
2
+
2
A
λ
+
B
λ
(
λ
3
+
A
λ
2
+
B
λ
+
C
)
+
D
λ
(
D
λ
+
E
)
-
τ
λ
.
Thus,
(27)
Re
(
d
λ
d
τ
)
-
1
|
λ
=
i
w
0
=
Re
(
-
B
-
3
w
0
2
+
i
2
A
w
0
i
w
0
[
(
C
-
A
w
0
2
)

+

i
w
0
(
B
-
w
0
2
)
]

+

D
i
w
0
(
i
D
w
0
+
E
)
)
=
(
2
D
2
w
0
6
+
(
3
E
2
+
A
2
D
2
-
2
B
D
2
)
w
0
4
+
(
2
A
2
-
4
B
)
E
2
w
0
2
+
(
B
2
E
2
-
C
2
D
2
-
2
A
C
E
2
)
)
×
(
[
w
0
2
(
B
-
w
0
2
)
2
+
(
C
-
A
w
0
2
)
2
]
[
(
D
w
0
)
2
+
E
2
]
)
-
1
.
We can rewrite the numerator as follows. Let
(28)
V
=
w
2
,
f
(
w
)
=
2
D
2
w
6
+
(
3
E
2
+
A
2
D
2
-
2
B
D
2
)
w
4
+
(
2
A
2
-
4
B
)
E
2
w
2
+
(
B
2
E
2
-
C
2
D
2
-
2
A
C
E
2
)
,
and then
(29)
G
(
V
)
=
f
(
w
)
=
2
D
2
V
3
+
(
3
E
2
+
A
2
D
2
-
2
B
D
2
)
V
2
+
(
2
A
2
-
4
B
)
E
2
V
+
(
B
2
E
2
-
C
2
D
2
-
2
A
C
E
2
)
,
G
′
=
2
[
3
D
2
V
2
+
(
3
E
2
+
(
A
2
-
2
B
)
D
2
)
V
+
(
A
2
-
2
B
)
E
2
]
.
For G
′,
(30)
Δ
=
(
3
E
2
+
A
2
D
2
-
2
B
D
2
)
2
-
12
(
A
2
-
2
B
)
D
2
E
2
=
[
3
E
2
-
(
A
2
-
2
B
)
D
2
]
2
≥
0
.
G
′ has two real roots, which take the form
(31)
V
1
=
-
(
3
E
2
+
(
A
2
-
2
B
)
D
2
)
+
Δ
6
D
2
<
0
,
V
2
=
-
(
3
E
2
+
(
A
2
-
2
B
)
D
2
)
-
Δ
6
D
2
<
0
.
Then we know that G
(
V
) monotonously increases in (
V
1
,
+
∞
), that is to say, that f
(
w
) monotonously increases in (
0
,
+
∞
). And as we know f
(
0
)
=
B
2
E
2
-
C
2
D
2
-
2
A
C
E
2
>
0, we have f
(
w
)
>
0 for w
>
0. Then we obtain
(32)
sign
⁡
d
(
Re
λ
(
τ
)
)
d
τ
|
τ
=
τ
n
=
signRe
(
d
λ
d
τ
)
-
1
|
λ
=
i
w
0
>
0
.
Therefore, the transversality condition holds and hence Hopf bifurcation occurs. For (16), when τ
=
0, the characteristic equation is
(33)
λ
3
+
A
λ
2
+
(
B
+
D
)
λ
+
C
+
E
=
0
.
Under the conditions of the theorem, A
(
B
+
D
)
>
C
+
E and Routh-Hurwitz criterion, we know that all roots of (16) have negative real part; that is to say, the equilibrium E
* is locally stable for τ
=
0, while τ
0 is the minimum τ
n at which the real parts of these roots are zero. So, E
* is locally asymptotically stable if τ
∈
[
0
,
τ
0
) and unstable if τ
>
τ
0.


## 4. Direction and Stability of the Hopf Bifurcation

In the previous section, we obtain the conditions that a family periodic solutions bifurcate from the steady stateE
* at the critical value τ
n. Throughout this section, we assume that these conditions hold. As pointed by Hassard et al. [12], it is interesting to determine the direction, stability, and period of these periodic solutions bifurcating from the steady state. In this section, we will follow the idea of Ross [1] to derive the explicit formulas determining these factors. Let u
1
=
S
-
S
*, u
2
=
I
-
I
*, and u
3
=
V
-
V
*. Equation (3) becomes
(34)
S
˙
(
t
)
=
b
1
-
λ
1
(
u
1
+
S
*
)
(
u
2
+
I
*
)
-
λ
2
(
u
1
+
S
*
)
(
u
3
+
V
*
)
-
μ
1
(
u
1
+
S
*
)
,
I
˙
(
t
)
=
λ
1
(
u
1
+
S
*
)
(
u
2
+
I
*
)
+
λ
2
(
u
1
+
S
*
)
(
u
3
+
V
*
)
-
γ
(
u
2
+
I
*
)
-
μ
1
(
u
2
+
I
*
)
,
V
˙
(
t
)
=
λ
3
(
b
2
μ
2
-
u
3
-
V
*
)
(
u
2
(
t
-
τ
)
+
I
*
)
-
μ
2
(
u
3
+
V
*
)
.The linearization of (34) at u
=
(
u
1
,
u
2
,
u
3
)
=
(
0,0
,
0
) is
(35)
u
˙
1
=
-
l
1
u
1
-
l
2
u
2
-
l
3
u
3
,
u
˙
2
=
l
1
u
1
+
m
2
u
2
+
m
3
u
3
,
u
˙
3
=
-
n
2
u
2
(
t
-
τ
)
-
n
3
u
3
,
where l
1
=
λ
1
I
*
+
λ
2
V
*
+
μ
1, l
2
=
λ
1
S
*, m
2
=
λ
1
S
*
-
γ
-
μ
1, m
3
=
λ
2
S
*, n
2
=
λ
3
(
b
2
/
μ
2
-
V
*
), and n
3
=
λ
3
I
*
-
μ
2. Let x
i
(
t
)
=
u
i
(
t
τ
) and τ
=
τ
n
+
μ
(
n
=
0,1
,
2
,
…
) and τ
n is defined in (22) and μ
∈
R, and the system (34) can be written as FDE in C
=
C
(
[
-
1,0
]
,
R
3
) as
(36)
x
˙
(
t
)
=
L
μ
(
x
t
)
+
F
(
μ
,
x
t
)
,
where x
t
(
θ
)
=
x
(
t
+
θ
)
∈
C,L
μ
:
C
→
R,F
:
R
×
C
→
R are given, respectively, by
(37)
L
μ
(
ϕ
)
=
(
τ
n
+
μ
)
(
-
l
1
-
l
2
-
l
3
l
1
m
2
m
3
0
0
-
n
3
)
(
ϕ
1
(
0
)
ϕ
2
(
0
)
ϕ
3
(
0
)
)
+
(
τ
n
+
μ
)
(
0
0
0
0
0
0
0
-
n
2
0
)
(
ϕ
1
(
-
1
)
ϕ
2
(
-
1
)
ϕ
3
(
-
1
)
)
,
F
(
μ
,
x
t
)
=
(
τ
n
+
μ
)
(
-
λ
1
ϕ
1
(
0
)
ϕ
2
(
0
)
-
λ
2
ϕ
1
(
0
)
ϕ
3
(
0
)
λ
1
ϕ
1
(
0
)
ϕ
2
(
0
)
+
λ
2
ϕ
1
(
0
)
ϕ
3
(
0
)
-
λ
3
ϕ
2
(
-
1
)
ϕ
3
(
0
)
)
,
where ϕ
(
θ
)
=
(
ϕ
1
(
θ
)
,
ϕ
2
(
θ
)
,
ϕ
3
(
θ
)
)
T
∈
C. By the Riezs representation theorem, there exists a function η
(
θ
,
μ
) of bounded variation for θ
∈
[
-
1,0
], such that
(38)
L
μ
ϕ
=
∫
-
1
0
d
η
(
θ
,
0
)
ϕ
(
θ
)
,
ϕ
∈
C
.
In fact, we can choose
(39)
η
(
θ
,
μ
)
=
(
τ
n
+
μ
)
(
-
l
1
-
l
2
-
l
3
l
1
m
2
m
3
0
0
-
n
3
)
δ
(
θ
)
-
(
τ
n
+
μ
)
(
0
0
0
0
0
0
0
-
n
2
0
)
δ
(
θ
+
1
)
,
where δ is Dirac delta function.Forϕ
∈
C
1
(
[
-
1,0
]
,
R
3
), define
(40)
A
(
μ
)
ϕ
=
{
d
ϕ
(
θ
)
d
θ
θ
∈
[
-
1,0
)
,
∫
-
1
0
d
η
(
s
,
μ
)
ϕ
(
s
)
θ
=
0
,
R
(
μ
)
ϕ
=
{
0
,
θ
∈
[
-
1,0
)
,
F
(
μ
,
ϕ
)
,
θ
=
0
.
Then system (36) is equivalent to
(41)
x
˙
t
=
A
(
μ
)
x
t
+
R
(
μ
)
x
t
,
where
(42)
x
t
(
θ
)
=
x
(
t
+
θ
)
θ
∈
[
-
1,0
)
.
For φ
∈
C
1
(
[
0,1
]
,
(
R
3
)
*
), define
(43)
A
*
ψ
(
s
)
=
{
-
d
ψ
(
s
)
d
s
,
s
∈
[
0,1
)
,
∫
-
1
0
‍
d
η
T
(
t
,
0
)
ψ
(
-
t
)
,
s
=
0
,
and bilinear inner product
(44)
〈
ψ
(
s
)
,
ϕ
(
s
)
〉
=
ψ
T
¯
(
0
)
ϕ
(
0
)
-
∫
-
1
0
‍
∫
ξ
=
0
θ
‍
ψ
T
¯
(
ξ
-
θ
)
d
η
(
θ
)
ϕ
(
ξ
)
d
ξ
,
where η
(
θ
)
=
η
(
θ
,
0
). Then A
(
0
) and A
* are adjoint operators. From Section 3, we know that ±
i
w
0 are eigenvalues of A
(
0
). Thus they are eigenvalues of A
*. We need to compute the eigenvector of A
(
0
) and A
* corresponding to i
w
0 and -
i
w
0, respectively.Supposeq
(
θ
)
=
(
1
,
α
,
β
)
T
e
i
θ
w
0
τ
n is the eigenvector of A
(
0
) corresponding to i
w
0. Then A
(
0
)
q
(
θ
)
=
i
w
0
τ
n
q
(
θ
). It follows from the definition of A
(
0
) and η
(
θ
,
μ
) that
(45)
τ
n
(
-
l
1
-
i
w
0
-
l
2
-
l
3
l
1
m
2
-
i
w
0
m
3
0
-
n
2
e
-
i
w
0
τ
n
-
n
3
-
i
w
0
)
q
(
0
)
=
(
0
0
0
)
.
Then we can get
(46)
q
(
0
)
=
(
1
,
α
,
β
)
=
(
1
,
l
1
n
3
l
3
n
2
e
-
i
w
0
τ
n
-
l
2
n
3
,
-
l
1
n
2
l
3
n
2
e
-
i
w
0
τ
n
-
l
2
n
3
e
-
i
w
0
τ
n
)
T
.
We can suppose that q
*
(
s
)
=
D
(
1
,
α
*
,
β
*
)
e
i
s
w
0
τ
n is the eigenvector of A
* corresponding to -
i
w
0
τ
n, and similarly we can obtain
(47)
τ
n
(
-
l
1
+
i
w
0
l
1
0
-
l
2
m
2
+
i
w
0
-
n
2
e
i
w
0
τ
n
-
l
3
m
3
-
n
3
+
i
w
0
)
q
*
(
0
)
=
(
0
0
0
)
,
α
*
=
l
1
-
i
w
0
l
1
,
β
*
=
l
1
m
2
+
w
0
2
-
l
1
l
2
-
i
w
0
(
m
2
-
l
1
)
l
1
n
2
e
i
w
0
τ
n
.
By (44), we get
(48)
〈
q
*
(
s
)
,
q
(
θ
)
〉
=
D
¯
(
1
,
α
¯
*
,
β
¯
*
)
(
1
,
α
,
β
)
T
-
∫
-
1
0
‍
∫
ξ
=
0
θ
‍
D
¯
(
1
,
α
¯
*
,
β
¯
*
)
e
-
i
w
0
(
ξ
-
θ
)
d
η
(
θ
)
(
1
,
α
,
β
)
T
e
i
w
0
ξ
d
ξ
=
D
¯
{
1
+
α
¯
*
α
+
β
¯
*
β
-
∫
-
1
0
‍
(
1
,
α
¯
*
,
β
¯
*
)
θ
e
i
w
0
θ
d
η
(
θ
)
(
1
,
α
,
β
)
T
}
=
D
¯
{
1
+
α
¯
*
α
+
β
¯
*
β
-
τ
n
n
2
α
β
¯
*
e
-
i
w
0
τ
n
}
.
Then we choose
(49)
D
¯
=
1
1
+
α
¯
*
α
+
β
¯
*
β
-
τ
n
n
2
α
β
¯
*
e
-
i
w
0
τ
n
,
such that 〈
q
*
(
s
)
,
q
(
θ
)
〉
=
1
,
〈
q
*
(
s
)
,
q
¯
(
θ
)
〉
=
0.In the following, we use the ideas by Ross [1] to compute the coordinates describing center manifold C
0 at μ
=
0. Define
(50)
z
(
t
)
=
〈
q
*
,
x
t
〉
,
W
(
t
,
θ
)
=
x
t
-
2
Re
z
(
t
)
q
(
θ
)
.
On the center manifold C
0, we have
(51)
W
(
t
,
θ
)
=
W
(
z
(
t
)
,
z
¯
(
t
)
,
θ
)
=
W
20
(
θ
)
z
2
2
+
W
11
(
θ
)
z
z
¯
+
W
02
(
θ
)
z
¯
2
2
+
⋯
,
where z and z
¯ are local coordinates for center manifold C
0 in the direction of q
* and q
*
¯. Note that W is real if x
t is real. We can only consider real solutions. For the solution x
t
∈
C
0, since μ
=
0 and (41), we have
(52)
z
˙
=
i
w
0
τ
n
z
+
〈
q
*
(
θ
)
,
F
(
0
,
W
(
z
¯
,
z
,
θ
)
+
2
Re
z
q
(
θ
)
)
〉
=
i
w
0
τ
n
z
+
q
¯
*
(
0
)
F
(
0
,
W
(
z
¯
,
z
,
0
)
+
2
Re
z
q
(
0
)
)
=
def
i
w
0
τ
n
z
+
q
¯
*
(
0
)
F
0
(
z
,
z
¯
)
=
i
w
0
τ
n
z
+
g
(
z
,
z
¯
)
,
where
(53)
g
(
z
,
z
¯
)
=
q
¯
*
(
0
)
F
0
(
z
,
z
¯
)
=
g
20
z
2
2
+
g
11
z
z
¯
+
g
02
z
¯
2
2
+
g
21
z
2
z
¯
2
+
⋯
.
From (50), we have x
t
=
(
x
1
t
(
θ
)
,
x
2
t
(
θ
)
,
x
3
t
(
θ
)
)
=
W
(
t
,
θ
)
+
z
q
(
θ
)
+
z
q
(
θ
)
¯ and q
(
θ
)
=
(
1
,
α
,
β
)
T
e
i
w
0
τ
n, and then
(54)
x
1
t
(
0
)
=
z
+
z
¯
+
W
20
(
1
)
(
0
)
z
2
2
+
W
11
(
1
)
(
0
)
z
z
¯
+
W
02
(
1
)
(
0
)
z
¯
2
2
+
O
(
|
(
z
,
z
¯
)
|
3
)
,
x
2
t
(
0
)
=
z
α
+
z
α
¯
+
W
20
(
2
)
(
0
)
z
2
2
+
W
11
(
2
)
(
0
)
z
z
¯
+
W
02
(
2
)
(
0
)
z
¯
2
2
+
O
(
|
(
z
,
z
¯
)
|
3
)
,
x
3
t
(
0
)
=
z
β
+
z
β
¯
+
W
20
(
3
)
(
0
)
z
2
2
+
W
11
(
3
)
(
0
)
z
z
¯
+
W
02
(
3
)
(
0
)
z
¯
2
2
+
O
(
|
(
z
,
z
¯
)
|
3
)
,
x
2
t
(
-
1
)
=
z
α
e
-
i
w
0
τ
n
+
z
α
¯
e
i
w
0
τ
n
+
W
20
(
2
)
(
-
1
)
z
2
2
+
W
11
(
2
)
(
-
1
)
z
z
¯
+
W
02
(
2
)
(
-
1
)
z
¯
2
2
+
O
(
|
(
z
,
z
¯
)
|
3
)
.
From the definition of F
(
μ
,
x
t
), we have
(55)
g
(
z
,
z
¯
)
=
q
¯
*
(
0
)
F
0
(
z
,
z
¯
)
=
D
¯
τ
n
(
1
,
α
¯
*
,
β
¯
*
)
(
-
λ
1
x
1
t
(
0
)
x
2
t
(
0
)

-

λ
2
x
1
t
(
0
)
x
3
t
(
0
)
λ
1
x
1
t
(
0
)
x
2
t
(
0
)

+

λ
2
x
1
t
(
0
)
x
3
t
(
0
)
-
λ
3
x
2
t
(
-
1
)
x
3
t
(
0
)
)
=
D
¯
τ
n
{
+
α
W
11
(
3
)
(
0
)
e
-
i
w
0
τ
n
β
W
11
(
2
)
(
-
1
)
+
β
¯
2
W
20
(
2
)
(
-
1
)
)
λ
1
(
α
¯
*
-
1
)
(
W
11
(
2
)
(
0
)
+
1
2
W
20
(
2
)
(
0
)
]
+
⋯
z
2
[
α
(
α
¯
*
-
1
)
(
λ
1
+
β
λ
2
)
-
λ
3
α
β
β
¯
*
e
-
i
w
0
τ
n
]
+
2
z
z
¯
[
λ
3
β
¯
*
R
e
{
β
¯
α
e
-
i
w
0
τ
n
}
λ
1
Re
{
α
}
(
α
¯
*
-
1
)
+
λ
2
Re
{
α
β
¯
}
(
α
¯
*
-
1
)
-
λ
3
β
¯
*
Re
{
β
¯
α
e
-
i
w
0
τ
n
}
]
+
z
¯
2
[
λ
1
α
¯
(
α
¯
*

-

1
)

+

λ
2
α
β
¯
(
α
¯
*

-

1
)

-

λ
3
α
β
β
*
¯
e
i
w
0
τ
n
]
+
z
2
z
¯
[
+
α
W
11
(
3
)
(
0
)
e
-
i
w
0
τ
n
β
W
11
(
2
)
(
-
1
)
+
β
¯
2
W
20
(
2
)
(
-
1
)
)
λ
1
(
α
¯
*
-
1
)
(
+
α
¯
2
W
20
(
1
)
(
0
)
+
α
W
11
(
1
)
(
0
)
)
W
11
(
2
)
(
0
)
+
1
2
W
20
(
2
)
(
0
)
+
α
¯
2
W
20
(
1
)
(
0
)
+
α
W
11
(
1
)
(
0
)
)
+
λ
2
(
α
¯
*
-
1
)
(
+
β
¯
2
W
20
(
2
)
(
0
)
+
β
W
11
(
2
)
(
0
)
α
W
11
(
3
)
(
0
)
+
α
¯
2
W
20
(
3
)
(
0
)
+
β
¯
2
W
20
(
2
)
(
0
)
+
β
W
11
(
2
)
(
0
)
)
-
λ
3
β
¯
*
(
β
W
11
(
2
)
(
-
1
)
+
β
¯
2
W
20
(
2
)
(
-
1
)
+
α
¯
2
W
20
(
3
)
(
0
)
e
i
w
0
τ
n
+
α
W
11
(
3
)
(
0
)
e
-
i
w
0
τ
n
β
W
11
(
2
)
(
-
1
)
+
β
¯
2
W
20
(
2
)
(
-
1
)
)
λ
1
(
α
¯
*
-
1
)
(
W
11
(
2
)
(
0
)
+
1
2
W
20
(
2
)
(
0
)
]
+
⋯
}
.
Comparing the coefficients with (53), we obtain
(56)
g
20
=
2
D
¯
τ
n
[
α
(
α
¯
*
-
1
)
(
λ
1
+
β
λ
2
)
-
λ
3
α
β
β
¯
*
e
-
i
w
0
τ
n
]
,
g
11
=
2
D
¯
τ
n
[
λ
3
β
¯
*
R
e
{
β
¯
α
e
-
i
w
0
τ
n
}
λ
1
Re
{
α
}
(
α
¯
*
-
1
)
+
λ
2
Re
{
α
β
¯
}
(
α
¯
*
-
1
)
-
λ
3
β
¯
*
Re
{
β
¯
α
e
-
i
w
0
τ
n
}
]
,
g
02
=
2
D
¯
τ
n
[
λ
1
α
¯
(
α
¯
*

-

1
)

+

λ
2
α
β
¯
(
α
¯
*

-

1
)

-

λ
3
α
β
β
*
¯
e
i
w
0
τ
n
]
,
g
21
=
2
D
¯
τ
n
[
∫
λ
1
(
α
¯
*
-
1
)
(
W
11
(
2
)
(
0
)

+

1
2
W
20
(
2
)
(
0
)
+
α
¯
2
W
20
(
1
)
(
0
)
+
α
W
11
(
1
)
(
0
)
W
11
(
2
)
(
0
)

+

1
2
W
20
(
2
)
(
0
)
+
α
¯
2
W
20
(
1
)
(
0
)
)
+
λ
2
(
α
¯
*
-
1
)
(
β
¯
2
W
20
(
2
)
(
0
)
+
β
W
11
(
2
)
(
0
)
α
W
11
(
3
)
(
0
)
+
α
¯
2
W
20
(
3
)
(
0
)
+
β
¯
2
W
20
(
2
)
(
0
)
+
β
W
11
(
2
)
(
0
)
)
-
λ
3
β
¯
*
(
β
W
11
(
2
)
(
-
1
)
+
β
¯
2
W
20
(
2
)
(
-
1
)
+
α
¯
2
W
20
(
3
)
(
0
)
e
i
w
0
τ
n
+
α
W
11
(
3
)
(
0
)
e
-
i
w
0
τ
n
β
W
11
(
2
)
(
-
1
)
+
β
¯
2
W
20
(
2
)
(
-
1
)
)
]
.
In order to determine g
21, we need to compute W
20
(
θ
) and W
11
(
θ
). From (41) and (50), we have
(57)
W
˙
=
x
˙
t
-
z
˙
q
-
z
¯
˙
q
¯
=
{
A
W
-
2
Re
{
q
*
¯
(
0
)
F
0
q
(
θ
)
}
,
θ
∈
[
-
1,0
)
,
A
W
-
2
Re
{
q
*
¯
(
0
)
F
0
q
(
θ
)
}
+
F
0
,
θ
=
0
,
=
def
A
W
+
H
(
z
,
z
¯
,
θ
)
,
where
(58)
H
(
z
,
z
¯
,
θ
)
=
H
20
(
θ
)
z
2
2
+
H
11
(
θ
)
z
z
¯
+
H
02
z
¯
2
2
+
⋯
.
Note that, on the center manifold C
0 near to the origin,
(59)
W
˙
=
W
z
z
˙
+
W
z
¯
z
¯
˙
.
Thus we obtain
(60)
(
A
-
2
i
w
0
τ
n
)
W
20
(
θ
)

=

-
H
20
(
θ
)
,
A
W
11
(
θ
)

=

-
H
11
(
θ
)
.
By (57), we know that, for θ
∈
[
-
1,0
),
(61)
H
(
z
,
z
¯
,
θ
)
=
-
q
¯
*
(
0
)
F
0
q
(
θ
)
-
q
*
(
0
)
F
¯
0
q
¯
(
θ
)
=
-
g
q
(
θ
)
-
g
q
¯
(
θ
)
.
Comparing the coefficients with (58), we can get
(62)
H
20
(
θ
)
=
-
g
20
q
(
θ
)
-
g
¯
02
q
¯
(
θ
)
,
H
11
(
θ
)
=
-
g
11
q
(
θ
)
-
g
¯
11
q
¯
(
θ
)
.
From (60), (62), and the definition of A, we have
(63)
W
˙
20
(
θ
)
=
2
i
w
0
τ
n
W
20
(
θ
)
+
g
20
q
(
θ
)
+
g
¯
02
q
¯
(
θ
)
.
Noticing q
(
θ
)
=
q
(
0
)
e
i
θ
w
0
τ
n, we have
(64)
W
20
(
θ
)
=
i
g
20
τ
n
w
0
q
(
0
)
e
i
τ
n
w
0
θ

+

i
g
¯
02
3
τ
n
w
0
q
¯
(
0
)
e
-
i
τ
n
w
0
θ

+

E
1
e
2
i
τ
n
w
0
θ
,
where E
1
=
(
E
1
(
1
)
,
E
1
(
2
)
,
E
1
(
3
)
)
∈
R
3 is a constant vector. Similarly, we can have
(65)
W
11
(
θ
)
=
-
i
g
11
τ
n
w
0
q
(
0
)
e
i
τ
n
w
0
θ
+
i
g
¯
11
τ
n
w
0
q
¯
(
0
)
e
-
i
τ
n
w
0
θ
+
E
2
,
where E
2
=
(
E
2
(
1
)
,
E
2
(
2
)
,
E
2
(
3
)
)
∈
R
3 is a constant vector.In the following, we wiiill find outE
1 and E
2. From the definition of A and (60), we can obtain
(66)
∫
-
1
0
d
η
(
θ
)
W
20
(
θ
)
=
2
i
w
0
τ
n
W
20
(
0
)
-
H
20
(
0
)
,
(67)
∫
-
1
0
d
η
(
θ
)
W
11
(
θ
)
=
-
H
11
(
0
)
,
where η
(
θ
)
=
η
(
0
,
θ
). From (57) and (58), we have
(68)
H
20
(
0
)
=
-
g
20
q
(
0
)
-
g
¯
02
q
¯
(
0
)
+
2
τ
n
(
-
λ
1
α
-
λ
2
α
β
λ
1
α
+
λ
2
α
β
-
λ
3
α
β
e
-
i
w
0
τ
n
)
,
(69)
H
11
(
0
)
=
-
g
11
q
(
0
)
-
g
¯
11
q
¯
(
0
)
+
2
τ
n
(
-
2
λ
1
Re
(
α
)
-
2
λ
2
Re
(
α
β
¯
)
2
λ
1
Re
(
α
)
+
2
λ
2
Re
(
α
β
¯
)
-
2
λ
3
Re
(
α
β
¯
e
-
i
w
0
τ
n
)
)
.
Substituting (68) into (66) and noticing that
(70)
(
i
w
0
τ
n
I
-
∫
-
1
0
e
i
θ
w
0
τ
n
d
η
(
θ
)
)
q
(
0
)
=
0
,
(
-
i
w
0
τ
n
I
-
∫
-
1
0
e
-
i
θ
w
0
τ
n
d
η
(
θ
)
)
q
(
0
)
=
0
,
we then obtain
(71)
(
2
i
τ
w
0
I
-
∫
-
1
0
e
2
i
θ
τ
w
0
d
η
(
θ
)
)
E
1
=
2
τ
n
(
-
λ
1
α
-
λ
2
α
β
λ
1
α
+
λ
2
α
β
-
λ
3
α
β
e
-
i
w
0
τ
n
)
,
which is
(72)
(
2
i
w
0
+
l
1
l
2
l
3
-
l
1
2
i
w
0
-
m
2
-
m
3
0
n
2
e
-
2
i
w
0
τ
n
2
i
w
0
+
n
3
)

E
1

=

2

(
-
λ
1
α
-
λ
2
α
β
λ
1
α
+
λ
2
α
β
-
λ
3
α
β
e
-
i
w
0
τ
n
)

.
Solving the equation, we get
(73)
E
1
(
1
)
=
2
R
|
-
λ
1
α
-
λ
2
α
β
l
2
l
3
λ
1
α
+
λ
2
α
β
2
i
w
0
-
m
2
-
m
3
-
λ
3
α
β
e
-
i
w
0
τ
n
n
2
e
-
2
i
w
0
τ
n
2
i
w
0
+
n
3
|
,
E
1
(
2
)
=
2
R
|
2
i
w
0
+
l
1
-
λ
1
α
-
λ
2
α
β
l
3
-
l
1
λ
1
α
+
λ
2
α
β
-
m
3
0
-
λ
3
α
β
e
-
i
w
0
τ
n
2
i
w
0
+
n
3
|
,
E
1
(
3
)
=
2
R
|
2
i
w
0
+
l
1
l
2
-
λ
1
α
-
λ
2
α
β
-
l
1
2
i
w
0
-
m
2
λ
1
α
+
λ
2
α
β
0
n
2
e
-
2
i
w
0
τ
n
-
λ
3
α
β
e
-
i
w
0
τ
n
|
,
where
(74)
R
=
|
-
λ
1
α
-
λ
2
α
β
l
2
l
3
λ
1
α
+
λ
2
α
β
2
i
w
0
-
m
2
-
m
3
-
λ
3
α
β
e
-
i
w
0
τ
n
n
2
e
-
2
i
w
0
τ
n
2
i
w
0
+
n
3
|
.
Similarly, we can get E
2:
(75)
(
-
l
1
-
l
2
-
l
3
l
1
m
2
m
3
0
-
n
2
-
n
3
)
E
2
=
4
(
λ
1
Re
(
α
)
+
λ
2
Re
(
α
β
¯
)
-
λ
1
Re
(
α
)
-
λ
2
Re
(
α
β
¯
)
λ
3
Re
(
α
¯
β
e
i
w
0
τ
n
)
)
,
E
2
(
1
)
=
4
S
|
λ
1
Re
(
α
)
+
λ
2
Re
(
α
β
¯
)
-
l
2
-
l
3
-
λ
1
Re
(
α
)
-
λ
2
Re
(
α
β
¯
)
m
2
m
3
λ
3
Re
(
α
¯
β
e
i
w
0
τ
n
)
-
n
2
-
n
3
|
,
E
2
(
2
)
=
4
S
|
-
l
1
λ
1
Re
(
α
)
+
λ
2
Re
(
α
β
¯
)
-
l
3
l
1
-
λ
1
Re
(
α
)
-
λ
2
Re
(
α
β
¯
)
m
3
0
λ
3
Re
(
α
¯
β
e
i
w
0
τ
n
)
-
n
3
|
,
E
2
(
3
)
=
4
S
|
-
l
1
-
l
2
λ
1
Re
(
α
)
+
λ
2
Re
(
α
β
¯
)
l
1
m
2
-
λ
1
Re
(
α
)
-
λ
2
Re
(
α
β
¯
)
0
-
n
2
λ
3
Re
(
α
¯
β
e
i
w
0
τ
n
)
|
,
where
(76)
S
=
|
-
l
1
-
l
2
-
l
3
l
1
m
2
m
3
0
-
n
2
-
n
3
|
.
Therefor, all g
i
j have been expressed in terms of parameters. And we can compute the following values:
(77)
c
1
(
0
)
=
i
2
τ
n
w
0
(
g
11
g
20
-
2
|
g
11
|
2
-
|
g
02
|
2
3
)
+
g
21
2
,
μ
2
=
-
Re
{
c
1
(
0
)
}
Re
{
λ
′
(
τ
n
)
}
,
β
2
=
2
Re
{
c
1
(
0
)
}
,
T
2
=
-
Im
⁡
{
c
1
(
0
)
}
+
μ
2
Im
⁡
λ
′
(
τ
n
)
τ
n
w
0
,
which determine the quantities of bifurcating periodic solutions in the center manifold at the critical value τ
n; μ
2 determines the direction of the Hopf bifurcation. If μ
2
>
0
(
μ
2
<
0
), then the Hopf bifurcation is supercritical (subcritical); β
2 determines the stability of the bifurcating periodic solutions: the periodic solutions are stable (unstable) if β
2
<
0
(
β
2
>
0
); T
2 determines the period of the bifurcating solutions: the periodic increases (decreases) if T
2
>
0
(
<
0
).



---

*Source: 804204-2014-05-20.xml*