# Hopf Bifurcation and Stability of Periodic Solutions for Delay Differential Model of HIV Infection of CD4+ T-cells

Authors: P. Balasubramaniam; M. Prakash; Fathalla A. Rihan; S. Lakshmanan
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This paper deals with stability and Hopf bifurcation analyses of a mathematical model of HIV infection ofCD
4
+ T-cells. The model is based on a system of delay differential equations with logistic growth term and antiretroviral treatment with a discrete time delay, which plays a main role in changing the stability of each steady state. By fixing the time delay as a bifurcation parameter, we get a limit cycle bifurcation about the infected steady state. We study the effect of the time delay on the stability of the endemically infected equilibrium. We derive explicit formulae to determine the stability and direction of the limit cycles by using center manifold theory and normal form method. Numerical simulations are presented to illustrate the results.


---

## Body

## 1. Introduction

Since 1980, the human immunodeficiency virus (HIV) or the associated syndrome of opportunistic infections that causes acquired immunodeficiency syndrome (AIDS) has been considered as one of the most serious global public health menaces. When HIV enters the body, its main target is the CD4 lymphocytes, also called CD4 T-cells (includingCD
4
+ T-cells). When a CD4 cell is infected with HIV, the virus goes through multiple steps to reproduce itself and create many more virus particles. The AIDS term, which is known as the late stage of HIV, covers the range of infections and illnesses which can result from a weakened immune system caused by HIV. Based on the clinical studies, it is known that, for a normal person, the CD
4
+ T-cells count is around 1000
mm
-
3 and for HIV infected patient it gradually decreases to 200
mm
-
3 or below, which leads to AIDS. However, this may take several years for the number of CD4 T-cells to reduce to a level where the immune system is weakened [1–6].Mathematical models, usingdelay differential equations (DDEs), have provided insights in understanding the dynamics of HIV infection. Discrete or continuous time delays have been introduced to the models to describe the time between infection of aCD
4
+ T-cell and the emission of viral particles on a cellular level [7–13]. In general, DDEs exhibit much more complicated dynamics than ODEs since the time delay could cause a stable equilibrium to become unstable and cause the populations to fluctuate [14–16]. In studying the viral clearance rates, Perelson et al. [17] assumed that there are two types of delays that occur between the administration of drug and the observed decline in viral load: a pharmacological delay that occurs between the ingestion of drug and its appearance within cells and an intracellular delay that is between initial infection of a cell by HIV and the release of new virion. In this paper, we incorporate an intracellular delay to the model to describe the time between infection of a CD
4
+ T-cell and the emission of viral particles on a cellular level [18]. We study the impact of the presence of such time delay on the dynamics of the model.The outline of the present paper is as follows. In Section2, we describe the model. In Section 3, we study the qualitative behavior of the model via stability of the steady states and Hopf bifurcation when time delay is considered as a bifurcation parameter. In Section 4, we provide an explicit formula to determine the direction of bifurcating periodic solution by applying center manifold theory and normal form method. We provide some numerical simulations to demonstrate the effectiveness of the analysis in Section 5 and we conclude in Section 6.


## 2. Description of the Model

Let us start the analysis with some basic models of the dynamics of target (uninfected) cells and infectedCD
4
+ T-cells by HIV. As a first approximation, the dynamics between HIV and the macrophage population was described by the simplest model of infection dynamics presented in [19–21]. Denoting uninfected cells by x
(
t
) and infected cells by y
(
t
) and assuming that viruses are transmitted mainly by cell to cell contact, the model is given by
(1)
x
˙
(
t
)
=
Λ
-
δ
1
x
(
t
)
-
β
x
(
t
)
y
(
t
)
,
y
˙
(
t
)
=
β
x
(
t
)
y
(
t
)
-
δ
2
y
(
t
)
.
The target (uninfected) CD
4
+ T-cells are produced at a rate Λ, die at a rate δ
1, and become infected by virus at a rate β. The infected host cells die at a rate δ
2. The basic reproductive ratio of the virus is then given by R
0
=
Λ
β
/
δ
1
δ
2. If there is no infection or if R
0
<
1, there is only trivial equilibrium (E
0
=
(
Λ
/
δ
1
,
0
)) with no virus-producing cells. Whereas if R
0
>
1, the virus can establish an infection and the system converges to the equilibrium with both uninfected cells and infected cells, E
1
=
(
δ
2
/
β
,
Λ
/
δ
2
-
δ
1
/
β
).However, in most viral infections, the CTL response plays a crucial part in antiviral defence by attacking viral infected cells [22, 23]. As the the cytotoxic T-lymphocyte (CTL) immune response is necessary to eliminate or control the viral infection, we incorporated the antiviral CTL immune response into the basic model (1). Therefore, if we add CTL response, which is denoted by z
(
t
), into model (1) (see [19]), then the extended model is
(2)
x
˙
(
t
)
=
Λ
-
δ
1
x
(
t
)
-
β
x
(
t
)
y
(
t
)
,
y
˙
(
t
)
=
β
x
(
t
)
y
(
t
)
-
δ
2
y
(
t
)
-
p
y
(
t
)
z
(
t
)
,
z
˙
(
t
)
=
c
q
y
(
t
)
z
(
t
)
-
h
z
(
t
)
.
Thus, CTLs proliferate in response to antigen at a rate c, die at a rate h, and lyse infected cells at a rate p. We assume that the CTL pool consists of two populations: the precursors w
(
t
) and the effectors z
(
t
). In other words, we assume that there are primary and secondary responses to viral infections. Then, the model (2) becomes
(3)
x
˙
(
t
)
=
Λ
-
δ
1
x
(
t
)
-
β
x
(
t
)
y
(
t
)
,
y
˙
(
t
)
=
β
x
(
t
)
y
(
t
)
-
δ
2
y
(
t
)
-
p
y
(
t
)
z
(
t
)
,
w
˙
(
t
)
=
c
(
1
-
q
)
y
(
t
)
w
(
t
)
-
b
w
(
t
)
,
z
˙
(
t
)
=
c
q
y
(
t
)
w
(
t
)
-
h
z
(
t
)
.
The infected cells are killed by CTL effector cells at a rate p
y
z. Upon contact with antigen, CTLp proliferate at a rate c
y
(
t
)
w
(
t
) and differentiate into effector cells CTLe at a rate c
q
y
(
t
)
w
(
t
). CTL precursors die at a rate b
w, and effectors die at a rate h
z
(
t
); see Figure 1.Figure 1
A simplified model of virus-CTL interaction. The virus dynamics is described by the basic model of Nowak and Bangham [19]. The uninfected target cells are produced at a rate Λ and die at a rate δ
1
x. They become infected by the virus at a rate β
x
y. The infected cells produce new virus particle and die at a rate δ
2
y. When CT
L
p recognize antigen on the surface of infected cells, they become activated and expand at a rate c
y
w, decay at a rate b
w, and differentaite into efector cells at a rate c
q
w
y. The effector cells lyse the infected cells at a rate p
y
z.Since the proliferation ofCD
4
+ T-cells is density dependent, that is, the rate of proliferation decreases as T-cells increase and reach the carrying capacity, we then extend the above basic viral infection model to include the density dependent growth of the CD
4
+ T-cell population (see [24–26]). It is also known that HIV infection leads to low levels of CD
4
+ T-cells via three main mechanisms: direct viral killing of infected cells, increased rates of apoptosis in infected cells, and killing of infected CD
4
+ T-cells by cytotoxic T-lymphocytes [26]. Hence, it is reasonable to include apoptosis of infected cells. An average of 1
0
10 viral particles is produced by infected cells per day. The treatment with single antiviral drug is considered to be failed, so that the combination of antiviral drugs is needed for the better treatment [25]. Therefore, in the below revised model, we combine the antiretroviral drugs, namely, reverse transcriptase inhibitor (RTI) and protease inhibitor (PI) to make the model realistic (see [27–29]). RTIs can block the infection of target T-cells by infectious virus, and PIs cause infected cells to produce noninfectious virus particles. The modified model takes the form
(4)
x
˙
(
t
)
=
Λ
-
δ
1
x
(
t
)
+
r
(
1
-
x
(
t
)
+
y
(
t
)
T
max
⁡
)
x
(
t
)
-
(
1
-
ϵ
)
(
1
-
η
)
β
x
(
t
)
y
(
t
)
,
y
˙
(
t
)
=
(
1
-
ϵ
)
(
1
-
η
)
β
x
(
t
)
y
(
t
)
-
δ
2
y
(
t
)
-
e
1
y
(
t
)
-
p
y
(
t
)
z
(
t
)
,
w
˙
(
t
)
=
c
y
(
t
)
w
(
t
)
-
c
q
y
(
t
)
w
(
t
)
-
b
w
(
t
)
,
z
˙
(
t
)
=
c
q
y
(
t
)
w
(
t
)
-
h
z
(
t
)
.
The first equation of model (4) represents the rate of change in the count of healthy CD
4
+ T-cells that produced at rate Λ and become infected at rate β, with the mortality δ
1. We assume that the uninfected CD
4
+ T-cells proliferate logistically, thus the growth rate r is multiplied by the term (
1
-
(
x
+
y
)
/
T
max
⁡
) and this term approaches zero when the total number of T-cells approaches the carrying capacity T
max
⁡. The effects of combination of RTI and PI antiviral drugs are represented by the term (
1
-
ϵ
)
(
1
-
η
)
β
x
y, where (
1
-
ϵ
), 0
<
ϵ
<
1, represents the effects of RTI and (
1
-
η
), 0
<
η
<
1, represents the effects of PI. The second equation of model (4) denotes the rate of change in the count of infected CD
4
+ T-cells. The infected CD
4
+ T-cells decay at a rate δ
2 and e
1 denotes apoptosis rate of infected cell; infected cells are killed by CTL effectors at a rate p. The third equation of the model denotes the rate of change in the CTLp population; proliferation rate of the CTLp is given by c and is proportional to the infected cells y; CTLp die at a rate b and differentiate into CTL effectors at a rate c
q. The last equation of the model represents the concentration of CTL effectors, which die at a rate h. In reality, the specific immune system is not immediately effective following invasion by a novel pathogen. There may be an explicit time delay between infection and immune initiation and there may be a gradual build-up in immune efficacy during which the immune response develops, before reaching maximal specificity to the pathogen ([8, 30, 31]). In order to make model (4) more realistic, time delay in the immune response should be included in the following model:
(5)
x
˙
(
t
)
=
Λ
-
(
1
-
ϵ
)
(
1
-
η
)
β
x
(
t
)
y
(
t
)
+
r
(
1
-
x
(
t
)
+
y
(
t
)
T
max
⁡
)
x
(
t
)
-
δ
1
x
(
t
)
,
y
˙
(
t
)
=
(
1
-
ϵ
)
(
1
-
η
)
β
x
(
t
)
y
(
t
)
-
(
δ
2
+
e
1
)
y
(
t
)
-
p
y
(
t
)
z
(
t
)
,
w
˙
(
t
)
=
c
(
1
-
q
)
y
(
t
-
τ
)
w
(
t
-
τ
)
-
b
w
(
t
)
z
˙
(
t
)
=
c
q
y
(
t
-
τ
)
w
(
t
-
τ
)
-
h
z
(
t
)
.
The range of parameter values of the model are given in Table 1.Table 1
Parameter definitions and estimations used in the underlying model.
Parameter
Notes
Estimated Value
Range
Source
Λ
Source of uninfected CD4+ T-cells
10
0–10
[26]
β
Rate of infection
0.1
0.00001–0.5
[26]
T
max
⁡
Total carrying capacity
1500
1500
[26]
r
Logistic growth term
0.03
0.03–3
[26]
δ
1
Mortality rate of CD4+ T-cells
0.06
0.007–0.1
[26]
ϵ
Antiretroviral (RTI) therapy
0.9
0-1
see text
δ
2
Infected cells died out naturally
0.3
0.2–1.4
[26]
e
1
Apoptosis rate of infected cells
0.2
0.2
[26]
p
Clearance rate of infected cells
1
0.001–1
[26]
η
Protease inhibitor therapy
0.9
[0, 1]
see text
q
Rate of differentiation of CTLs
0.02
Assumed
—
b
Death rate of CTL precursors
0.02
0.005–0.15
[26]
c
Proliferation of CTLs responsiveness
0.1
0.001–1
[26]
h
Mortality rate or CTL effectors
0.1
0.005–0.15
[26]We start our analysis by presenting some notations that will be used in the sequel. LetC
=
C
(
[
-
τ
,
0
]
,
R
+
4
) be the Banach space of continuous functions mapping the interval [
-
τ
,
0
] into R
+
4, where R
+
4
=
(
x
,
y
,
w
,
z
); the initial conditions are given by
(6)
x
(
θ
)
=
φ
1
(
θ
)
≥
0
,
y
(
θ
)
=
φ
2
(
θ
)
≥
0
,
w
(
θ
)
=
φ
3
(
θ
)
≥
0
,
z
(
θ
)
=
φ
4
(
θ
)
≥
0
,
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
θ
∈
[
-
τ
,
0
]
,
where φ
i
(
θ
)
∈
C
1 are smooth functions, for all i
=
1,2
,
3,4.   From the fundamental theory of functional differential equations (see [32, 33]), it is easy to see that the solutions (
x
(
t
)
,
y
(
t
)
,
w
(
t
)
,
z
(
t
)
) of system (5) with the initial conditions as stated above exist for all t
≥
0 and are unique. It can be shown that these solutions exist for all t
>
0 and stay nonnegative. In fact, if x
(
0
)
>
0, then x
(
t
)
>
0 for all t
>
0. The same argument is true for the y,  w, and z components. Hence, the interior R
+
4 is invariant for system (5).


## 3. Steady States

We can obtain the steady state values by settingx
˙
=
y
˙
=
w
˙
=
z
˙
=
0. The steady state value of the infection-free steady sate E
0 is given by E
0
=
(
(
T
max
⁡
/
2
r
)
(
r
-
δ
1
+
(
r
-
δ
1
)
2
+
4
r
Λ
/
T
max
⁡
)
,
0,0
,
0
), while the infected steady state E
+
=
(
x
*
,
y
*
,
w
*
,
z
*
) is given by
(7)
y
*
=
b
c
(
1
-
q
)
,
w
*
=
h
(
1
-
q
)
z
*
q
b
,
z
*
=
(
1
-
ϵ
)
(
1
-
η
)
β
x
*
-
(
δ
2
+
e
1
)
p
,
and x
* is given by the following quadratic equation:
(8)
c
1
x
2
+
c
2
x
-
c
3
=
0
,
where c
1
=
c
(
1
-
q
)
r, c
2
=
T
max
⁡
b
β
(
1
-
ϵ
)
(
1
-
η
)
+
b
r
-
c
(
1
-
q
)
T
max
⁡
(
r
-
δ
1
), c
3
=
c
(
1
-
q
)
Λ
T
max
⁡.


### 3.1. Stability and Hopf Bifurcation Analysis of Infected Steady StateE
+

In order to study full dynamics of model (4) by using time delay as a bifurcation parameter, we need to linearize the model around the steady state E
+ and determine the characteristic equation of the Jacobian matrix. The roots of the characteristic equation determine the asymptotic stability and existence of Hopf bifurcation for the model. The characteristic equation of the linearized system is given by(9)
|
-
A
1
y
*
+
r
-
2
r
T
max
⁡
x
*
-
r
T
max
⁡
y
*
-
δ
1
-
λ
-
A
1
x
*
-
r
T
max
⁡
x
*
0
0
A
1
y
*
A
1
x
*
-
(
δ
2
+
e
1
)
-
p
z
*
-
λ
0
-
p
y
*
0
c
(
1
-
q
)
e
-
λ
τ
w
*
c
(
1
-
q
)
e
-
λ
τ
y
*
-
b
-
λ
0
0
c
q
e
-
λ
τ
w
*
c
q
e
-
λ
τ
y
*
-
h
-
λ
|
=
0
, which is equivalent to the equation
(10)
λ
4
+
p
1
λ
3
+
p
2
λ
2
+
p
3
λ
+
p
4
λ
4
+
e
-
λ
τ
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
=
0
,
where A
1
=
(
1
-
ϵ
)
(
1
-
η
)
β and
(11)
p
1
=
-
a
1
-
a
4
-
a
8
-
a
11
,
p
2
=
a
1
a
8
+
a
8
a
11
+
a
1
a
11
+
a
4
a
8
+
a
4
a
11
+
a
1
a
4
-
a
2
a
3
,
p
3
=
a
2
a
3
a
8
+
a
2
a
3
a
11
-
a
1
a
8
a
11
-
a
4
a
8
a
11
-
a
1
a
4
a
8
-
a
1
a
4
a
11
,
p
4
=
a
1
a
4
a
8
a
11
-
a
2
a
3
a
8
a
11
,
q
1
=
-
a
7
,
q
2
=
a
1
a
7
+
a
7
a
11
+
a
4
a
7
-
a
5
a
9
,
q
3
=
a
5
a
8
a
9
+
a
1
a
5
a
9
+
a
2
a
3
a
7
-
a
1
a
7
a
11
-
a
4
a
7
a
11
-
a
1
a
4
a
7
,
q
4
=
a
1
a
4
a
7
a
11
-
a
1
a
5
a
8
a
9
-
a
2
a
3
a
7
a
11
,
a
1
=
-
(
1
-
ϵ
)
(
1
-
η
)
β
y
*
+
r
-
2
r
x
*
T
max
⁡
-
r
y
*
T
max
⁡
-
δ
1
,
a
2
=
-
(
1
-
ϵ
)
(
1
-
η
)
β
x
*
-
r
x
*
T
max
⁡
,
a
3
=
(
1
-
ϵ
)
(
1
-
η
)
β
y
*
,
a
4
=
(
1
-
ϵ
)
(
1
-
η
)
β
x
*
-
(
δ
2
+
e
1
)
-
p
z
*
,
a
5
=
-
p
y
*
,
a
6
=
c
(
1
-
q
)
w
*
,
a
7
=
c
(
1
-
q
)
y
*
,
a
8
=
-
b
,
a
9
=
c
q
w
*
,
a
10
=
c
q
y
*
,
a
11
=
-
h
.
Let us consider the following equation:
(12)
φ
(
λ
,
τ
)
=
λ
4
+
p
1
λ
3
+
p
2
λ
2
+
p
3
λ
+
p
4
+
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
e
-
λ
τ
.
For the nondelayed model (say τ
=
0), from (10), we have
(13)
λ
4
+
D
1
λ
3
+
D
2
λ
2
+
D
3
λ
+
D
4
=
0
,
where
(14)
D
1
=
p
1
+
q
1
,
D
2
=
p
2
+
q
2
,
D
3
=
p
3
+
q
3
,
D
4
=
p
4
+
q
4
.Lemma 1.
Forτ
=
0, the unique nontrivial equilibrium is locally asymptotically stable if the real parts of all the roots of (13) are negative.Proof.
The proof of the above lemma is based on holding the following conditions:D
1
>
0,  D
3
>
0,  D
4
>
0, and D
1
D
2
D
3
>
D
1
2
D
4
+
D
3
2, as proposed by Routh-Hurwitz criterion. We conclude that equilibrium E
+ is locally asymptotically stable if and only if all the roots of the characteristic equation (13) have negative real parts which depends on the numerical values of parameters that are shown in the numerical exploration.


### 3.2.  Existence of Hopf Bifurcation

We here study the impact of the time-delay parameter on the stability of HIV infection ofCD
4
+ T-cells. We deduce criteria that ensure the asymptotic stability of infected steady state E
+, for all τ
>
0. We arrive at the following theorem.Theorem 2.
Necessary and sufficient conditions for the infected equilibriumE
+ to be asymptotically stable for all delay τ
≥
0 are as follows (i)
the real parts of all the roots ofφ
(
λ
,
τ
)
=
0 are negative;
(ii)
for allω and τ
≥
0, φ
(
i
ω
,
τ
)
≠
0, where i
=
-
1.Proof.
Assume that Lemma1 is true. Now, for ω
=
0, we have
(15)
φ
(
0
,
τ
)
=
D
4
=
p
4
+
q
4
≠
0
.
Substituting λ
=
i
ω
(
ω
>
0
) into (5) and separating the real and imaginary parts of the equations yields
(16)
(
ω
4
-
p
2
ω
2
+
p
4
)
+
(
-
q
2
ω
2
+
q
4
)
cos
⁡
(
ω
τ
)
+
(
-
q
1
ω
3
+
q
3
ω
)
sin
(
ω
τ
)
=
0
,
(
-
p
1
ω
3
+
p
3
ω
)
+
(
-
q
1
ω
3
+
q
3
ω
)
cos
⁡
(
ω
τ
)
-
(
-
q
2
ω
2
+
q
4
)
sin
(
ω
τ
)
=
0
.
After some mathematical manipulations, we obtain the following equations
(17)
cos
⁡
(
ω
τ
)
=
(
(
q
2
-
p
1
q
1
)
ω
6
+
(
p
3
q
1
-
q
4
-
p
2
q
2
+
p
1
q
3
)
ω
4
i
i
i
i
i
i
i
i
i
+
(
p
2
q
4
+
p
4
q
2
-
p
3
q
3
)
ω
2
-
p
4
q
4
)
×
(
q
1
2
ω
6
+
(
q
2
2
-
2
q
1
q
3
)
ω
4
+
(
q
3
2
-
2
q
2
q
4
)
ω
2
+
q
4
2
)
-
1
,
sin
(
ω
τ
)
=
(
2009
q
1
ω
7
+
(
p
1
q
2
-
q
3
-
p
2
q
1
)
ω
5
iiiiiiiii
+
(
p
2
q
3
+
p
4
q
1
-
p
3
q
2
-
p
1
q
4
)
ω
3
iiiiiiii
+
(
p
3
q
4
-
p
4
q
3
)
ω
2009
)
×
(
2009
q
1
2
ω
6
+
(
q
2
2
-
2
q
1
q
3
)
ω
4
+
(
q
3
2
-
2
q
2
q
4
)
ω
2
+
q
4
2
2009
)
-
1
.
Let
(18)
b
1
=
q
2
-
p
1
q
1
,
b
2
=
p
3
q
1
-
q
4
-
p
2
q
2
+
p
1
q
3
,
b
3
=
p
2
q
4
+
p
4
q
2
-
p
3
q
3
,
b
4
=
-
p
4
q
4
,
b
5
=
q
1
2
,
b
6
=
q
2
2
-
2
q
1
q
3
,
b
7
=
q
3
2
-
2
q
2
q
4
,
b
8
=
q
4
2
,
b
9
=
q
1
,
b
10
=
p
1
q
2
-
q
3
-
p
2
q
1
,
b
11
=
p
2
q
3
+
p
4
q
1
-
p
3
q
2
-
p
1
q
4
,
b
12
=
p
3
q
4
-
p
4
q
3
.
From (16), we have
(19)
ω
8
+
c
1
ω
6
+
c
2
ω
4
+
c
3
ω
2
+
c
4
=
0
,
where
(20)
c
1
=
p
1
2
-
2
p
2
-
q
1
2
,
c
2
=
p
2
2
-
2
p
1
p
3
+
2
q
1
q
3
+
2
p
4

-

q
2
2
,
c
3
=
p
3
2
-
2
p
2
p
4
+
2
q
2
q
4
-
q
3
2
,
c
4
=
p
4
2
-
q
4
2
.
The conditions (
i
) and (
ii
) of Theorem 2 hold if and only if (19) has no real positive root.Letm
=
ω
2; then (19) takes the form
(21)
m
4
+
c
1
m
3
+
c
2
m
2
+
c
3
m
+
c
4
=
0
.
If c
4
<
0, then (19) has at least one positive root. In the case when (19) has four positive roots, we have
(22)
ω
1
=
m
1
,
ω
2
=
m
2
,
ω
3
=
m
3
,
ω
4
=
m
4
.
From (16), we have
(23)
τ
k
(
j
)
=
1
ω
k
{
arcsin
b
9
ω
k
7
+
b
10
ω
k
5
+
b
11
ω
k
3
+
b
12
ω
k
b
5
ω
k
6
+
b
6
ω
k
4
+
b
7
ω
k
2
+
b
8
+
2
j
π
}
,
where k
=
1,2
,
3,4 and j
=
0,1
,
2
,
…; we choose τ
0
=
min
⁡
(
τ
k
(
j
)
).To establish Hopf bifurcation atτ
=
τ
0, we need to show that
(24)
R
(
d
λ
d
τ
)
τ
=
τ
0
≠
0
.
By differentiating (10) with respect to τ, we can get
(25)
d
λ
d
τ
=
λ
e
-
λ
τ
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
×
(
e
-
λ
τ
(
4
λ
3
+
3
p
1
λ
2
+
2
p
2
λ
+
p
3
)
+
e
-
λ
τ
i
i
i
i
i
i
i
i
i
×
[
(
3
q
1
λ
2
+
2
q
2
λ
+
q
3
)
iiiiiiiiiiiiii
-
τ
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
]
e
-
λ
τ
)
-
1
.
It follows that
(26)
(
d
λ
d
τ
)
-
1
=
(
(
4
λ
3
+
3
p
1
λ
2
+
2
p
2
λ
+
p
3
)
+
e
-
λ
τ
i
i
i
i
i
i
×
[
(
3
q
1
λ
2
+
2
q
2
λ
+
q
3
)
i
i
i
i
i
i
i
i
i
i
-
τ
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
]
)
×
(
λ
e
-
λ
τ
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
)
-
1
.
Then, by combining (10), we get
(27)
(
d
λ
d
τ
)
-
1
=
(
(
4
λ
3
+
3
p
1
λ
2
+
2
p
2
λ
+
p
3
)
i
i
i
i
i
+
e
-
λ
τ
(
3
q
1
λ
2
+
2
q
2
λ
+
q
3
)
)
×
(
λ
e
-
λ
τ
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
)
-
1
-
τ
λ
.
Substituting λ
=
i
ω
0 in (27) (where ω
0
>
0 and i
=
-
1) yields
(28)
(
d
λ
d
τ
)
-
1
|
τ
=
τ
0
=
d
1
+
i
d
2
d
3
+
i
d
4
-
τ
λ
,
where
(29)
d
1
=
(
p
3
-
3
p
1
ω
0
2
)
+
(
q
3
-
3
q
1
ω
0
2
)
cos
⁡
(
ω
0
τ
0
)
+
2
q
2
ω
0
sin
(
ω
0
τ
0
)
,
d
2
=
(
2
p
2
ω
0
-
4
ω
3
)
+
2
q
2
ω
0
cos
⁡
(
ω
0
τ
0
)
-
(
q
3
-
3
q
1
ω
0
2
)
sin
(
ω
0
τ
0
)
,
d
3

=

(
q
1
ω
0
4
-
q
3
ω
0
2
)
cos
⁡
(
ω
0
τ
0
)
+
(
q
4
ω
0
-
q
2
ω
0
3
)
sin
(
ω
0
τ
0
)
,
d
4

=

(
q
4
ω
0
-
q
2
ω
0
3
)
cos
⁡
(
ω
0
τ
0
)
-
(
q
1
ω
0
4
-
q
3
ω
0
2
)
sin
(
ω
0
τ
0
)
.
Thus,
(30)
R
(
d
λ
d
τ
)
-
1
|
τ
=
τ
0
=
d
1
d
3
+
d
2
d
4
d
3
2
+
d
4
2
.
Notice that
(31)
sign
⁡
(
R
d
λ
(
t
)
d
τ
)
|
τ
=
τ
0
=
sign
⁡
(
R
(
d
λ
d
τ
)
-
1
)
|
τ
=
τ
0
.
By summarizing the above analysis, we arrive at the following theorem.Theorem 3.
The infected equilibriumE
+ of the system (5) is asymptotically stable for τ
∈
[
0
,
τ
0
) and it undergoes Hopf bifurcation at τ
=
τ
0.


## 3.1. Stability and Hopf Bifurcation Analysis of Infected Steady StateE
+

In order to study full dynamics of model (4) by using time delay as a bifurcation parameter, we need to linearize the model around the steady state E
+ and determine the characteristic equation of the Jacobian matrix. The roots of the characteristic equation determine the asymptotic stability and existence of Hopf bifurcation for the model. The characteristic equation of the linearized system is given by(9)
|
-
A
1
y
*
+
r
-
2
r
T
max
⁡
x
*
-
r
T
max
⁡
y
*
-
δ
1
-
λ
-
A
1
x
*
-
r
T
max
⁡
x
*
0
0
A
1
y
*
A
1
x
*
-
(
δ
2
+
e
1
)
-
p
z
*
-
λ
0
-
p
y
*
0
c
(
1
-
q
)
e
-
λ
τ
w
*
c
(
1
-
q
)
e
-
λ
τ
y
*
-
b
-
λ
0
0
c
q
e
-
λ
τ
w
*
c
q
e
-
λ
τ
y
*
-
h
-
λ
|
=
0
, which is equivalent to the equation
(10)
λ
4
+
p
1
λ
3
+
p
2
λ
2
+
p
3
λ
+
p
4
λ
4
+
e
-
λ
τ
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
=
0
,
where A
1
=
(
1
-
ϵ
)
(
1
-
η
)
β and
(11)
p
1
=
-
a
1
-
a
4
-
a
8
-
a
11
,
p
2
=
a
1
a
8
+
a
8
a
11
+
a
1
a
11
+
a
4
a
8
+
a
4
a
11
+
a
1
a
4
-
a
2
a
3
,
p
3
=
a
2
a
3
a
8
+
a
2
a
3
a
11
-
a
1
a
8
a
11
-
a
4
a
8
a
11
-
a
1
a
4
a
8
-
a
1
a
4
a
11
,
p
4
=
a
1
a
4
a
8
a
11
-
a
2
a
3
a
8
a
11
,
q
1
=
-
a
7
,
q
2
=
a
1
a
7
+
a
7
a
11
+
a
4
a
7
-
a
5
a
9
,
q
3
=
a
5
a
8
a
9
+
a
1
a
5
a
9
+
a
2
a
3
a
7
-
a
1
a
7
a
11
-
a
4
a
7
a
11
-
a
1
a
4
a
7
,
q
4
=
a
1
a
4
a
7
a
11
-
a
1
a
5
a
8
a
9
-
a
2
a
3
a
7
a
11
,
a
1
=
-
(
1
-
ϵ
)
(
1
-
η
)
β
y
*
+
r
-
2
r
x
*
T
max
⁡
-
r
y
*
T
max
⁡
-
δ
1
,
a
2
=
-
(
1
-
ϵ
)
(
1
-
η
)
β
x
*
-
r
x
*
T
max
⁡
,
a
3
=
(
1
-
ϵ
)
(
1
-
η
)
β
y
*
,
a
4
=
(
1
-
ϵ
)
(
1
-
η
)
β
x
*
-
(
δ
2
+
e
1
)
-
p
z
*
,
a
5
=
-
p
y
*
,
a
6
=
c
(
1
-
q
)
w
*
,
a
7
=
c
(
1
-
q
)
y
*
,
a
8
=
-
b
,
a
9
=
c
q
w
*
,
a
10
=
c
q
y
*
,
a
11
=
-
h
.
Let us consider the following equation:
(12)
φ
(
λ
,
τ
)
=
λ
4
+
p
1
λ
3
+
p
2
λ
2
+
p
3
λ
+
p
4
+
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
e
-
λ
τ
.
For the nondelayed model (say τ
=
0), from (10), we have
(13)
λ
4
+
D
1
λ
3
+
D
2
λ
2
+
D
3
λ
+
D
4
=
0
,
where
(14)
D
1
=
p
1
+
q
1
,
D
2
=
p
2
+
q
2
,
D
3
=
p
3
+
q
3
,
D
4
=
p
4
+
q
4
.Lemma 1.
Forτ
=
0, the unique nontrivial equilibrium is locally asymptotically stable if the real parts of all the roots of (13) are negative.Proof.
The proof of the above lemma is based on holding the following conditions:D
1
>
0,  D
3
>
0,  D
4
>
0, and D
1
D
2
D
3
>
D
1
2
D
4
+
D
3
2, as proposed by Routh-Hurwitz criterion. We conclude that equilibrium E
+ is locally asymptotically stable if and only if all the roots of the characteristic equation (13) have negative real parts which depends on the numerical values of parameters that are shown in the numerical exploration.


## 3.2.  Existence of Hopf Bifurcation

We here study the impact of the time-delay parameter on the stability of HIV infection ofCD
4
+ T-cells. We deduce criteria that ensure the asymptotic stability of infected steady state E
+, for all τ
>
0. We arrive at the following theorem.Theorem 2.
Necessary and sufficient conditions for the infected equilibriumE
+ to be asymptotically stable for all delay τ
≥
0 are as follows (i)
the real parts of all the roots ofφ
(
λ
,
τ
)
=
0 are negative;
(ii)
for allω and τ
≥
0, φ
(
i
ω
,
τ
)
≠
0, where i
=
-
1.Proof.
Assume that Lemma1 is true. Now, for ω
=
0, we have
(15)
φ
(
0
,
τ
)
=
D
4
=
p
4
+
q
4
≠
0
.
Substituting λ
=
i
ω
(
ω
>
0
) into (5) and separating the real and imaginary parts of the equations yields
(16)
(
ω
4
-
p
2
ω
2
+
p
4
)
+
(
-
q
2
ω
2
+
q
4
)
cos
⁡
(
ω
τ
)
+
(
-
q
1
ω
3
+
q
3
ω
)
sin
(
ω
τ
)
=
0
,
(
-
p
1
ω
3
+
p
3
ω
)
+
(
-
q
1
ω
3
+
q
3
ω
)
cos
⁡
(
ω
τ
)
-
(
-
q
2
ω
2
+
q
4
)
sin
(
ω
τ
)
=
0
.
After some mathematical manipulations, we obtain the following equations
(17)
cos
⁡
(
ω
τ
)
=
(
(
q
2
-
p
1
q
1
)
ω
6
+
(
p
3
q
1
-
q
4
-
p
2
q
2
+
p
1
q
3
)
ω
4
i
i
i
i
i
i
i
i
i
+
(
p
2
q
4
+
p
4
q
2
-
p
3
q
3
)
ω
2
-
p
4
q
4
)
×
(
q
1
2
ω
6
+
(
q
2
2
-
2
q
1
q
3
)
ω
4
+
(
q
3
2
-
2
q
2
q
4
)
ω
2
+
q
4
2
)
-
1
,
sin
(
ω
τ
)
=
(
2009
q
1
ω
7
+
(
p
1
q
2
-
q
3
-
p
2
q
1
)
ω
5
iiiiiiiii
+
(
p
2
q
3
+
p
4
q
1
-
p
3
q
2
-
p
1
q
4
)
ω
3
iiiiiiii
+
(
p
3
q
4
-
p
4
q
3
)
ω
2009
)
×
(
2009
q
1
2
ω
6
+
(
q
2
2
-
2
q
1
q
3
)
ω
4
+
(
q
3
2
-
2
q
2
q
4
)
ω
2
+
q
4
2
2009
)
-
1
.
Let
(18)
b
1
=
q
2
-
p
1
q
1
,
b
2
=
p
3
q
1
-
q
4
-
p
2
q
2
+
p
1
q
3
,
b
3
=
p
2
q
4
+
p
4
q
2
-
p
3
q
3
,
b
4
=
-
p
4
q
4
,
b
5
=
q
1
2
,
b
6
=
q
2
2
-
2
q
1
q
3
,
b
7
=
q
3
2
-
2
q
2
q
4
,
b
8
=
q
4
2
,
b
9
=
q
1
,
b
10
=
p
1
q
2
-
q
3
-
p
2
q
1
,
b
11
=
p
2
q
3
+
p
4
q
1
-
p
3
q
2
-
p
1
q
4
,
b
12
=
p
3
q
4
-
p
4
q
3
.
From (16), we have
(19)
ω
8
+
c
1
ω
6
+
c
2
ω
4
+
c
3
ω
2
+
c
4
=
0
,
where
(20)
c
1
=
p
1
2
-
2
p
2
-
q
1
2
,
c
2
=
p
2
2
-
2
p
1
p
3
+
2
q
1
q
3
+
2
p
4

-

q
2
2
,
c
3
=
p
3
2
-
2
p
2
p
4
+
2
q
2
q
4
-
q
3
2
,
c
4
=
p
4
2
-
q
4
2
.
The conditions (
i
) and (
ii
) of Theorem 2 hold if and only if (19) has no real positive root.Letm
=
ω
2; then (19) takes the form
(21)
m
4
+
c
1
m
3
+
c
2
m
2
+
c
3
m
+
c
4
=
0
.
If c
4
<
0, then (19) has at least one positive root. In the case when (19) has four positive roots, we have
(22)
ω
1
=
m
1
,
ω
2
=
m
2
,
ω
3
=
m
3
,
ω
4
=
m
4
.
From (16), we have
(23)
τ
k
(
j
)
=
1
ω
k
{
arcsin
b
9
ω
k
7
+
b
10
ω
k
5
+
b
11
ω
k
3
+
b
12
ω
k
b
5
ω
k
6
+
b
6
ω
k
4
+
b
7
ω
k
2
+
b
8
+
2
j
π
}
,
where k
=
1,2
,
3,4 and j
=
0,1
,
2
,
…; we choose τ
0
=
min
⁡
(
τ
k
(
j
)
).To establish Hopf bifurcation atτ
=
τ
0, we need to show that
(24)
R
(
d
λ
d
τ
)
τ
=
τ
0
≠
0
.
By differentiating (10) with respect to τ, we can get
(25)
d
λ
d
τ
=
λ
e
-
λ
τ
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
×
(
e
-
λ
τ
(
4
λ
3
+
3
p
1
λ
2
+
2
p
2
λ
+
p
3
)
+
e
-
λ
τ
i
i
i
i
i
i
i
i
i
×
[
(
3
q
1
λ
2
+
2
q
2
λ
+
q
3
)
iiiiiiiiiiiiii
-
τ
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
]
e
-
λ
τ
)
-
1
.
It follows that
(26)
(
d
λ
d
τ
)
-
1
=
(
(
4
λ
3
+
3
p
1
λ
2
+
2
p
2
λ
+
p
3
)
+
e
-
λ
τ
i
i
i
i
i
i
×
[
(
3
q
1
λ
2
+
2
q
2
λ
+
q
3
)
i
i
i
i
i
i
i
i
i
i
-
τ
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
]
)
×
(
λ
e
-
λ
τ
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
)
-
1
.
Then, by combining (10), we get
(27)
(
d
λ
d
τ
)
-
1
=
(
(
4
λ
3
+
3
p
1
λ
2
+
2
p
2
λ
+
p
3
)
i
i
i
i
i
+
e
-
λ
τ
(
3
q
1
λ
2
+
2
q
2
λ
+
q
3
)
)
×
(
λ
e
-
λ
τ
(
q
1
λ
3
+
q
2
λ
2
+
q
3
λ
+
q
4
)
)
-
1
-
τ
λ
.
Substituting λ
=
i
ω
0 in (27) (where ω
0
>
0 and i
=
-
1) yields
(28)
(
d
λ
d
τ
)
-
1
|
τ
=
τ
0
=
d
1
+
i
d
2
d
3
+
i
d
4
-
τ
λ
,
where
(29)
d
1
=
(
p
3
-
3
p
1
ω
0
2
)
+
(
q
3
-
3
q
1
ω
0
2
)
cos
⁡
(
ω
0
τ
0
)
+
2
q
2
ω
0
sin
(
ω
0
τ
0
)
,
d
2
=
(
2
p
2
ω
0
-
4
ω
3
)
+
2
q
2
ω
0
cos
⁡
(
ω
0
τ
0
)
-
(
q
3
-
3
q
1
ω
0
2
)
sin
(
ω
0
τ
0
)
,
d
3

=

(
q
1
ω
0
4
-
q
3
ω
0
2
)
cos
⁡
(
ω
0
τ
0
)
+
(
q
4
ω
0
-
q
2
ω
0
3
)
sin
(
ω
0
τ
0
)
,
d
4

=

(
q
4
ω
0
-
q
2
ω
0
3
)
cos
⁡
(
ω
0
τ
0
)
-
(
q
1
ω
0
4
-
q
3
ω
0
2
)
sin
(
ω
0
τ
0
)
.
Thus,
(30)
R
(
d
λ
d
τ
)
-
1
|
τ
=
τ
0
=
d
1
d
3
+
d
2
d
4
d
3
2
+
d
4
2
.
Notice that
(31)
sign
⁡
(
R
d
λ
(
t
)
d
τ
)
|
τ
=
τ
0
=
sign
⁡
(
R
(
d
λ
d
τ
)
-
1
)
|
τ
=
τ
0
.
By summarizing the above analysis, we arrive at the following theorem.Theorem 3.
The infected equilibriumE
+ of the system (5) is asymptotically stable for τ
∈
[
0
,
τ
0
) and it undergoes Hopf bifurcation at τ
=
τ
0.


## 4. Direction and Stability of Bifurcating Periodic Solutions

In the previous section, we obtained conditions for Hopf bifurcation to occur whenτ
0
=
τ
k
(
j
), j
=
0,1
,
2
,
…. It is also important to derive explicit formulae from which we can determine the direction, stability, and period of periodic solutions bifurcating around the infected equilibrium E
+ at the critical value τ
0. We use the cafeteria of normal forms and center manifold proposed by Hassard [34]. We assume that the model (5) undergoes Hopf bifurcation at the infected equilibrium E
+ when τ
0
=
τ
k
(
j
), j
=
0,1
,
2
,
…, and then ±
i
ω
0 are the corresponding purely imaginary roots of the characteristic equation at the infected equilibrium E
+. Assume also that
(32)
(
X
1
(
t
)
,
X
2
(
t
)
,
X
3
(
t
)
,
X
4
(
t
)
)
T
=
(
z
*
(
t
)
x
(
t
)
-
x
*
,
y
(
t
)
-
y
*
(
t
)
,
i
i
i
i
i
i
i
i
w
(
t
)
-
w
*
(
t
)
,
z
(
t
)
-
z
*
(
t
)
)
T
;
then using Taylors expansion for system (3) at the equilibrium point yields
(33)
X
˙
1
=
k
11
X
1
(
t
)
+
k
12
X
2
(
t
)
+
k
13
X
1
(
t
)
X
1
(
t
)
+
k
14
X
1
(
t
)
X
2
(
t
)
,
X
˙
2
=
k
21
X
1
(
t
)
+
k
22
X
2
(
t
)
+
k
23
X
4
(
t
)
+
k
24
X
1
(
t
)
X
2
(
t
)
+
k
25
X
2
(
t
)
X
4
(
t
)
,
X
˙
3
=
k
31
X
3
(
t
)
+
k
32
X
2
(
t
-
τ
)
+
k
33
X
3
(
t
-
τ
)
+
k
34
X
2
(
t
-
τ
)
X
3
(
t
-
τ
)
,
X
˙
4
=
k
41
X
4
(
t
)
+
k
42
X
2
(
t
-
τ
)
+
k
43
X
3
(
t
-
τ
)
+
k
44
X
2
(
t
-
τ
)
X
3
(
t
-
τ
)
.
Here,
(34)
k
11
=
-
A
1
y
*
+
r
-
2
r
x
*
T
max
⁡
-
r
y
*
T
max
⁡
-
δ
1
,
k
12
=
-
A
1
x
*
-
r
x
*
T
max
⁡
,
k
13
=
-
2
r
T
max
⁡
,
k
14
=
-
r
T
max
⁡
-
A
1
,
k
21
=
A
1
y
*
,
k
22
=
A
1
x
*
-
A
2
-
p
z
*
,
k
23
=
-
p
y
*
,
k
24
=
A
1
,
k
25
=
-
p
,
k
31
=
-
b
,
k
32
=
c
(
1
-
q
)
w
*
,
k
33
=
c
(
1
-
q
)
y
*
,
k
34
=
c
(
1
-
q
)
,
k
41
=
-
h
,
k
42
=
c
q
w
*
,
k
43
=
c
q
y
*
,
k
44
=
c
q
.
For convenience, let τ
=
τ
0
+
μ and u
t
(
θ
)
=
u
(
t
+
θ
) for θ
∈
[
-
τ
,
0
]. Denote C
k
(
[
-
τ
,
0
]
,
R
4
)
=
{
ϕ
∣
ϕ
:
[
-
τ
,
0
]
→
R
4
}; ϕ has k-order continuous derivative. For initial conditions ϕ
(
θ
)
=
(
ϕ
1
(
θ
)
,
ϕ
2
(
θ
)
,
ϕ
3
(
θ
)
,
ϕ
4
(
θ
)
)
T
∈
C
(
[
-
τ
,
0
]
,
R
4
), (33) can be rewritten as
(35)
u
˙
(
t
)
=
L
μ
(
u
t
)
+
F
(
u
t
,
μ
)
,
where u
(
t
)
=
(
u
1
(
t
)
,
u
2
(
t
)
,
u
3
(
t
)
,
u
4
(
t
)
)
T
∈
C,  L
μ
:
C
→
R
4, and F
:
C
→
R
4 are given, respectively, by
(36)
L
μ
ϕ
=
(
τ
0
+
μ
)
G
1
ϕ
(
0
)
+
(
τ
0
+
μ
)
G
2
ϕ
(
-
τ
)
,
F
(
ϕ
,
μ
)
=
(
τ
0
+
μ
)
(
F
1
,
F
2
,
F
3
,
F
4
)
T
.
L
μ is one parameter family of bounded linear operators in C and
(37)
G
1
=
(
k
11
k
12
0
0
k
21
k
22
0
k
24
0
0
k
31
0
0
0
0
k
41
)
,
G
2
=
(
0
0
0
0
0
0
0
0
0
k
32
k
33
0
0
k
42
k
43
0
)
,
F
=
(
k
13
ϕ
1
(
0
)
ϕ
1
(
0
)
+
k
14
ϕ
1
(
0
)
ϕ
2
(
0
)
k
24
ϕ
1
(
0
)
ϕ
2
(
0
)
+
k
25
ϕ
2
(
0
)
ϕ
4
(
0
)
k
34
ϕ
2
(
-
τ
)
ϕ
3
(
-
τ
)
k
44
ϕ
2
(
-
τ
)
ϕ
3
(
-
τ
)
)
.
From the discussion in the above section, we know that if μ
=
0, then model (5) undergoes a Hopf bifurcation at the infected equilibrium E
+, and the associated characteristic equation of model (5) has a pair of purely imaginary roots ±
i
τ
0
ω
0. By Reisz representation, there exists a function η
(
θ
,
μ
) of bounded variation for θ
∈
[
-
τ
,
0
] such that
(38)
L
μ
ϕ
=
∫
-
τ
0
d
η
(
θ
,
μ
)
ϕ
(
θ
)
.
In fact, we can choose
(39)
η
(
θ
,
μ
)
=
(
τ
0
+
μ
)
G
1
δ
(
θ
)
+
(
τ
0
+
μ
)
G
2
δ
(
θ
+
τ
)
,
where δ
(
θ
) is Dirac delta function. Next, for ϕ
∈
C
1
(
[
-
τ
,
0
]
,
R
4
), define
(40)
A
(
μ
)
ϕ
=
{
d
ϕ
d
θ
,
θ
∈
[
-
τ
,
0
)
∫
-
τ
0
‍
d
η
(
θ
,
μ
)
ϕ
(
θ
)
,
θ
=
0
,
(41)
R
(
μ
)
ϕ
=
{
0
,
θ
=
[
-
τ
,
0
)
F
(
ϕ
,
μ
)
,
θ
=
0
.
Since u
˙
(
t
)
=
u
˙
t
(
θ
), (35) can be written as
(42)
u
˙
t
=
A
(
μ
)
u
t
+
R
(
μ
)
u
t
,
where u
t
=
u
(
t
+
θ
),  θ
∈
[
-
τ
,
0
]. For ψ
∈
C
1
(
[
0
,
τ
]
,
R
4
), the adjoint operator A
* of A can be defined as
(43)
A
*
ψ
(
s
)
ϕ
=
{
-
d
ψ
(
s
)
d
s
,
s
∈
(
-
τ
,
0
]
∫
-
τ
0
‍
d
η
(
θ
,
μ
)
ϕ
(
θ
)
,
s
=
0
.
For ϕ
∈
C
1
(
[
-
τ
,
0
]
,
R
4
) and ψ
∈
C
1
(
[
0
,
τ
]
,
R
4
), in order to normalize the eigenvalues of operator A and adjoint operator A
*, the following bilinear form is defined by
(44)
〈
ψ
,
ϕ
〉
=
ψ
-
(
0
)
ϕ
(
0
)
-
∫
θ
=
-
τ
0
‍
∫
ξ
=
0
θ
‍
ψ
-
(
ξ
-
θ
)
[
d
η
(
θ
)
]
ϕ
(
ξ
)
d
ξ
,
where η
(
θ
)
=
η
(
θ
,
0
) and ψ
- is complex conjugate of ψ. It can verify that A
* and A
(
0
) are adjoint operators with respect to this bilinear form.We assume that±
i
ω
0 are eigenvalues of A
(
0
) and the other eigenvalues have strictly negative real parts. Thus, they are also eigenvalues of A
*. Now we compute the eigenvector q of A corresponding to the eigenvalue i
ω
0 and the eigenvector q
* of A
* corresponding to the eigenvalue -
i
ω
0. Suppose that q
(
θ
)
=
(
1
,
p
1
,
p
2
,
p
3
)
T
e
i
ω
0
θ is eigenvector of A
(
0
) associated with i
ω
0; then, A
(
0
)
q
(
θ
)
=
i
ω
0
q
(
θ
). It follows from the definition of A
(
0
) and (36), (38), and (40) that(45)
(
k
11
-
i
ω
0
k
12
0
0
k
21
k
22
-
i
ω
0
0
k
23
0
k
32
e
-
i
ω
0
τ
0
k
31
+
k
33
e
-
i
ω
0
τ
0
-
i
ω
0
0
0
k
42
e
-
i
ω
0
τ
0
k
43
e
-
i
ω
0
τ
0
k
41
-
i
ω
0
)
q
(
0
)
=
(
0
0
0
0
)
.Solving (45), we can easily obtain q
(
0
)
=
(
1
,
p
1
,
p
2
,
p
3
)
T, where
(46)
p
1
=
i
ω
0
-
k
11
k
12
,
p
2
=
k
32
(
k
11
-
i
ω
0
)
e
-
i
ω
0
τ
0
k
12
(
k
31
+
k
33
e
-
i
ω
0
τ
0
-
i
ω
0
)
,
p
3
=
(
k
11
-
i
ω
0
)
(
k
22
-
i
ω
0
)
-
k
12
k
21
k
12
k
23
.
Similarly, suppose that the eigenvector q
* of A
* corresponding to -
i
ω
0 is q
*
(
s
)
=
(
1
/
D
)
(
1
,
p
1
*
,
p
2
*
,
p
3
*
)
T
e
i
ω
0
s, s
∈
[
0
,
τ
]. By the definition of A
* and (36), (38), and (40), one gets (47)
(
k
11
+
i
ω
0
k
21
0
0
k
12
k
22
+
i
ω
0
k
32
e
-
i
ω
0
τ
0
k
42
e
-
i
ω
0
τ
0
0
0
k
31
+
k
33
e
-
i
ω
0
τ
0
+
i
ω
0
k
43
e
-
i
ω
0
τ
0
0
k
23
0
k
41
+
i
ω
0
)
q
*
(
0
)
=
(
0
0
0
0
)
.Solving (47), we easily obtain q
*
(
0
)
=
(
1
/
D
)
(
1
,
p
1
*
,
p
2
*
,
p
3
*
)
T, where
(48)
p
1
*
=
-
k
11
+
i
ω
0
k
21
,
p
2
*
=
-
k
23
k
43
(
k
11
+
i
ω
0
)
e
-
i
ω
0
τ
0
k
21
(
k
41
+
i
ω
0
)
(
k
31
+
k
33
e
-
i
ω
0
τ
0
+
i
ω
0
)
,
p
3
*
=
k
23
(
k
11
+
i
ω
0
)
k
21
(
k
41
+
i
ω
0
)
.
In order to assure that 〈
q
*
,
q
〉
=
1, we need to determine the value of D. From (44), one gets
(49)
〈
q
*
,
q
〉
=
q
*
-
T
(
0
)
q
(
0
)
-
∫
θ
=
-
τ
0
0
‍
∫
ξ
=
0
θ
‍
q
*
-
T
(
ξ
-
θ
)
[
d
η
(
θ
)
]
q
(
ξ
)
d
(
ξ
)
=
1
D
-
(
1
+
p
1
p
1
-
*
+
p
2
p
2
-
*
+
p
3
p
3
-
*
)
-
∫
-
τ
0
0
∫
ξ
=
0
θ
1
D
-
(
1
,
p
1
-
*
,
p
2
-
*
p
3
-
*
)
e
-
i
ω
0
(
ξ
-
θ
)
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
×
[
d
η
(
θ
)
]
(
1
,
p
1
,
p
2
,
p
3
)
T
e
i
ω
0
ξ
d
ξ
=
1
D
-
(
1
+
p
1
p
1
-
*
+
p
2
p
2
-
*
+
p
3
p
3
-
*
)
-
∫
-
τ
0
0
1
D
-
(
1
,
p
1
-
*
,
p
2
-
*
,
p
3
-
*
)
θ
e
i
ω
0
θ
i
i
i
i
i
i
i
i
i
i
i
×
[
d
η
(
θ
)
]
(
1
,
p
1
,
p
2
,
p
3
)
T
=
1
D
-
(
(
1
,
p
1
,
p
2
,
p
3
)
T
(
1
+
p
1
p
1
-
*
+
p
2
p
2
-
*
+
p
3
p
3
-
*
)
i
i
i
i
i
i
i
i
i
i
+
τ
0
e
-
i
ω
0
τ
0
(
1
,
p
1
-
*
,
p
2
-
*
,
p
3
-
*
)
i
i
i
i
i
i
i
i
i
i
×
G
2
(
1
,
p
1
,
p
2
,
p
3
)
T
)
=
1
D
-
(
e
-
i
ω
0
τ
0
(
1
+
p
1
p
1
-
*
+
p
2
p
2
-
*
+
p
3
p
3
-
*
)
+
τ
0
e
-
i
ω
0
τ
0
i
i
i
i
i
i
i
i
i
i
×
(
(
k
32
p
2
-
*
+
k
42
p
3
-
*
)
p
1
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
+
(
k
33
p
2
-
*
+
k
43
p
3
-
*
)
p
2
)
e
-
i
ω
0
τ
0
)
;
D
-
=
(
1
+
p
1
p
1
-
*
+
p
2
p
2
-
*
+
p
3
p
3
-
*
)
+
τ
0
e
-
i
ω
0
τ
0
(
(
k
32
p
2
-
*
+
k
42
p
3
-
*
)
p
1
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
+
(
k
33
p
2
-
*
+
k
43
p
3
-
*
)
p
2
)
.
Let
(50)
v
(
t
)
=
〈
q
*
,
u
t
〉
,
W
(
t
,
θ
)
=
u
t
-
v
q
-
v
q
-
=
u
t
-
2
Re
(
v
(
t
)
q
(
θ
)
)
.
On the center manifold Ω
0, we have
(51)
W
(
t
,
θ
)
=
W
(
v
(
t
)
,
v
-
(
t
)
,
θ
)
,
where
(52)
W
(
v
,
v
-
,
θ
)
=
W
20
(
θ
)
v
2
2
+
W
11
(
θ
)
v
v
-
+
W
02
(
θ
)
v
-
v
2
2
+
⋯
.
v and v
- are local coordinates of the center manifold Ω
0 in the direction of q
* and q
-
*, respectively. Note that W is real if u
t is real. So we only consider real solutions. From (50), we obtain
(53)
〈
q
*
,
W
〉
=
〈
q
*
,
u
t
-
v
q
-
v
q
-
〉
=
〈
q
*
,
u
t
〉
-
v
(
t
)
〈
q
*
,
q
〉
-
v
-
(
t
)
〈
q
*
,
q
-
〉
.
For the solution u
t
∈
Ω
0 of (35), from (41) and (44), since μ
=
0, we have
(54)
v
˙
(
t
)
=
〈
q
*
,
u
˙
t
〉
=
〈
q
*
,
A
(
0
)
u
t
+
R
(
0
)
u
t
〉
=
〈
q
*
,
A
(
0
)
u
t
〉
+
〈
q
*
,
R
(
0
)
u
t
〉
=
〈
A
*
q
*
,
u
t
〉
+
q
*
-
T
(
0
)
F
(
u
t
,
0
)
=
i
ω
0
v
(
t
)
+
q
*
-
T
(
0
)
f
0
(
v
,
v
-
)
.
Rewrite (54) as
(55)
v
˙
(
t
)
=
i
ω
0
v
(
t
)
+
g
(
v
,
v
-
)
,
where
(56)
g
(
v
,
v
-
)
=
q
*
-
T
(
0
)
f
0
(
v
,
v
-
)
=
q
*
-
T
(
0
)
F
(
W
(
v
,
v
-
,
θ
)
+
2
Re
{
v
(
t
)
q
(
θ
)
,
0
}
)
=
g
20
v
2
2
+
g
11
v
v
-
+
g
02
v
-
2
2
+
g
21
v
2
v
-
2
⋯
.
Substituting (42) and (54) into (50) yields
(57)
W
˙
=
u
˙
(
t
)
-
v
˙
q
-
v
-
˙
q
-
=
A
u
t
+
R
u
t
-
(
i
ω
0
v
+
q
*
-
T
(
0
)
f
0
(
v
,
v
-
)
)
q
-
(
i
ω
0
v
-
+
q
*
-
T
(
0
)
f
0
(
v
,
v
-
)
¯
)
q
-
=
A
u
t
+
R
u
t
-
A
v
q
-
A
v
-
q
-
-
2
Re
(
q
*
-
T
(
0
)
f
0
(
v
,
v
-
)
q
)
,
(58)
W
˙

=

{

A
W

-

2
Re
(
q
*
-
T
(
0
)
f
0
(
v
,
v
-
)
q
)
,
θ

∈

[
-
τ
,
0
)

A
W

-

2
Re
(
q
*
-
T
(
0
)
f
0
(
v
,
v
-
)
q
)

+

f
0
(
v
,
v
-
)
,
θ

=

0
,
which can be written as
(59)
W
˙
=
A
W
+
H
(
v
,
v
-
,
θ
)
,
where
(60)
H
(
v
,
v
-
,
θ
)
=
H
20
(
θ
)
v
2
2
+
H
11
(
θ
)
v
v
-
+
H
02
(
θ
)
v
-
2
2
+
⋯
.
On the center manifold Ω
0, we have
(61)
W
˙
=
W
v
v
˙
+
W
v
-
v
-
˙
.
Substituting (52) and (55) into (61), one obtains
(62)
W
˙
=
(
W
20
v
+
W
11
v
-
+
⋯
)
(
i
ω
0
v
+
g
)
+
(
W
11
v
+
W
02
v
-
+
⋯
)
(
-
i
ω
0
v
-
+
g
-
)
.
Substituting (52) and (60) into (59) yields
(63)
W
˙
=
(
A
W
20
+
H
20
)
v
2
2
+
(
A
W
11
+
H
11
)
v
v
-
+
(
A
W
02
+
H
02
)
v
-
2
2
+
⋯
.
Comparing the coefficients of (62) and (63), one gets
(64)
(
A
-
i
2
ω
0
)
W
20
(
θ
)
=
-
H
20
(
θ
)
,
A
W
11
(
θ
)
=
-
H
11
(
θ
)
,
(
A
+
i
2
ω
0
)
W
02
(
θ
)
=
-
H
02
(
θ
)
.
Since u
t
=
u
(
t
+
θ
)
=
W
(
v
,
v
-
,
θ
)
+
v
q
+
v
-
q, then we have
(65)
u
t
=
(
u
1
(
t
+
θ
)
u
2
(
t
+
θ
)
u
3
(
t
+
θ
)
u
4
(
t
+
θ
)
)
=
(
W
(
1
)
(
v
,
v
-
,
θ
)
W
(
2
)
(
v
,
v
-
,
θ
)
W
(
3
)
(
v
,
v
-
,
θ
)
W
(
4
)
(
v
,
v
-
,
θ
)
)
+
v
(
1
p
1
p
2
p
3
)
e
i
ω
0
θ
+
v
-
(
1
p
-
1
p
-
2
p
-
3
)
e
-
i
ω
0
θ
.
Thus, we obtain
(66)
u
1
(
t
+
θ
)
=
W
(
1
)
(
v
,
v
-
,
θ
)
+
v
e
i
ω
0
θ
+
v
-
e
-
i
ω
0
θ
=
(
W
20
(
1
)
(
θ
)
v
2
2
+
W
11
(
1
)
(
θ
)
v
v
-
+

W
02
(
1
)
(
θ
)
v
-
2
2
+
⋯
)
+
v
e
i
ω
0
θ
+
v
-
e
-
i
ω
0
θ
,
u
2
(
t
+
θ
)
=
W
(
2
)
(
v
,
v
-
,
θ
)
+
v
p
1
e
i
ω
0
θ
+
v
-
p
-
1
e
-
i
ω
0
θ
=
(
W
20
(
2
)
(
θ
)
v
2
2
+
W
11
(
2
)
(
θ
)
v
v
-
+
W
02
(
2
)
(
θ
)
v
-
2
2
+
⋯
)
+
v
p
1
e
i
ω
0
θ
+
v
-
p
-
1
e
-
i
ω
0
θ
,
u
3
(
t
+
θ
)
=
W
(
3
)
(
v
,
v
-
,
θ
)
+
v
p
2
e
i
ω
0
θ
+
v
-
p
-
2
e
-
i
ω
0
θ
=
(

W
20
(
3
)
(
θ
)
v
2
2

+

W
11
(
3
)
(
θ
)
v
v
-

+

W
02
(
3
)
(
θ
)
v
-
2
2

+
⋯

)
+
v
p
2
e
i
ω
0
θ
+
v
-
p
-
2
e
-
i
ω
0
θ
,
u
4
(
t
+
θ
)
=
W
(
4
)
(
v
,
v
-
,
θ
)
+
v
p
3
e
i
ω
0
θ
+
v
-
p
-
3
e
-
i
ω
0
θ
=
(

W
20
(
4
)
(
θ
)
v
2
2

+

W
11
(
4
)
(
θ
)
v
v
-

+

W
02
(
4
)
(
θ
)
v
-
2
2

+
⋯

)
+
v
p
3
e
i
ω
0
θ
+
v
-
p
-
3
e
-
i
ω
0
θ
.
It is obvious that
(67)
ϕ
1
(
0
)
=
v
+
v
-
+
W
20
(
1
)
(
0
)
v
2
2
+
W
11
(
1
)
(
0
)
v
v
-
+
W
02
(
1
)
(
0
)
v
-
2
2
+
⋯
,
ϕ
2
(
0
)
=
v
p
1
+
v
-
p
-
1
+
W
20
(
2
)
(
0
)
v
2
2
+
W
11
(
2
)
(
0
)
v
v
-
+
W
02
(
2
)
(
0
)
v
-
2
2
+
⋯
,
ϕ
4
(
0
)
=
v
p
3
+
v
-
p
-
3
+
W
20
(
4
)
(
0
)
v
2
2
+
W
11
(
4
)
(
0
)
v
v
-
+
W
02
(
4
)
(
0
)
v
-
2
2
+
⋯
.
So
(68)
ϕ
1
(
0
)
ϕ
1
(
0
)
=
v
2
+
v
-
2
+
2
v
v
-
+
1
2
(
4
W
11
(
1
)
(
0
)
+
2
W
20
(
1
)
(
0
)
)
v
2
v
-
+
⋯
,
ϕ
1
(
0
)
ϕ
2
(
0
)
=
p
1
v
2
+
p
-
1
v
-
2
+
(
p
1
+
p
-
1
)
v
v
-
+
1
2
(
2
W
11
(
2
)
(
0
)
+
W
20
(
2
)
(
0
)
+
W
20
(
1
)
(
0
)
p
-
1
i
i
i
i
i
i
i
i
i
i
+
2
W
11
(
1
)
(
0
)
p
1
)
v
2
v
-
+
⋯
,
ϕ
2
(
0
)
ϕ
4
(
0
)
=
p
1
p
3
v
2
+
p
-
1
p
-
3
v
-
2
+
[
p
1
p
-
3
+
p
-
1
p
3
]
v
v
-
+
1
2
(
2
W
11
(
4
)
(
0
)
p
1
+
W
20
(
4
)
(
0
)
p
-
1
i
i
i
i
i
i
i
i
i
i
i
+
W
20
(
0
)
(
0
)
p
-
3
+
2
W
11
(
2
)
(
0
)
p
3
)
v
2
v
-
⋯
;
also
(69)
ϕ
2
(
-
τ
)
=
v
p
1
e
-
i
ω
0
τ
+
v
-
p
-
1
e
i
ω
0
τ
+
W
20
(
2
)
(
-
τ
)
v
2
2
+
W
11
(
2
)

(
-
τ
)
v
v
-
+
W
02
(
2
)

(
-
τ
)
v
-
2
2

+

⋯
,
ϕ
3
(
-
τ
)
=
v
p
2
e
-
i
ω
0
τ
+
v
-
p
-
2
e
i
ω
0
τ
+
W
20
(
3
)
(
-
τ
)
v
2
2
+
W
11
(
3
)
(
-
τ
)
v
v
-
+
W
02
(
3
)
(
-
τ
)
v
-
2
2
+
⋯
and hence
(70)
ϕ
2
(
-
τ
)
ϕ
3
(
-
τ
)
=
p
1
p
2
e
-
2
i
ω
0
τ
0
v
2
+
p
-
1
p
-
2
e
2
i
ω
0
τ
0
v
-
2
+
(
p
1
p
-
2
+
p
-
1
p
2
)
v
v
-
+
1
2
(
2
p
1
e
-
i
ω
0
τ
0
W
11
(
3
)

(
-
τ
)

+

p
-
1
e
i
ω
0
τ
W
20
(
3
)

(
-
τ
)

iiiiiiiiiii
+
2
p
2
e
-
i
ω
0
τ
0
W
11
(
2
)
(
-
τ
)
)
v
2
v
-
+
⋯
.
It follows from (54) that
(71)
f
0
(
v
,
v
-
)
=
(
k
13
ϕ
1
(
0
)
ϕ
1
(
0
)
+
k
14
ϕ
1
(
0
)
ϕ
2
(
0
)
k
24
ϕ
1
(
0
)
ϕ
2
(
0
)
+
k
25
ϕ
2
(
0
)
ϕ
4
(
0
)
k
34
ϕ
2
(
-
τ
)
ϕ
3
(
-
τ
)
k
44
ϕ
2
(
-
τ
)
ϕ
3
(
-
τ
)
)
=
(
F
11
v
2
+
F
12
v
-
2
+
F
13
v
v
-
+
F
14
v
2
v
-
F
21
v
2
+
F
22
v
-
2
+
F
23
v
v
-
+
F
24
v
2
v
-
F
31
v
2
+
F
32
v
-
2
+
F
33
v
v
-
+
F
34
v
2
v
-
F
41
v
2
+
F
42
v
-
2
+
F
43
v
v
-
+
F
44
v
2
v
-
)
,
where
(72)
F
11
=
k
13
+
k
14
p
1
,
F
12
=
k
13
+
k
14
p
-
1
,
F
13
=
2
k
13
+
k
14
(
p
1
+
p
-
1
)
,
F
14
=
k
13
(
2
W
11
(
1
)
(
0
)
+
W
20
(
1
)
(
0
)
)
F
14
+
1
2
k
14
(
2
W
11
(
2
)
(
0
)
+
W
20
(
2
)
(
0
)
iiii
iiiiiimiiiiiiii
+
W
20
(
1
)
(
0
)
p
-
1
+
2
W
11
(
1
)
(
0
)
p
1
)
,
F
21
=
k
24
p
1
+
k
25
p
1
p
3
,
F
22
=
k
24
p
-
1
+
k
25
p
-
1
p
-
3
,
F
23
=
k
24
(
p
1
+
p
-
1
)
+
k
25
(
p
1
p
-
3
+
p
-
1
p
3
)
,
F
24
=
1
2
k
24
(
2
W
11
(
2
)
(
0
)
+
W
20
(
2
)
(
0
)
+
W
20
(
1
)
(
0
)
p
-
1
i
i
i
i
i
i
n
m
i
i
i
i
i
i
+
2
W
11
(
2
)
(
0
)
p
1
)
F
24
+
1
2
k
25
(
2
W
11
(
4
)
(
0
)
p
1
+
W
20
(
4
)
(
0
)
p
-
1
+
W
20
(
2
)
(
0
)
p
-
3
i
i
i
i
i
i
i
i
i
m
n
i
i
i
i
i
i
+
2
W
11
(
2
)
(
0
)
p
3
)
,
F
31
=
k
34
(
p
1
p
2
e
-
2
i
ω
0
τ
0
)
,
F
32
=
k
34
(
p
-
1
p
-
2
e
2
i
ω
0
τ
0
)
,
F
33
=
k
34
(
p
1
p
-
2
+
p
-
1
p
2
)
,
F
34
=
1
2
k
34
(
2
p
1
e
-
i
ω
0
τ
0
W
11
(
3
)
(
-
τ
)
+
p
-
1
e
i
ω
0
τ
0
W
20
(
3
)
(
-
τ
)
i
i
i
i
i
i
i
i
m
i
i
i
i
i
i
+
2
p
2
e
-
i
ω
0
τ
0
W
11
(
2
)
(
-
τ
)
)
,
F
41
=
k
44
(
p
1
p
2
e
-
2
i
ω
0
τ
0
)
,
F
42
=
k
44
(
p
-
1
p
-
2
e
2
i
ω
0
τ
0
)
,
F
43
=
k
44
(
p
1
p
-
2
+
p
-
1
p
2
)
,
F
44
=
1
2
k
44
(
2
p
1
e
-
i
ω
0
τ
0
W
11
(
3
)
(
-
τ
)
+
p
-
1
e
i
ω
0
τ
0
W
20
(
3
)
(
-
τ
)
m
m
m
m
m
m
+
2
p
2
e
-
i
ω
0
τ
0
W
11
(
2
)
(
-
τ
)
)
.
Since q
-
*
(
0
)
=
(
1
/
D
-
)
(
1
,
p
-
1
*
,
p
-
2
*
,
p
-
3
*
)
T, we have
(73)
g
(
v
,
v
-
)
=
q
-
*
(
0
)
T
f
0
(
v
,
v
-
)
=
1
D
-
(
1
,
p
-
1
*
,
p
-
2
*
,
p
-
3
*
)
×
(
F
11
v
2
+
F
12
v
-
2
+
F
13
v
v
-
+
F
14
v
2
v
-
F
21
v
2
+
F
22
v
-
2
+
F
23
v
v
-
+
F
24
v
2
v
-
F
31
v
2
+
F
32
v
-
2
+
F
33
v
v
-
+
F
34
v
2
v
-
F
41
v
2
+
F
42
v
-
2
+
F
43
v
v
-
+
F
44
v
2
v
-
)
=
1
D
-
(
v
2
v
-
(
F
11
+
F
21
p
-
1
*
+
F
31
p
-
2
*
+
F
41
p
-
3
*
)
v
2
iiiiiiiii
+
(
F
12
+
F
22
p
-
1
*
+
F
32
p
-
2
*
+
F
42
p
-
3
*
)
v
-
2
iiiiiiiii
+
(
F
13
+
F
23
p
-
1
*
+
F
33
p
-
2
*
+
F
43
p
-
3
*
)
v
v
-
iiiiiiii
+
(
F
14
+
F
24
p
-
1
*
+
F
34
p
-
2
*
+
F
44
p
-
3
*
)
v
2
v
-
)
.
Comparing the coefficients of the above equation with those in (61), we have
(74)
g
20
=
2
D
-
(
F
11
+
F
21
p
-
1
*
+
F
31
p
-
2
*
+
F
41
p
-
3
*
)
,
g
11
=
1
D
-
(
F
13
+
F
23
p
-
1
*
+
F
33
p
-
2
*
+
F
43
p
-
3
*
)
,
g
02
=
2
D
-
(
F
12
+
F
22
p
-
1
*
+
F
32
p
-
2
*
+
F
42
p
-
3
*
)
,
g
21
=
2
D
-
(
F
14
+
F
24
p
-
1
*
+
F
34
p
-
2
*
+
F
44
p
-
3
*
)
.
We need to compute W
20
(
θ
) and W
11
(
θ
) for θ
∈
[
-
τ
,
0
). Equations (62) and (63) imply that
(75)
H
(
v
,
v
-
,
θ
)
=
-
2
Re
{
q
*
-
T
(
0
)
f
0
(
v
,
v
-
)
q
(
θ
)
}
=
-
2
Re
{
g
(
v
,
v
-
)
q
(
θ
)
}
=
-
g
(
v
,
v
-
)
q
(
θ
)
-
g
-
(
v
,
v
-
)
q
-
(
θ
)
,
H
(
v
,
v
-
,
θ
)
=
-
(
g
20
v
2
2
+
g
11
v
v
-
+
g
02
v
-
2
2
+
g
21
v
2
v
-
2
⋯
)
q
(
θ
)
-

(

g
-
20
v
-
2
2
+
g
-
11
v
v
-
+
g
-
02
v
2
2
+
g
-
21
v
-
2
v
2
⋯

)
q
-
(
θ
)
.
Comparing the coefficients of the above equation with (60), we have
(76)
H
20
(
θ
)
=
-
g
20
q
(
θ
)
-
g
-
02
q
-
(
θ
)
,
H
11
(
θ
)
=
-
g
11
q
(
θ
)
-
g
-
11
q
-
(
θ
)
,
H
02
(
θ
)
=
-
g
02
q
(
θ
)
-
g
-
20
q
-
(
θ
)
.
It follows from (40) and (64) that
(77)
W
˙
(
θ
)
=
A
W
20
=
2
i
ω
0
W
20
(
θ
)
-
H
20
(
θ
)
=
2
i
ω
0
W
20
(
θ
)
+
g
20
q
(
0
)
e
i
ω
0
θ
+
g
-
02
q
-
(
0
)
e
-
i
ω
0
θ
.
By solving the above equation for W
20
(
θ
) and for W
11
(
θ
), one obtains
(78)
W
20
(
θ
)
=
i
g
20
ω
0
q
(
0
)
e
i
ω
0
θ
+
i
g
-
02
3
ω
0
q
-
(
0
)
e
-
i
ω
0
θ
+
E
1
e
2
i
ω
0
θ
,
W
11
(
θ
)
=
-
i
g
11
ω
0
q
(
0
)
e
i
ω
0
θ
+
i
g
-
11
ω
0
q
-
(
0
)
e
-
i
ω
0
θ
+
E
2
,
where E
1 and E
2 can be determined by setting θ
=
0 in H
(
v
,
v
-
,
θ
).In fact, we have(79)
H
(
v
,
v
-
,
0
)
=
-
2
Re
{
q
*
-
T
(
0
)
f
0
(
v
,
v
-
q
)
}
+
f
0
(
v
,
v
-
)
=
-
(
g
20
v
2
2
+
g
11
v
v
-
+
g
02
v
-
2
2
+
g
21
v
2
v
-
2
⋯
)
q
(
0
)
-

(

g
-
20
v
-
2
2
+
g
-
11
v
v
-
+
g
-
02
v
2
2
+
g
-
20
v
-
2
v
2
+
⋯

)

q
-
(
0
)
+
(
F
11
v
2
+
F
12
v
-
2
+
F
13
v
v
-
+
F
14
v
2
v
-
F
21
v
2
+
F
22
v
-
2
+
F
23
v
v
-
+
F
24
v
2
v
-
F
31
v
2
+
F
32
v
-
2
+
F
33
v
v
-
+
F
34
v
2
v
-
F
41
v
2
+
F
42
v
-
2
+
F
43
v
v
-
+
F
44
v
2
v
-
)
;
comparing the coefficients of the above equations with those in (61), it follows that
(80)
H
20
(
0
)
=
-
g
20
q
(
0
)
-
g
-
02
q
-
(
0
)
+
(
F
11
,
F
21
,
F
31
,
F
41
)
T
,
H
11
(
0
)
=
-
g
11
q
(
0
)
-
g
-
11
q
-
(
0
)
+
(
F
13
,
F
23
,
F
33
,
F
43
)
T
.
By the definition of A and (40) and (64), we get
(81)
∫
-
τ
0
0
d
η
(
θ
)
W
20
(
θ
)
=
A
W
20
(
0
)
=
2
i
ω
0
W
20
(
0
)
-
H
20
(
0
)
,
∫
-
τ
0
0
d
η
(
θ
)
W
11
(
θ
)
=
A
W
11
(
0
)
=
-
H
11
(
0
)
.
One can notice that
(82)
(
i
ω
0
I
-
∫
-
τ
0
0
‍
e
i
ω
0
θ
d
η
(
θ
)
)
q
(
0
)
=
0
,
(
-
i
ω
0
I
-
∫
-
τ
0
0
‍
e
-
i
ω
0
θ
d
η
(
θ
)
)
q
-
(
0
)
=
0
.
Thus, we obtain
(83)
(
2
i
ω
0
I
-
∫
-
τ
0
0
e
2
i
w
0
θ
d
η
(
θ
)
)
E
1
=
(
F
11
,
F
21
,
F
31
,
F
41
)
T
(
∫
-
τ
0
0
d
η
(
θ
)
)
E
2
=
-
(
F
13
,
F
23
,
F
33
,
F
43
)
T
,
where E
1
=
(
E
1
(
1
)
,
E
1
(
2
)
,
E
1
(
3
)
,
E
1
(
4
)
)
T, E
2
=
(
E
2
(
1
)
,
E
2
(
2
)
,
E
2
(
3
)
,
E
2
(
4
)
)
T; the above equation can be written as(84)
(
2
i
ω
0
-
k
11
-
k
12
0
0
-
k
21
2
i
ω
0
-
k
22
0
-
k
23
0
-
k
32
e
-
i
w
0
τ
0
2
i
ω
0
-
k
31
-
k
33
e
-
i
ω
0
τ
0
0
0
-
k
42
e
-
i
ω
0
τ
0
-
k
43
e
-
i
ω
0
τ
0
2
i
ω
0
-
k
41
)
E
1
=
(
F
11
F
21
F
31
F
41
)
,
(
k
11
k
12
0
0
k
21
k
22
0
k
23
0
k
32
k
31
0
0
k
42
k
43
k
41
)
E
2
=
(
F
13
F
23
F
33
F
43
)
.From (78), (84), we can calculate g
21, and we can derive the following parameters:
(85)
C
1
(
0
)
=
i
2
ω
0
(
1
3
g
20
g
11
-
2
|
g
11
|
2
-
1
3
|
g
02
|
2
)
+
g
21
2
,
μ
2
=
-
Re
(
C
1
(
0
)
)
Re
(
λ
′
(
τ
0
)
)
,
β
2
=
2
Re
C
1
(
0
)
,
T
2
=
-
Im
⁡
{
C
1
(
0
)
}
+
μ
2
Im
⁡
λ
(
τ
0
)
ω
0
.
We arrive at the following theorem.Theorem 4.
The periodic solution is supercritical (subcritical) ifμ
2
>
0  (
μ
2
<
0
); the bifurcating periodic solutions are orbitally asymptotically stable with asymptotical phase (unstable) if β
2
<
0  (
β
2
>
0
); the period of the bifurcating periodic solution increases (decreases) if T
2
>
0  (
T
2
<
0
).


## 5. Numerical Simulations

In this section, we provide some simulations of model (4) to exhibit the impact of discrete time delay in the model. We consider the parameters values: Λ
=
10, δ
1
=
0.06, δ
2
=
0.3, e
1
=
0.2, β
=
0.1, p
=
1, c
=
0.1, b
=
0.02, q
=
0.02, η
∈
[
0,1
], h
=
0.1, r
=
0.03, ϵ
∈
[
0,1
], and T
max
⁡
=
1500. According to the given parameters’ values, the threshold critical value τ
0
=
0.4957 from the formula (21) exists. The steady state E
+ exists and is asymptotically stable (see Figure 1). We may notice that the solution converges to the equilibrium E
+ with damping oscillations as the value of τ increases. Once the delay τ crosses the critical value τ
0, then the model shows the existence of Hopf bifurcation which is depicted in the Figure 2. In Figure 3, we consider the efficacy of antiretroviral value is 0.9, which may be responsible for the loss of stability. The asymptotic behavior to the infection-free steady state, when we consider antiviral treatment (with ϵ
=
0.9,  η
=
0.9, and time delay τ
=
15), is shown in Figure 4. According to Theorem 4, the parameters C
1
=
-
2.1108
e
+
004
+
1.1224
e
+
005
i, λ
′
=
-
12.1371
-
0.6438
i, μ
2
=
-
1.7391
e
+
003, β
2
=
-
4.2215
e
+
004, and T
2
=
-
2.8052
e
+
005 are estimated. Based on these values one can conclude that bifurcating periodic solutions are unstable and decreases in the period of bifurcating periodic solutions. The existence of periodic solution is subcritical. For numerical treatment of DDEs and related issues; we refer the readers to [35, 36].Each panel (from (a) to (h)) shows the time evolution and trajectory of model (4) when τ
(
=
0.4
)
<
τ
0 (critical value) and the effect of therapies is considered to be ϵ
=
0.9 and η
=
0.2. It shows that the endemic steady state E
+ of model is asymptotically stable.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)It shows the numerical simulations of model (4), when the time delay of immune activation exceeds the critical value, τ
=
0.5
>
τ
0. The endemic steady state E
+ of the model undergoes Hopf bifurcation; stability switch and periodic solutions appear.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)It shows the numerical simulations of model (4), when the efficacy rate of antiretroviral treatments is considered to be low; that is, ϵ
=
0.2 and η
=
0.2. It shows that the equilibrium E
+ of the model undergoes Hopf bifurcation with oscillatory behavior in solutions even though the delay value is less than the critical value (
τ
=
0.4
<
τ
0
).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Several packages and types of software are available for the numerical integration and/or the study of bifurcations in delay differential equations (see, e.g., [37, 38]. In this paper we utilizeMIDDE code [39]) which is suitable to simulate stiff and nonstiff delay differential equations and Volterra delay integrodifferential equations, using monoimplicit RK methods.


## 6. Concluding Remarks

In this manuscript, we provided a conceptualCD
4
+ T-cell infection model which includes the logistic growth term along with two different types of antiretroviral drug therapies. The model includes a discrete time delay in the immune activation response, which plays an important role in the dynamics of the model. The infection-free and endemic steady states of the model are determined (Figure 5). The stability of steady states is analyzed. We deduced a formula that determines the critical value (branch value) τ
0. Necessary and sufficient conditions for the equilibrium to be asymptotically stable for all positive delay values are proved. We have seen that if the time delay exceeds the critical value τ
0, model (4) undergoes a Hopf bifurcation. The direction and stability of bifurcating periodic solutions are deduced in explicit formulae, using center manifold and normal forms. We also presented some numerical simulations to the underlying model to investigate the obtained results and theory. We have seen also that the antiretroviral treatments help to increase the level of uninfected CD
4
+ T-cells. The theoretical results that were confirmed by the numerical simulations show that the delayed CTL response can lead to complex bifurcations, and, in particular, the coexistence of multiple stable periodic solutions. When the time delay exceeds the critical (threshold) value, we may get subcritical behaviour that leads to a loss of uninfected CD
4
+ T-cells.It shows the numerical simulations model (4) when the efficacy rate of antiretroviral treatment is at expected level, ϵ
=
0.9 and η
=
0.9, and the delay value exceeds the critical value τ
=
15
>
τ
0. The solution always lies within the feasible region and the infection-free steady state E
0 is asymptotically stable.
(a)
(b)
(c)
(d)



---

*Source: 838396-2014-08-31.xml*