# Global Stability of Humoral Immunity HIV Infection Models with Chronically Infected Cells and Discrete Delays

Authors: A. M. Elaiw; N. A. Alghamdi
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

We study the global stability of three HIV infection models with humoral immune response. We consider two types of infected cells: the first type is the short-lived infected cells and the second one is the long-lived chronically infected cells. In the three HIV infection models, we modeled the incidence rate by bilinear, saturation, and general forms. The models take into account two types of discrete-time delays to describe the time between the virus entering into an uninfected CD4+ T cell and the emission of new active viruses. The existence and stability of all equilibria are completely established by two bifurcation parameters, R
0 and R
1. The global asymptotic stability of the steady states has been proven using Lyapunov method. In case of the general incidence rate, we have presented a set of sufficient conditions which guarantee the global stability of model. We have presented an example and performed numerical simulations to confirm our theoretical results.


---

## Body

## 1. Introduction

During last decades, many researchers have developed and analyzed several mathematical models which describe human immunodeficiency virus (HIV) dynamics (see, e.g., [1–12]). HIV mainly targets the CD4+ T cells, leading to Acquired Immunodeficiency Syndrome (AIDS). Most of the HIV mathematical models presented in the literature consider only one type of infected cells called short-lived infected cells. However, it was shown that there is another source for the virus which is called long-lived chronically infected cells. This type of cells generates smaller number of viruses than the short-lived infected cells, but it lives longer [4]. The basic HIV dynamics model with long-lived chronically infected cells presented in [4] is given by(1)
T
˙
=
π
-
1
-
ε
k
T
V
-
d
T
,
T
˙
∗
=
1
-
ξ
1
-
ε
k
T
V
-
b
T
∗
,
C
˙
∗
=
ξ
1
-
ε
k
T
V
-
a
C
∗
,
V
˙
=
N
T
b
T
∗
+
N
C
a
C
∗
-
c
V
.Here T
,
T
∗
,
C
∗, and V are the concentrations of the uninfected CD4+ T cells, short-lived infected CD4+ T cells, long-lived chronically infected CD4+ T cells, and free virus particles, respectively. π represents birth rate constant of the uninfected CD4+ T cells. k is the infection rate constant. Parameters d
,
b
,
a, and c are the death rate constants of uninfected CD4+ T cells, short-lived infected CD4+ T cells, long-lived chronically infected CD4+ T cells, and free viruses, respectively. The fractions ξ and (
1
-
ξ
) with 0
<
ξ
<
1 are the probabilities that, upon infection, an uninfected CD4+ T cell will become either long-lived chronically infected or short-lived infected. N
T and N
C denote the average numbers of free virus particles produced in the lifetime of the short-lived infected and long-lived chronically infected cells, respectively. Model (1) incorporates reverse transcriptase inhibitor drugs with drug efficacy ε and 0
≤
ε
<
1.In model (1), the immune response has not been modeled. The immune response plays an important role in controlling the diseases. In reality, the immune response needs indispensable components to do its job such as antibodies, cytokines, natural killer cells, and T cells. The antibody immune response is a part of the adaptive system in which the body responds to pathogens by primarily using the antibodies which are generated by the B cells, while the other part is the Cytotoxic T Lymphocytes (CTL) immune response where the CTL attacks and kills the infected cells [3]. In malaria disease, the humoral immune response is more effective than the CTL immune response. In the virus dynamics literature, several models have considered the effect of CTL immune response [3, 13] or the humoral immune response [14–16]. Obaid and Elaiw [15] proposed the following model which takes into consideration the humoral immune response:(2)
T
˙
=
π
-
1
-
ε
k
T
V
-
d
T
,
T
˙
∗
=
1
-
ξ
1
-
ε
k
T
V
-
b
T
∗
,
C
˙
∗
=
ξ
1
-
ε
k
T
V
-
a
C
∗
,
V
˙
=
N
T
b
T
∗
+
N
C
a
C
∗
-
r
V
Z
-
c
V
,
Z
˙
=
g
V
Z
-
μ
Z
.Here, the variable Z represents the concentration of B cells. The HIV are attacked by the B cells at rate r
V
Z. The terms g
V
Z and μ
Z represent the proliferation and death rates of the B cells, respectively. In model (2), it is assumed that once the HIV contacts the CD4+ T cell, it becomes infected producing new viruses. Actually, there exists an intracellular time delay between the time the HIV contacts an uninfected CD4+ T cell and the time it becomes actively infected CD4+ T cell [17]. In the literature, several papers have proposed various HIV models with time delays [17–21].Our aim in this work is to propose three HIV dynamics models with two types of infected cells, two types of intracellular delays, and humoral immunity. Bilinear and saturated incidences have been proposed in the first and second model, respectively, while a general nonlinear incidence rate is proposed in the third model. For each model, we derive two bifurcation parameters,R
0 and R
1, and establish the global stability using Lyapunov functional.


## 2. Model with Bilinear Incidence

We propose the following HIV infection model with humoral immunity, two types of infected cells and two types of intracellular delays:(3)
T
˙
t
=
π
-
1
-
ε
k
T
t
V
t
-
d
T
t
,
(4)
T
˙
∗
t
=
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
b
T
∗
t
,
(5)
C
˙
∗
t
=
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
a
C
∗
t
,
(6)
V
˙
t
=
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
t
Z
t
-
c
V
t
,
(7)
Z
˙
t
=
g
V
t
Z
t
-
μ
Z
t
.Parameter τ
1 represents the time between HIV contact with an uninfected CD4+ T cell and the cell becoming infected but not yet producer cell. The factor e
-
m
1
τ
1 represents the loss of CD4+ T cells during the interval [
t
-
τ
1
,
t
]. The parameters τ
2 and τ
3 represent the time necessary for producing new infectious viruses from the short-lived infected and long-lived chronically infected cells, respectively. The factors e
-
m
2
τ
2 and e
-
m
3
τ
3 represent the loss of short-lived and long-lived chronically infected cells during the intervals [
t
-
τ
2
,
t
] and [
t
-
τ
3
,
t
], respectively. Here, m
1
,
m
2, and m
3 are positive constants.The initial conditions for system (2) are given by(8)
T
u
=
φ
1
u
,
T
∗
u
=
φ
2
u
,
C
∗
u
=
φ
3
u
,
V
u
=
φ
4
u
,
Z
u
=
φ
5
u
,
φ
j
u
⩾
0
,
u
∈
-
τ
,
0
,
j
=
1,2
,
…
,
5
,
φ
j
0
>
0
,
j
=
1,2
,
…
,
5
,where φ
1
u
,
φ
2
u
,
…
,
φ
5
u
∈
C
(
[
-
τ
.
0
)
,
R
+
5
), where C is the Banach space of continuous functions and τ
=
max
⁡
τ
1
,
τ
2
,
τ
3. We note that system (3)–(7) with initial conditions (8) has a unique solution satisfying t [22].


### 2.1. Positive Invariance

Proposition 1.
The solution(
T
(
t
)
,
T
∗
(
t
)
,
C
∗
(
t
)
,
V
(
t
)
,
Z
(
t
)
) of (3)–(7) with initial conditions (8) is nonnegative for t
⩾
0 and ultimately bounded.Proof.
Assume thatT
(
t
)
≤
0 on [
0
,
ρ
] for some constant ρ and T
(
t
∗
)
=
0, t
∗
∈
[
0
,
ρ
]. From (3), we get T
˙
(
t
∗
)
=
π
>
0 and hence T
(
t
)
>
0, for some t
∈
(
t
∗
,
t
∗
+
ϱ
), and sufficiently small ϱ
>
0. This leads to a contradiction; therefore, T
(
t
)
>
0, for all t
⩾
0. From (4), (5), and (6), we have(9)
T
∗
t
=
φ
2
0
e
-
b
t
+
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
∫
0
t
e
-
b
t
-
η
T
η
-
τ
1
V
η
-
τ
1
d
η
,
C
∗
t
=
φ
3
0
e
-
a
t
+
e
-
m
1
τ
1
ξ
1
-
ε
k
∫
0
t
e
-
a
t
-
η
T
η
-
τ
1
V
η
-
τ
1
d
η
,
V
t
=
φ
4
0
e
-
∫
0
t
c
+
r
Z
ς
d
ς
+
e
-
m
2
τ
2
N
T
b
∫
0
t
e
-
∫
η
t
c
+
r
Z
ς
d
ς
T
∗
η
-
τ
2
d
η
+
e
-
m
3
τ
3
N
C
a
∫
0
t
e
-
∫
η
t
c
+
r
Z
ς
d
ς
C
∗
η
-
τ
3
d
η
.This confirms that T
∗
(
t
)
⩾
0, C
∗
(
t
)
⩾
0, and V
(
t
)
⩾
0, for all t
∈
[
0
,
τ
]. By a recursive argument, we get that T
∗
(
t
)
⩾
0, C
∗
(
t
)
⩾
0, and V
(
t
)
⩾
0, for all t
⩾
0. Moreover, from (7), we obtain(10)
Z
t
=
φ
5
0
e
-
∫
0
t
μ
-
g
V
ς
d
ς
. Clearly, Z
(
t
)
⩾
0, t
⩾
0. Now, we let G
1
(
t
)
=
e
-
m
1
τ
1
T
(
t
-
τ
1
)
+
T
∗
(
t
)
+
C
∗
(
t
), G
2
(
t
)
=
V
(
t
)
+
r
/
g
Z
(
t
), s
1
=
min
⁡
d
,
a
,
b, s
2
=
min
⁡
c
,
μ, L
1
=
π
/
s
1, L
2
=
(
N
T
b
+
N
C
a
)
L
1
/
s
2, and L
3
=
g
L
2
/
r; then,(11)
G
˙
1
t
=
e
-
m
1
τ
1
π
-
d
T
t
-
τ
1
-
b
T
∗
t
-
a
C
∗
t
≤
π
e
-
m
1
τ
1
-
s
1
e
-
m
1
τ
1
T
t
-
τ
1
+
T
∗
t
+
C
∗
t
=
π
e
-
m
1
τ
1
-
s
1
G
1
t
<
π
-
s
1
G
1
t
.Hence, l
i
m
s
u
p
t
→
∞
G
1
(
t
)
≤
L
1. Since T
(
t
)
>
0, T
∗
(
t
)
≥
0, and C
∗
(
t
)
≥
0, then l
i
m
s
u
p
t
→
∞
T
∗
(
t
)
≤
L
1 and l
i
m
s
u
p
t
→
∞
C
∗
(
t
)
≤
L
1. Also,(12)
G
˙
2
t
=
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
c
V
t
-
r
μ
g
Z
t
≤
e
-
m
2
τ
2
N
T
b
+
e
-
m
3
τ
3
N
C
a
L
1
-
s
2
V
t
+
r
g
Z
t
<
N
T
b
+
N
C
a
L
1
-
s
2
G
2
t
.Hence, l
i
m
s
u
p
t
→
∞
G
2
(
t
)
≤
L
2. Since V
(
t
) and Z
(
t
) are nonnegative, then l
i
m
s
u
p
t
→
∞
V
(
t
)
≤
L
2 and l
i
m
s
u
p
t
→
∞
Z
(
t
)
≤
L
3. Therefore, T
(
t
)
,
T
∗
(
t
)
,
C
∗
(
t
)
,
V
(
t
), and Z
(
t
) are ultimately bounded.


### 2.2. Steady States

System (3)–(7) always admits an uninfected steady state S
0
=
(
π
/
d
,
0,0
,
0,0
). Let T
0
=
π
/
d
. Now, we define the basic reproduction number for system (3)–(7) as(13)
R
0
=
1
-
ε
k
T
0
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
c
. The system has two other steady states, infected steady state without humoral immune response S
1
=
(
T
1
,
T
1
∗
,
C
1
∗
,
V
1
,
0
) and infected steady state with humoral immune response S
2
=
(
T
2
,
T
2
∗
,
C
2
∗
,
V
2
,
Z
2
):(14)
T
1
=
T
0
R
0
,
T
1
∗
=
e
-
m
1
τ
1
1
-
ξ
π
b
R
0
-
1
R
0
,
C
1
∗
=
e
-
m
1
τ
1
ξ
π
a
R
0
-
1
R
0
,
V
1
=
d
1
-
ε
k
R
0
-
1
,
T
2
=
π
g
g
d
+
1
-
ε
k
μ
,
T
2
∗
=
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
π
μ
b
d
g
+
1
-
ε
k
μ
,
C
2
∗
=
e
-
m
1
τ
1
ξ
1
-
ε
k
π
μ
a
d
g
+
1
-
ε
k
μ
,
V
2
=
μ
g
,
Z
2
=
c
r
d
g
R
0
d
g
+
1
-
ε
k
μ
-
1
.We note that T
1, T
1
∗, C
1
∗, and V
1 are positive when R
0
>
1 and where T
2, T
2
∗, C
2
∗, and V
2
>
0 and Z
2
>
0 when d
g
R
0
/
d
g
+
(
1
-
ε
)
k
μ
>
1. Now, we define humoral immune response reproduction number R
1 as(15)
R
1
=
R
0
1
+
1
-
ε
k
μ
/
d
g
.Clearly, R
0
>
R
1. From above, we can state the following lemma.Lemma 2.
For system (3)–(7), one has the following:(i)
IfR
0
≤
1, then the system has only one positive steady state S
0.
(ii)
IfR
1
≤
1
<
R
0, then the system has two positive steady states S
0 and S
1.
(iii)
IfR
1
>
1, then the system has three positive steady states S
0, S
1, and S
2.


### 2.3. Global Stability Analysis

We establish the global stability of all the steady states of system (3)–(7) employing the method of Lyapunov function. Let us define(16)
H
θ
=
θ
-
1
-
ln
⁡
θ
.Theorem 3.
For system (3)–(7), if R
0
≤
1, then S
0 is GAS.Proof.
Define(17)
U
0
=
T
0
H
T
T
0
+
η
1
T
∗
+
η
2
C
∗
+
η
3
V
+
η
4
Z
+
η
5
∫
0
τ
1
T
t
-
θ
V
t
-
θ
d
θ
+
η
6
∫
0
τ
2
T
∗
t
-
θ
d
θ
+
η
7
∫
0
τ
3
C
∗
t
-
θ
d
θ
,where η
i
,
i
=
1
,
…
,
7, are positive constants satisfying the following equations:(18)
e
-
m
1
τ
1
1
-
ε
k
1
-
ξ
η
1
+
ξ
η
2
=
η
5
,
η
6
-
b
η
1
=
0
,
η
7
-
a
η
2
=
0
,
η
6
-
η
3
e
-
m
2
τ
2
N
T
b
=
0
,
η
7
-
η
3
e
-
m
3
τ
3
N
C
a
=
0
,
g
η
4
-
r
η
3
=
0
,
η
5
=
1
-
ε
k
.The solution of (18) is given by(19)
η
1
=
N
T
e
m
1
τ
1
+
m
3
τ
3
e
m
3
τ
3
1
-
ξ
N
T
+
e
m
2
τ
2
ξ
N
C
,
η
2
=
N
C
e
m
1
τ
1
+
m
2
τ
2
e
m
3
τ
3
1
-
ξ
N
T
+
e
m
2
τ
2
ξ
N
C
,
η
3
=
e
m
1
τ
1
+
m
2
τ
2
+
m
3
τ
3
e
m
3
τ
3
1
-
ξ
N
T
+
e
m
2
τ
2
ξ
N
C
,
η
4
=
r
e
m
1
τ
1
+
m
2
τ
2
+
m
3
τ
3
g
e
m
3
τ
3
1
-
ξ
N
T
+
e
m
2
τ
2
ξ
N
C
,
η
5
=
1
-
ε
k
,
η
6
=
b
N
T
e
m
1
τ
1
+
m
3
τ
3
e
m
3
τ
3
1
-
ξ
N
T
+
e
m
2
τ
2
ξ
N
C
,
η
7
=
a
N
C
e
m
1
τ
1
+
m
2
τ
2
e
m
3
τ
3
1
-
ξ
N
T
+
e
m
2
τ
2
ξ
N
C
.The values of η
i
,
i
=
1
,
…
,
7, will be used through the paper. Calculating the derivative of U
0 along the solutions of system (3)–(7) and applying π
=
T
0
d, we obtain(20)
d
U
0
d
t
=
1
-
T
0
T
π
-
1
-
ε
k
T
V
-
d
T
+
η
1
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
b
T
∗
+
η
2
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
a
C
∗
+
η
3
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
g
V
Z
-
μ
Z
+
η
5
T
V
-
η
5
T
t
-
τ
1
V
t
-
τ
1
+
η
6
T
∗
-
η
6
T
∗
t
-
τ
2
+
η
7
C
∗
-
η
7
C
∗
t
-
τ
3
.Then,(21)
d
U
0
d
t
=
1
-
T
0
T
π
-
d
T
+
1
-
ε
k
T
0
V
-
η
3
c
V
-
η
4
μ
Z
=
-
d
T
-
T
0
2
T
+
c
η
3
R
0
-
1
V
-
η
4
μ
Z
.Therefore, if R
0
≤
1, then d
U
0
/
d
t
≤
0, for all T
,
V
,
Z
>
0. The solutions of system (3)–(7) limit M, the largest invariant subset of d
U
0
/
d
t
=
0 [22]. We note that d
U
0
/
d
t
=
0 if and only if T
(
t
)
=
T
0, V
(
t
)
=
0, and Z
(
t
)
=
0. For each element of M, we have V
(
t
)
=
0 and Z
(
t
)
=
0,; then, V
˙
(
t
)
=
0 and(22)
0
=
V
˙
t
=
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
.Since T
∗
(
t
)
,
C
∗
(
t
)
≥
0, then T
∗
(
t
)
=
C
∗
(
t
)
=
0. Hence, d
U
0
/
d
t
=
0 if and only if T
(
t
)
=
T
0, T
∗
(
t
)
=
0, C
∗
(
t
)
=
0, V
(
t
)
=
0, and Z
(
t
)
=
0. It follows from LaSalle’s invariance principle (LIP) that S
0 is GAS when R
0
≤
1.Theorem 4.
For system (3)–(7), assume that R
1
≤
1
<
R
0; then, S
1 is GAS.Proof.
Define(23)
U
1
=
T
1
H
T
T
1
+
η
1
T
1
∗
H
T
∗
T
1
∗
+
η
2
C
1
∗
H
C
∗
C
1
∗
+
η
3
V
1
H
V
V
1
+
η
4
Z
+
η
5
T
1
V
1
∫
0
τ
1
H
T
t
-
θ
V
t
-
θ
T
1
V
1
d
θ
+
η
6
T
1
∗
∫
0
τ
2
H
T
∗
t
-
θ
T
1
∗
d
θ
+
η
7
C
1
∗
∫
0
τ
3
H
C
∗
t
-
θ
C
1
∗
d
θ
.Then, d
U
1
/
d
t is given by(24)
d
U
1
d
t
=
1
-
T
1
T
π
-
1
-
ε
k
T
V
-
d
T
+
η
1
1
-
T
1
∗
T
∗
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
b
T
∗
+
η
2
1
-
C
1
∗
C
∗
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
a
C
∗
+
η
3
1
-
V
1
V
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
g
V
Z
-
μ
Z
+
η
5
T
V
-
T
t
-
τ
1
V
t
-
τ
1
+
T
1
V
1
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
6
T
∗
-
T
∗
t
-
τ
2
+
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
7
C
∗
-
C
∗
t
-
τ
3
+
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Equation (24) can be simplified as(25)
d
U
1
d
t
=
1
-
T
1
T
π
-
d
T
-
e
m
1
τ
1
η
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
T
∗
+
η
1
b
T
1
∗
-
e
-
m
1
τ
1
η
2
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
C
∗
+
η
2
a
C
1
∗
-
b
η
1
T
∗
t
-
τ
2
V
1
V
-
a
η
2
C
∗
t
-
τ
3
V
1
V
+
c
η
3
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
+
η
5
T
1
V
1
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Applying π
=
d
T
1
+
(
1
-
ε
)
k
T
1
V
1, we get(26)
d
U
1
d
t
=
1
-
T
1
T
d
T
1
-
d
T
+
1
-
ε
k
T
1
V
1
1
-
T
1
T
-
e
-
m
1
τ
1
η
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
T
∗
+
η
1
b
T
1
∗
-
e
-
m
1
τ
1
η
2
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
C
∗
+
η
2
a
C
1
∗
-
b
η
1
T
∗
t
-
τ
2
V
1
V
-
a
η
2
C
∗
t
-
τ
3
V
1
V
+
c
η
3
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
+
T
1
V
1
e
-
m
1
τ
1
η
1
1
-
ξ
1
-
ε
k
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
T
1
V
1
e
-
m
1
τ
1
η
2
ξ
1
-
ε
k
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the steady state conditions for S
1,(27)
1
-
ξ
1
-
ε
k
T
1
V
1
=
e
m
1
τ
1
b
T
1
∗
,
ξ
1
-
ε
k
T
1
V
1
=
e
m
1
τ
1
a
C
1
∗
,
c
V
1
=
e
-
m
2
τ
2
N
T
b
T
1
∗
+
e
-
m
3
τ
3
N
C
a
C
1
∗
,we get (
1
-
ε
)
k
T
1
V
1
=
η
1
b
T
1
∗
+
η
2
a
C
1
∗ and(28)
d
U
1
d
t
=
-
d
T
-
T
1
2
T
+
η
1
b
T
1
∗
1
-
T
1
T
+
η
2
a
C
1
∗
1
-
T
1
T
-
η
1
b
T
1
∗
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
T
1
V
1
T
∗
+
η
1
b
T
1
∗
-
η
2
a
C
1
∗
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
T
1
V
1
C
∗
+
η
2
a
C
1
∗
-
η
1
b
T
1
∗
V
1
T
∗
t
-
τ
2
V
T
1
∗
-
η
2
a
C
1
∗
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
η
1
b
T
1
∗
+
η
2
a
C
1
∗
+
r
η
3
V
1
Z
-
μ
η
4
Z
+
η
1
b
T
1
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
2
a
C
1
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Consider the following equalities:(29)
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
=
ln
⁡
T
j
T
+
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
j
∗
T
j
V
j
T
∗
+
ln
⁡
V
j
T
∗
V
T
j
∗
,
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
=
ln
⁡
T
j
T
+
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
C
j
∗
T
j
V
j
C
∗
+
ln
⁡
V
j
C
∗
V
C
j
∗
,
ln
⁡
T
∗
t
-
τ
2
T
∗
=
ln
⁡
T
∗
t
-
τ
2
V
j
V
T
j
∗
+
ln
⁡
V
T
j
∗
V
j
T
∗
,
ln
⁡
C
∗
t
-
τ
3
C
∗
=
ln
⁡
C
∗
t
-
τ
3
V
j
V
C
j
∗
+
ln
⁡
V
C
j
∗
V
j
C
∗
,
j
=
1,2
.Using (29) in case of j
=
1, we obtain(30)
d
U
1
d
t
=
-
d
T
-
T
1
2
T
-
η
1
b
T
1
∗
T
1
T
-
1
-
l
n
T
1
T
-
η
2
a
C
1
∗
T
1
T
-
1
-
ln
⁡
T
1
T
-
η
1
b
T
1
∗
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
T
1
V
1
T
∗
-
1
-
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
T
1
V
1
T
∗
-
η
2
a
C
1
∗
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
T
1
V
1
C
∗
-
1
-
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
T
1
V
1
C
∗
-
η
1
b
T
1
∗
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
1
-
ln
⁡
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
η
2
a
C
1
∗
C
∗
t
-
τ
3
V
1
V
C
1
∗
-
1
-
ln
⁡
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
r
η
3
V
1
-
μ
g
Z
=
-
d
T
-
T
1
2
T
-
η
1
b
T
1
∗
H
T
1
T
+
H
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
T
1
V
1
T
∗
+
H
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
η
2
a
C
1
∗
H
T
1
T
+
H
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
T
1
V
1
C
∗
+
H
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
r
η
3
d
g
+
μ
1
-
ε
k
g
1
-
ε
k
R
1
-
1
Z
.We have if R
0
>
1, then T
1
,
T
1
∗
,
C
1
∗
,
V
1
>
0. Since R
1
≤
1, then for all T
,
T
∗
,
C
∗
,
V
,
Z
>
0, we have d
U
1
/
d
t
≤
0. We note that d
U
1
/
d
t
=
0 at S
1. Then, from LIP, S
1 is GAS.Theorem 5.
For system (3)–(7), assume that R
1
>
1; then, S
2 is GAS.Proof.
We consider(31)
U
2
=
T
2
H
T
T
2
+
η
1
T
2
∗
H
T
∗
T
2
∗
+
η
2
C
2
∗
H
C
∗
C
2
∗
+
η
3
V
2
H
V
V
2
+
η
4
Z
2
H
Z
Z
2
+
η
5
T
2
V
2
∫
0
τ
1
H
T
t
-
θ
V
t
-
θ
T
2
V
2
d
θ
+
η
6
T
2
∗
∫
0
τ
2
H
T
∗
t
-
θ
T
2
∗
d
θ
+
η
7
C
2
∗
∫
0
τ
3
H
C
∗
t
-
θ
C
2
∗
d
θ
.Function U
2 along the trajectories of system (3)–(7) satisfies(32)
d
U
2
d
t
=
1
-
T
2
T
π
-
1
-
ε
k
T
V
-
d
T
+
η
1
1
-
T
2
∗
T
∗
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
b
T
∗
+
η
2
1
-
C
2
∗
C
∗
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
a
C
∗
+
η
3
1
-
V
2
V
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
1
-
Z
2
Z
g
V
Z
-
μ
Z
+
η
5
T
V
-
T
t
-
τ
1
V
t
-
τ
1
+
T
2
V
2
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
6
T
∗
-
T
∗
t
-
τ
2
+
T
2
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
7
C
∗
-
C
∗
t
-
τ
3
+
C
2
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the steady state conditions for S
2,(33)
π
=
d
T
2
+
1
-
ε
k
T
2
V
2
,
1
-
ξ
1
-
ε
k
T
2
V
2
=
e
m
1
τ
1
b
T
2
∗
,
ξ
1
-
ε
k
T
2
V
2
=
e
m
1
τ
1
a
C
2
∗
,
c
V
2
=
e
-
m
2
τ
2
N
T
b
T
2
∗
+
e
-
m
3
τ
3
N
C
a
C
2
∗
-
r
V
2
Z
2
,we get (
1
-
ε
)
k
T
2
V
2
=
η
1
b
T
2
∗
+
η
2
a
C
2
∗ and(34)
d
U
2
d
t
=
-
d
T
-
T
2
2
T
+
η
1
b
T
2
∗
1
-
T
2
T
+
η
2
a
C
2
∗
1
-
T
2
T
-
η
1
b
T
2
∗
T
t
-
τ
1
V
t
-
τ
1
T
2
∗
T
2
V
2
T
∗
+
η
1
b
T
2
∗
-
η
2
a
C
2
∗
T
t
-
τ
1
V
t
-
τ
1
C
2
∗
T
2
V
2
C
∗
+
η
2
a
C
2
∗
-
η
1
b
T
2
∗
V
2
T
∗
t
-
τ
2
V
T
2
∗
-
η
2
a
C
2
∗
V
2
C
∗
t
-
τ
3
V
C
2
∗
+
η
1
b
T
2
∗
+
η
2
a
C
2
∗
+
η
1
b
T
2
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
2
a
C
2
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
1
b
T
2
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
2
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using (29) in case of j
=
2, we obtain(35)
d
U
2
d
t
=
-
d
T
-
T
2
2
T
-
η
1
b
T
2
∗
H
T
2
T
+
H
T
t
-
τ
1
V
t
-
τ
1
T
2
∗
T
2
V
2
T
∗
+
H
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
η
2
a
C
2
∗
H
T
2
T
+
H
T
t
-
τ
1
V
t
-
τ
1
C
2
∗
T
2
V
2
C
∗
+
H
C
∗
t
-
τ
3
V
2
V
C
2
∗
.Since R
1
>
1, then T
2
,
T
2
∗
,
C
2
∗
,
V
2, and Z
2
>
0. It is observed that d
U
2
/
d
t
=
0 if and only if T
=
T
2, T
∗
=
T
2
∗, C
∗
=
C
2
∗, and V
=
V
2. Therefore, if V
=
V
2, then V
˙
=
0 and (6) becomes(36)
0
=
e
-
m
2
τ
2
N
T
b
T
2
∗
+
e
-
m
3
τ
3
N
C
a
C
2
∗
-
c
V
2
-
r
V
2
Z
,which gives Z
=
Z
2. Hence, d
U
2
/
d
t is equal to zero at S
2. The global stability of S
2 follows from LIP.


## 2.1. Positive Invariance

Proposition 1.
The solution(
T
(
t
)
,
T
∗
(
t
)
,
C
∗
(
t
)
,
V
(
t
)
,
Z
(
t
)
) of (3)–(7) with initial conditions (8) is nonnegative for t
⩾
0 and ultimately bounded.Proof.
Assume thatT
(
t
)
≤
0 on [
0
,
ρ
] for some constant ρ and T
(
t
∗
)
=
0, t
∗
∈
[
0
,
ρ
]. From (3), we get T
˙
(
t
∗
)
=
π
>
0 and hence T
(
t
)
>
0, for some t
∈
(
t
∗
,
t
∗
+
ϱ
), and sufficiently small ϱ
>
0. This leads to a contradiction; therefore, T
(
t
)
>
0, for all t
⩾
0. From (4), (5), and (6), we have(9)
T
∗
t
=
φ
2
0
e
-
b
t
+
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
∫
0
t
e
-
b
t
-
η
T
η
-
τ
1
V
η
-
τ
1
d
η
,
C
∗
t
=
φ
3
0
e
-
a
t
+
e
-
m
1
τ
1
ξ
1
-
ε
k
∫
0
t
e
-
a
t
-
η
T
η
-
τ
1
V
η
-
τ
1
d
η
,
V
t
=
φ
4
0
e
-
∫
0
t
c
+
r
Z
ς
d
ς
+
e
-
m
2
τ
2
N
T
b
∫
0
t
e
-
∫
η
t
c
+
r
Z
ς
d
ς
T
∗
η
-
τ
2
d
η
+
e
-
m
3
τ
3
N
C
a
∫
0
t
e
-
∫
η
t
c
+
r
Z
ς
d
ς
C
∗
η
-
τ
3
d
η
.This confirms that T
∗
(
t
)
⩾
0, C
∗
(
t
)
⩾
0, and V
(
t
)
⩾
0, for all t
∈
[
0
,
τ
]. By a recursive argument, we get that T
∗
(
t
)
⩾
0, C
∗
(
t
)
⩾
0, and V
(
t
)
⩾
0, for all t
⩾
0. Moreover, from (7), we obtain(10)
Z
t
=
φ
5
0
e
-
∫
0
t
μ
-
g
V
ς
d
ς
. Clearly, Z
(
t
)
⩾
0, t
⩾
0. Now, we let G
1
(
t
)
=
e
-
m
1
τ
1
T
(
t
-
τ
1
)
+
T
∗
(
t
)
+
C
∗
(
t
), G
2
(
t
)
=
V
(
t
)
+
r
/
g
Z
(
t
), s
1
=
min
⁡
d
,
a
,
b, s
2
=
min
⁡
c
,
μ, L
1
=
π
/
s
1, L
2
=
(
N
T
b
+
N
C
a
)
L
1
/
s
2, and L
3
=
g
L
2
/
r; then,(11)
G
˙
1
t
=
e
-
m
1
τ
1
π
-
d
T
t
-
τ
1
-
b
T
∗
t
-
a
C
∗
t
≤
π
e
-
m
1
τ
1
-
s
1
e
-
m
1
τ
1
T
t
-
τ
1
+
T
∗
t
+
C
∗
t
=
π
e
-
m
1
τ
1
-
s
1
G
1
t
<
π
-
s
1
G
1
t
.Hence, l
i
m
s
u
p
t
→
∞
G
1
(
t
)
≤
L
1. Since T
(
t
)
>
0, T
∗
(
t
)
≥
0, and C
∗
(
t
)
≥
0, then l
i
m
s
u
p
t
→
∞
T
∗
(
t
)
≤
L
1 and l
i
m
s
u
p
t
→
∞
C
∗
(
t
)
≤
L
1. Also,(12)
G
˙
2
t
=
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
c
V
t
-
r
μ
g
Z
t
≤
e
-
m
2
τ
2
N
T
b
+
e
-
m
3
τ
3
N
C
a
L
1
-
s
2
V
t
+
r
g
Z
t
<
N
T
b
+
N
C
a
L
1
-
s
2
G
2
t
.Hence, l
i
m
s
u
p
t
→
∞
G
2
(
t
)
≤
L
2. Since V
(
t
) and Z
(
t
) are nonnegative, then l
i
m
s
u
p
t
→
∞
V
(
t
)
≤
L
2 and l
i
m
s
u
p
t
→
∞
Z
(
t
)
≤
L
3. Therefore, T
(
t
)
,
T
∗
(
t
)
,
C
∗
(
t
)
,
V
(
t
), and Z
(
t
) are ultimately bounded.


## 2.2. Steady States

System (3)–(7) always admits an uninfected steady state S
0
=
(
π
/
d
,
0,0
,
0,0
). Let T
0
=
π
/
d
. Now, we define the basic reproduction number for system (3)–(7) as(13)
R
0
=
1
-
ε
k
T
0
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
c
. The system has two other steady states, infected steady state without humoral immune response S
1
=
(
T
1
,
T
1
∗
,
C
1
∗
,
V
1
,
0
) and infected steady state with humoral immune response S
2
=
(
T
2
,
T
2
∗
,
C
2
∗
,
V
2
,
Z
2
):(14)
T
1
=
T
0
R
0
,
T
1
∗
=
e
-
m
1
τ
1
1
-
ξ
π
b
R
0
-
1
R
0
,
C
1
∗
=
e
-
m
1
τ
1
ξ
π
a
R
0
-
1
R
0
,
V
1
=
d
1
-
ε
k
R
0
-
1
,
T
2
=
π
g
g
d
+
1
-
ε
k
μ
,
T
2
∗
=
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
π
μ
b
d
g
+
1
-
ε
k
μ
,
C
2
∗
=
e
-
m
1
τ
1
ξ
1
-
ε
k
π
μ
a
d
g
+
1
-
ε
k
μ
,
V
2
=
μ
g
,
Z
2
=
c
r
d
g
R
0
d
g
+
1
-
ε
k
μ
-
1
.We note that T
1, T
1
∗, C
1
∗, and V
1 are positive when R
0
>
1 and where T
2, T
2
∗, C
2
∗, and V
2
>
0 and Z
2
>
0 when d
g
R
0
/
d
g
+
(
1
-
ε
)
k
μ
>
1. Now, we define humoral immune response reproduction number R
1 as(15)
R
1
=
R
0
1
+
1
-
ε
k
μ
/
d
g
.Clearly, R
0
>
R
1. From above, we can state the following lemma.Lemma 2.
For system (3)–(7), one has the following:(i)
IfR
0
≤
1, then the system has only one positive steady state S
0.
(ii)
IfR
1
≤
1
<
R
0, then the system has two positive steady states S
0 and S
1.
(iii)
IfR
1
>
1, then the system has three positive steady states S
0, S
1, and S
2.


## 2.3. Global Stability Analysis

We establish the global stability of all the steady states of system (3)–(7) employing the method of Lyapunov function. Let us define(16)
H
θ
=
θ
-
1
-
ln
⁡
θ
.Theorem 3.
For system (3)–(7), if R
0
≤
1, then S
0 is GAS.Proof.
Define(17)
U
0
=
T
0
H
T
T
0
+
η
1
T
∗
+
η
2
C
∗
+
η
3
V
+
η
4
Z
+
η
5
∫
0
τ
1
T
t
-
θ
V
t
-
θ
d
θ
+
η
6
∫
0
τ
2
T
∗
t
-
θ
d
θ
+
η
7
∫
0
τ
3
C
∗
t
-
θ
d
θ
,where η
i
,
i
=
1
,
…
,
7, are positive constants satisfying the following equations:(18)
e
-
m
1
τ
1
1
-
ε
k
1
-
ξ
η
1
+
ξ
η
2
=
η
5
,
η
6
-
b
η
1
=
0
,
η
7
-
a
η
2
=
0
,
η
6
-
η
3
e
-
m
2
τ
2
N
T
b
=
0
,
η
7
-
η
3
e
-
m
3
τ
3
N
C
a
=
0
,
g
η
4
-
r
η
3
=
0
,
η
5
=
1
-
ε
k
.The solution of (18) is given by(19)
η
1
=
N
T
e
m
1
τ
1
+
m
3
τ
3
e
m
3
τ
3
1
-
ξ
N
T
+
e
m
2
τ
2
ξ
N
C
,
η
2
=
N
C
e
m
1
τ
1
+
m
2
τ
2
e
m
3
τ
3
1
-
ξ
N
T
+
e
m
2
τ
2
ξ
N
C
,
η
3
=
e
m
1
τ
1
+
m
2
τ
2
+
m
3
τ
3
e
m
3
τ
3
1
-
ξ
N
T
+
e
m
2
τ
2
ξ
N
C
,
η
4
=
r
e
m
1
τ
1
+
m
2
τ
2
+
m
3
τ
3
g
e
m
3
τ
3
1
-
ξ
N
T
+
e
m
2
τ
2
ξ
N
C
,
η
5
=
1
-
ε
k
,
η
6
=
b
N
T
e
m
1
τ
1
+
m
3
τ
3
e
m
3
τ
3
1
-
ξ
N
T
+
e
m
2
τ
2
ξ
N
C
,
η
7
=
a
N
C
e
m
1
τ
1
+
m
2
τ
2
e
m
3
τ
3
1
-
ξ
N
T
+
e
m
2
τ
2
ξ
N
C
.The values of η
i
,
i
=
1
,
…
,
7, will be used through the paper. Calculating the derivative of U
0 along the solutions of system (3)–(7) and applying π
=
T
0
d, we obtain(20)
d
U
0
d
t
=
1
-
T
0
T
π
-
1
-
ε
k
T
V
-
d
T
+
η
1
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
b
T
∗
+
η
2
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
a
C
∗
+
η
3
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
g
V
Z
-
μ
Z
+
η
5
T
V
-
η
5
T
t
-
τ
1
V
t
-
τ
1
+
η
6
T
∗
-
η
6
T
∗
t
-
τ
2
+
η
7
C
∗
-
η
7
C
∗
t
-
τ
3
.Then,(21)
d
U
0
d
t
=
1
-
T
0
T
π
-
d
T
+
1
-
ε
k
T
0
V
-
η
3
c
V
-
η
4
μ
Z
=
-
d
T
-
T
0
2
T
+
c
η
3
R
0
-
1
V
-
η
4
μ
Z
.Therefore, if R
0
≤
1, then d
U
0
/
d
t
≤
0, for all T
,
V
,
Z
>
0. The solutions of system (3)–(7) limit M, the largest invariant subset of d
U
0
/
d
t
=
0 [22]. We note that d
U
0
/
d
t
=
0 if and only if T
(
t
)
=
T
0, V
(
t
)
=
0, and Z
(
t
)
=
0. For each element of M, we have V
(
t
)
=
0 and Z
(
t
)
=
0,; then, V
˙
(
t
)
=
0 and(22)
0
=
V
˙
t
=
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
.Since T
∗
(
t
)
,
C
∗
(
t
)
≥
0, then T
∗
(
t
)
=
C
∗
(
t
)
=
0. Hence, d
U
0
/
d
t
=
0 if and only if T
(
t
)
=
T
0, T
∗
(
t
)
=
0, C
∗
(
t
)
=
0, V
(
t
)
=
0, and Z
(
t
)
=
0. It follows from LaSalle’s invariance principle (LIP) that S
0 is GAS when R
0
≤
1.Theorem 4.
For system (3)–(7), assume that R
1
≤
1
<
R
0; then, S
1 is GAS.Proof.
Define(23)
U
1
=
T
1
H
T
T
1
+
η
1
T
1
∗
H
T
∗
T
1
∗
+
η
2
C
1
∗
H
C
∗
C
1
∗
+
η
3
V
1
H
V
V
1
+
η
4
Z
+
η
5
T
1
V
1
∫
0
τ
1
H
T
t
-
θ
V
t
-
θ
T
1
V
1
d
θ
+
η
6
T
1
∗
∫
0
τ
2
H
T
∗
t
-
θ
T
1
∗
d
θ
+
η
7
C
1
∗
∫
0
τ
3
H
C
∗
t
-
θ
C
1
∗
d
θ
.Then, d
U
1
/
d
t is given by(24)
d
U
1
d
t
=
1
-
T
1
T
π
-
1
-
ε
k
T
V
-
d
T
+
η
1
1
-
T
1
∗
T
∗
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
b
T
∗
+
η
2
1
-
C
1
∗
C
∗
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
a
C
∗
+
η
3
1
-
V
1
V
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
g
V
Z
-
μ
Z
+
η
5
T
V
-
T
t
-
τ
1
V
t
-
τ
1
+
T
1
V
1
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
6
T
∗
-
T
∗
t
-
τ
2
+
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
7
C
∗
-
C
∗
t
-
τ
3
+
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Equation (24) can be simplified as(25)
d
U
1
d
t
=
1
-
T
1
T
π
-
d
T
-
e
m
1
τ
1
η
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
T
∗
+
η
1
b
T
1
∗
-
e
-
m
1
τ
1
η
2
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
C
∗
+
η
2
a
C
1
∗
-
b
η
1
T
∗
t
-
τ
2
V
1
V
-
a
η
2
C
∗
t
-
τ
3
V
1
V
+
c
η
3
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
+
η
5
T
1
V
1
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Applying π
=
d
T
1
+
(
1
-
ε
)
k
T
1
V
1, we get(26)
d
U
1
d
t
=
1
-
T
1
T
d
T
1
-
d
T
+
1
-
ε
k
T
1
V
1
1
-
T
1
T
-
e
-
m
1
τ
1
η
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
T
∗
+
η
1
b
T
1
∗
-
e
-
m
1
τ
1
η
2
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
C
∗
+
η
2
a
C
1
∗
-
b
η
1
T
∗
t
-
τ
2
V
1
V
-
a
η
2
C
∗
t
-
τ
3
V
1
V
+
c
η
3
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
+
T
1
V
1
e
-
m
1
τ
1
η
1
1
-
ξ
1
-
ε
k
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
T
1
V
1
e
-
m
1
τ
1
η
2
ξ
1
-
ε
k
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the steady state conditions for S
1,(27)
1
-
ξ
1
-
ε
k
T
1
V
1
=
e
m
1
τ
1
b
T
1
∗
,
ξ
1
-
ε
k
T
1
V
1
=
e
m
1
τ
1
a
C
1
∗
,
c
V
1
=
e
-
m
2
τ
2
N
T
b
T
1
∗
+
e
-
m
3
τ
3
N
C
a
C
1
∗
,we get (
1
-
ε
)
k
T
1
V
1
=
η
1
b
T
1
∗
+
η
2
a
C
1
∗ and(28)
d
U
1
d
t
=
-
d
T
-
T
1
2
T
+
η
1
b
T
1
∗
1
-
T
1
T
+
η
2
a
C
1
∗
1
-
T
1
T
-
η
1
b
T
1
∗
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
T
1
V
1
T
∗
+
η
1
b
T
1
∗
-
η
2
a
C
1
∗
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
T
1
V
1
C
∗
+
η
2
a
C
1
∗
-
η
1
b
T
1
∗
V
1
T
∗
t
-
τ
2
V
T
1
∗
-
η
2
a
C
1
∗
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
η
1
b
T
1
∗
+
η
2
a
C
1
∗
+
r
η
3
V
1
Z
-
μ
η
4
Z
+
η
1
b
T
1
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
2
a
C
1
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Consider the following equalities:(29)
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
=
ln
⁡
T
j
T
+
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
j
∗
T
j
V
j
T
∗
+
ln
⁡
V
j
T
∗
V
T
j
∗
,
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
=
ln
⁡
T
j
T
+
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
C
j
∗
T
j
V
j
C
∗
+
ln
⁡
V
j
C
∗
V
C
j
∗
,
ln
⁡
T
∗
t
-
τ
2
T
∗
=
ln
⁡
T
∗
t
-
τ
2
V
j
V
T
j
∗
+
ln
⁡
V
T
j
∗
V
j
T
∗
,
ln
⁡
C
∗
t
-
τ
3
C
∗
=
ln
⁡
C
∗
t
-
τ
3
V
j
V
C
j
∗
+
ln
⁡
V
C
j
∗
V
j
C
∗
,
j
=
1,2
.Using (29) in case of j
=
1, we obtain(30)
d
U
1
d
t
=
-
d
T
-
T
1
2
T
-
η
1
b
T
1
∗
T
1
T
-
1
-
l
n
T
1
T
-
η
2
a
C
1
∗
T
1
T
-
1
-
ln
⁡
T
1
T
-
η
1
b
T
1
∗
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
T
1
V
1
T
∗
-
1
-
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
T
1
V
1
T
∗
-
η
2
a
C
1
∗
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
T
1
V
1
C
∗
-
1
-
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
T
1
V
1
C
∗
-
η
1
b
T
1
∗
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
1
-
ln
⁡
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
η
2
a
C
1
∗
C
∗
t
-
τ
3
V
1
V
C
1
∗
-
1
-
ln
⁡
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
r
η
3
V
1
-
μ
g
Z
=
-
d
T
-
T
1
2
T
-
η
1
b
T
1
∗
H
T
1
T
+
H
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
T
1
V
1
T
∗
+
H
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
η
2
a
C
1
∗
H
T
1
T
+
H
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
T
1
V
1
C
∗
+
H
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
r
η
3
d
g
+
μ
1
-
ε
k
g
1
-
ε
k
R
1
-
1
Z
.We have if R
0
>
1, then T
1
,
T
1
∗
,
C
1
∗
,
V
1
>
0. Since R
1
≤
1, then for all T
,
T
∗
,
C
∗
,
V
,
Z
>
0, we have d
U
1
/
d
t
≤
0. We note that d
U
1
/
d
t
=
0 at S
1. Then, from LIP, S
1 is GAS.Theorem 5.
For system (3)–(7), assume that R
1
>
1; then, S
2 is GAS.Proof.
We consider(31)
U
2
=
T
2
H
T
T
2
+
η
1
T
2
∗
H
T
∗
T
2
∗
+
η
2
C
2
∗
H
C
∗
C
2
∗
+
η
3
V
2
H
V
V
2
+
η
4
Z
2
H
Z
Z
2
+
η
5
T
2
V
2
∫
0
τ
1
H
T
t
-
θ
V
t
-
θ
T
2
V
2
d
θ
+
η
6
T
2
∗
∫
0
τ
2
H
T
∗
t
-
θ
T
2
∗
d
θ
+
η
7
C
2
∗
∫
0
τ
3
H
C
∗
t
-
θ
C
2
∗
d
θ
.Function U
2 along the trajectories of system (3)–(7) satisfies(32)
d
U
2
d
t
=
1
-
T
2
T
π
-
1
-
ε
k
T
V
-
d
T
+
η
1
1
-
T
2
∗
T
∗
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
b
T
∗
+
η
2
1
-
C
2
∗
C
∗
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
-
a
C
∗
+
η
3
1
-
V
2
V
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
1
-
Z
2
Z
g
V
Z
-
μ
Z
+
η
5
T
V
-
T
t
-
τ
1
V
t
-
τ
1
+
T
2
V
2
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
6
T
∗
-
T
∗
t
-
τ
2
+
T
2
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
7
C
∗
-
C
∗
t
-
τ
3
+
C
2
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the steady state conditions for S
2,(33)
π
=
d
T
2
+
1
-
ε
k
T
2
V
2
,
1
-
ξ
1
-
ε
k
T
2
V
2
=
e
m
1
τ
1
b
T
2
∗
,
ξ
1
-
ε
k
T
2
V
2
=
e
m
1
τ
1
a
C
2
∗
,
c
V
2
=
e
-
m
2
τ
2
N
T
b
T
2
∗
+
e
-
m
3
τ
3
N
C
a
C
2
∗
-
r
V
2
Z
2
,we get (
1
-
ε
)
k
T
2
V
2
=
η
1
b
T
2
∗
+
η
2
a
C
2
∗ and(34)
d
U
2
d
t
=
-
d
T
-
T
2
2
T
+
η
1
b
T
2
∗
1
-
T
2
T
+
η
2
a
C
2
∗
1
-
T
2
T
-
η
1
b
T
2
∗
T
t
-
τ
1
V
t
-
τ
1
T
2
∗
T
2
V
2
T
∗
+
η
1
b
T
2
∗
-
η
2
a
C
2
∗
T
t
-
τ
1
V
t
-
τ
1
C
2
∗
T
2
V
2
C
∗
+
η
2
a
C
2
∗
-
η
1
b
T
2
∗
V
2
T
∗
t
-
τ
2
V
T
2
∗
-
η
2
a
C
2
∗
V
2
C
∗
t
-
τ
3
V
C
2
∗
+
η
1
b
T
2
∗
+
η
2
a
C
2
∗
+
η
1
b
T
2
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
2
a
C
2
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
T
V
+
η
1
b
T
2
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
2
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using (29) in case of j
=
2, we obtain(35)
d
U
2
d
t
=
-
d
T
-
T
2
2
T
-
η
1
b
T
2
∗
H
T
2
T
+
H
T
t
-
τ
1
V
t
-
τ
1
T
2
∗
T
2
V
2
T
∗
+
H
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
η
2
a
C
2
∗
H
T
2
T
+
H
T
t
-
τ
1
V
t
-
τ
1
C
2
∗
T
2
V
2
C
∗
+
H
C
∗
t
-
τ
3
V
2
V
C
2
∗
.Since R
1
>
1, then T
2
,
T
2
∗
,
C
2
∗
,
V
2, and Z
2
>
0. It is observed that d
U
2
/
d
t
=
0 if and only if T
=
T
2, T
∗
=
T
2
∗, C
∗
=
C
2
∗, and V
=
V
2. Therefore, if V
=
V
2, then V
˙
=
0 and (6) becomes(36)
0
=
e
-
m
2
τ
2
N
T
b
T
2
∗
+
e
-
m
3
τ
3
N
C
a
C
2
∗
-
c
V
2
-
r
V
2
Z
,which gives Z
=
Z
2. Hence, d
U
2
/
d
t is equal to zero at S
2. The global stability of S
2 follows from LIP.


## 3. HIV Dynamics Model with Saturated Incidence

We present an HIV infection model with saturated incidence:(37)
T
˙
t
=
π
-
1
-
ε
k
T
t
V
t
1
+
β
V
t
-
d
T
t
,
T
˙
∗
t
=
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
b
T
∗
t
,
C
˙
∗
t
=
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
a
C
∗
t
,
V
˙
t
=
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
t
Z
t
-
c
V
t
,
Z
˙
t
=
g
V
t
Z
t
-
μ
Z
t
,where β
>
0 is the saturation incidence rate constant. Similar to the previous section, one can show that the solutions of the model are nonnegative and bounded.


### 3.1. Steady States

System (37) always admits an uninfected steady state S
0
=
(
T
0
,
0,0
,
0,0
), where T
0
=
π
/
d. Now, we define the basic reproduction number for system (37) as(38)
R
0
=
1
-
ε
k
T
0
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
c
. The system has two other steady states S
1
(
T
1
,
T
1
∗
,
C
1
∗
,
V
1
,
0
) and S
2
(
T
2
,
T
2
∗
,
C
2
∗
,
V
2
,
Z
2
), where(39)
T
1
=
β
π
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
+
c
1
-
ε
k
+
d
β
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
,
T
1
∗
=
e
-
m
1
τ
1
1
-
ξ
c
d
b
1
-
ε
k
+
d
β
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
R
0
-
1
,
C
1
∗
=
e
-
m
1
τ
1
ξ
c
d
a
1
-
ε
k
+
d
β
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
R
0
-
1
,
V
1
=
d
1
-
ε
k
+
d
β
R
0
-
1
,
T
2
=
π
g
+
β
μ
g
d
+
1
-
ε
k
μ
+
d
β
μ
,
T
2
∗
=
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
π
μ
b
g
d
+
1
-
ε
k
μ
+
d
β
μ
,
C
2
∗
=
e
-
m
1
τ
1
ξ
1
-
ε
k
π
μ
a
g
d
+
1
-
ε
k
μ
+
d
β
μ
,
V
2
=
μ
g
,
Z
2
=
c
r
d
g
R
0
d
g
+
1
-
ε
k
μ
+
d
β
μ
-
1
.We note that T
1, T
1
∗, C
1
∗, and V
1 are positive when R
0
>
1. And T
2
,
T
2
∗
,
C
2
∗
,
V
2
>
0 and Z
2
>
0 when d
g
R
0
/
d
g
+
(
1
-
ε
)
k
μ
>
1. Now, we define another threshold parameter R
1 as(40)
R
1
=
R
0
1
+
1
-
ε
k
μ
+
d
β
μ
/
d
g
.Clearly, R
0
>
R
1. From above, we have the following result.Lemma 6.
For system (37), one has the following:(i)
IfR
0
≤
1, then the system has only one positive steady state S
0.
(ii)
IfR
1
≤
1
<
R
0, then the system has two positive steady states S
0 and S
1.
(iii)
IfR
1
>
1, then the system has three positive steady states S
0, S
1, and S
2.


### 3.2. Global Stability Analysis

In this subsection, we investigate the global stability of system (37) by constructing suitable Lyapunov functionals and applying LaSalle invariance principle.Theorem 7.
For system (37), if R
0
≤
1, then S
0 is GAS.Proof.
Define(41)
U
0
=
T
0
H
T
T
0
+
η
1
T
∗
+
η
2
C
∗
+
η
3
V
+
η
4
Z
+
η
5
∫
0
τ
1
T
t
-
θ
V
t
-
θ
1
+
β
V
t
-
θ
d
θ
+
η
6
∫
0
τ
2
T
∗
t
-
θ
d
θ
+
η
7
∫
0
τ
3
C
∗
t
-
θ
d
θ
.Calculating the derivative of U
0 along the solutions of system (37), we obtain(42)
d
U
0
d
t
=
1
-
T
0
T
π
-
1
-
ε
k
T
V
1
+
β
V
-
d
T
+
η
1
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
b
T
∗
+
η
2
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
a
C
∗
+
η
3
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
g
V
Z
-
μ
Z
+
η
5
T
V
1
+
β
V
-
η
5
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
+
η
6
T
∗
-
η
6
T
∗
t
-
τ
2
+
η
7
C
∗
-
η
7
C
∗
t
-
τ
3
.Collecting terms of (42), we obtain(43)
d
U
0
d
t
=
1
-
T
0
T
π
-
d
T
+
1
-
ε
k
T
0
V
1
+
β
V
-
η
3
c
V
-
η
4
μ
Z
=
-
d
T
-
T
0
2
T
+
η
3
c
β
R
0
V
2
1
+
β
V
+
η
4
μ
Z
+
c
η
3
R
0
-
1
V
.It follows that S
0 is GAS when R
0
≤
1.Theorem 8.
For system (37), assume that R
1
≤
1
<
R
0; then, S
1 is GAS.Proof.
Consider(44)
U
1
=
T
1
H
T
T
1
+
η
1
T
1
∗
H
T
∗
T
1
∗
+
η
2
C
1
∗
H
C
∗
C
1
∗
+
η
3
V
1
H
V
V
1
+
η
4
Z
+
η
5
T
1
V
1
1
+
β
V
1
∫
0
τ
1
H
T
t
-
θ
V
t
-
θ
1
+
β
V
1
T
1
V
1
1
+
β
V
t
-
θ
d
θ
+
η
6
T
1
∗
∫
0
τ
2
H
T
∗
t
-
θ
T
1
∗
d
θ
+
η
7
C
1
∗
∫
0
τ
3
H
C
∗
t
-
θ
C
1
∗
d
θ
.Then, d
U
1
/
d
t is given by(45)
d
U
1
d
t
=
1
-
T
1
T
π
-
d
T
-
1
-
ε
k
T
V
1
+
β
V
+
η
1
1
-
T
1
∗
T
∗
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
b
T
∗
+
η
2
1
-
C
1
∗
C
∗
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
a
C
∗
+
η
3
1
-
V
1
V
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
g
V
Z
-
μ
Z
+
η
5
T
V
1
+
β
V
-
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
+
T
1
V
1
1
+
β
V
1
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
6
T
∗
-
T
∗
t
-
τ
2
+
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
7
C
∗
-
C
∗
t
-
τ
3
+
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Simplifying (45), we get(46)
d
U
1
d
t
=
1
-
T
1
T
π
-
d
T
+
1
-
ε
k
T
1
V
1
+
β
V
-
e
m
1
τ
1
η
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
T
1
∗
T
∗
+
η
1
b
T
1
∗
-
e
-
m
1
τ
1
η
2
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
C
1
∗
C
∗
+
η
2
a
C
1
∗
-
b
η
1
T
∗
t
-
τ
2
V
1
V
-
a
η
2
C
∗
t
-
τ
3
V
1
V
-
c
η
3
V
+
c
η
3
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
+
η
5
T
1
V
1
1
+
β
V
1
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Applying π
=
d
T
1
+
(
1
-
ε
)
k
T
1
V
1
/
1
+
β
V
1, we get(47)
d
U
1
d
t
=
1
-
T
1
T
d
T
1
-
d
T
+
1
-
ε
k
T
1
V
1
1
+
β
V
1
1
-
T
1
T
+
1
-
ε
k
T
1
V
1
+
β
V
-
e
-
m
1
τ
1
η
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
T
1
∗
T
∗
+
η
1
b
T
1
∗
-
e
-
m
1
τ
1
η
2
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
C
1
∗
C
∗
+
η
2
a
C
1
∗
-
b
η
1
T
∗
t
-
τ
2
V
1
V
-
a
η
2
C
∗
t
-
τ
3
V
1
V
-
c
η
3
V
+
c
η
3
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
+
T
1
V
1
e
-
m
1
τ
1
η
1
1
-
ξ
1
-
ε
k
1
+
β
V
1
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
T
1
V
1
e
-
m
1
τ
1
η
2
ξ
1
-
ε
k
1
+
β
V
1
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the steady state conditions for S
1,(48)
1
-
ξ
1
-
ε
k
T
1
V
1
1
+
β
V
1
=
e
m
1
τ
1
b
T
1
∗
,
ξ
1
-
ε
k
T
1
V
1
1
+
β
V
1
=
e
m
1
τ
1
a
C
1
∗
,
c
V
1
=
e
-
m
2
τ
2
N
T
b
T
1
∗
+
e
-
m
3
τ
3
N
C
a
C
1
∗
,we obtain (
1
-
ε
)
k
T
1
V
1
/
1
+
β
V
1
=
η
1
b
T
1
∗
+
η
2
a
C
1
∗ and(49)
d
U
1
d
t
=
-
d
T
-
T
1
2
T
+
1
-
ε
k
T
1
V
1
1
+
β
V
1
1
+
β
V
1
V
1
+
β
V
V
1
+
η
1
b
T
1
∗
1
-
T
1
T
+
η
2
a
C
1
∗
1
-
T
1
T
-
η
1
b
T
1
∗
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
1
+
β
V
1
T
1
V
1
T
∗
1
+
β
V
t
-
τ
1
+
η
1
b
T
1
∗
-
η
2
a
C
1
∗
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
1
+
β
V
1
T
1
V
1
C
∗
1
+
β
V
t
-
τ
1
+
η
2
a
C
1
∗
-
η
1
b
T
1
∗
V
1
T
∗
t
-
τ
2
V
T
1
∗
-
η
2
a
C
1
∗
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
η
1
b
T
1
∗
+
η
2
a
C
1
∗
-
1
-
ε
k
T
1
V
1
1
+
β
V
1
V
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
+
η
1
b
T
1
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
2
a
C
1
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Consider the following equalities:(50)
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
=
ln
⁡
T
j
T
+
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
T
j
∗
T
j
V
j
1
+
β
V
t
-
τ
1
T
∗
+
ln
⁡
V
j
T
∗
V
T
j
∗
+
ln
⁡
1
+
β
V
1
+
β
V
j
,
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
=
ln
⁡
T
j
T
+
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
j
C
j
∗
T
j
V
j
1
+
β
V
t
-
τ
1
C
∗
+
ln
⁡
V
j
C
∗
V
C
j
∗
+
ln
⁡
1
+
β
V
1
+
β
V
j
,
ln
⁡
T
∗
t
-
τ
2
T
∗
=
ln
⁡
T
∗
t
-
τ
2
V
j
V
T
j
∗
+
ln
⁡
V
T
j
∗
V
j
T
∗
,
ln
⁡
C
∗
t
-
τ
3
C
∗
=
ln
⁡
C
∗
t
-
τ
3
V
j
V
C
j
∗
+
ln
⁡
V
C
j
∗
V
j
C
∗
,
j
=
1,2
.Using (50), when j
=
1, we get(51)
d
U
1
d
t
=
-
d
T
-
T
1
2
T
+
1
-
ε
k
T
1
V
1
1
+
β
V
1
1
+
β
V
1
V
1
+
β
V
V
1
-
V
V
1
-
1
+
1
+
β
V
1
+
β
V
1
-
η
1
b
T
1
∗
T
1
T
-
1
-
ln
⁡
T
1
T
-
η
2
a
C
1
∗
T
1
T
-
1
-
ln
⁡
T
1
T
-
η
1
b
T
1
∗
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
T
1
∗
T
1
V
1
1
+
β
V
t
-
τ
1
T
∗
-
1
-
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
T
1
∗
T
1
V
1
1
+
β
V
t
-
τ
1
T
∗
-
η
2
a
C
1
∗
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
C
1
∗
T
1
V
1
1
+
β
V
t
-
τ
1
C
∗
-
1
-
l
n
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
C
1
∗
T
1
V
1
1
+
β
V
t
-
τ
1
C
∗
-
η
1
b
T
1
∗
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
1
-
ln
⁡
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
η
2
a
C
1
∗
C
∗
t
-
τ
3
V
1
V
C
1
∗
-
1
-
ln
⁡
C
∗
t
-
τ
3
V
1
V
C
1
∗
-
η
1
b
T
1
∗
1
+
β
V
1
+
β
V
1
-
1
-
ln
⁡
1
+
β
V
1
+
β
V
1
-
η
2
a
C
1
∗
1
+
β
V
1
+
β
V
1
-
1
-
l
n
1
+
β
V
1
+
β
V
1
+
r
η
3
V
1
-
μ
g
Z
=
-
d
T
-
T
1
2
T
-
1
-
ε
k
T
1
V
1
1
+
β
V
1
β
V
-
V
1
2
1
+
β
V
1
1
+
β
V
V
1
-
η
1
b
T
1
∗
H
T
1
T
+
H
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
T
1
∗
T
1
V
1
1
+
β
V
t
-
τ
1
T
∗
+
H
T
∗
t
-
τ
2
V
1
V
T
1
∗
+
H
1
+
β
V
1
+
β
V
1
-
η
2
a
C
1
∗
H
T
1
T
+
H
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
C
1
∗
T
1
V
1
1
+
β
V
t
-
τ
1
C
∗
+
H
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
H
1
+
β
V
1
+
β
V
1
+
r
η
3
d
g
+
μ
1
-
ε
k
+
d
β
μ
g
1
-
ε
k
+
d
β
g
R
1
-
1
Z
.We have if R
1
≤
1
<
R
0, then d
U
1
/
d
t
≤
0, where equality occurs at S
1. LIP implies global stability of S
1.Theorem 9.
For system (37), assume that R
1
>
1; then, S
2 is GAS.Proof.
We consider(52)
U
2
=
T
2
H
T
T
2
+
η
1
T
2
∗
H
T
∗
T
2
∗
+
η
2
C
2
∗
H
C
∗
C
2
∗
+
η
3
V
2
H
V
V
2
+
η
4
Z
2
H
Z
Z
2
+
η
5
T
2
V
2
1
+
β
V
2
∫
0
τ
1
H
T
t
-
θ
V
t
-
θ
1
+
β
V
2
T
2
V
2
1
+
β
V
t
-
θ
d
θ
+
η
6
T
2
∗
∫
0
τ
2
H
T
∗
t
-
θ
T
2
∗
d
θ
+
η
7
C
2
∗
∫
0
τ
3
H
C
∗
t
-
θ
C
2
∗
d
θ
.Calculate d
U
2
/
d
t along the trajectories of system (37) as(53)
d
U
2
d
t
=
1
-
T
2
T
π
-
d
T
-
1
-
ε
k
T
V
1
+
β
V
+
η
1
1
-
T
2
∗
T
∗
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
b
T
∗
+
η
2
1
-
C
2
∗
C
∗
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
a
C
∗
+
η
3
1
-
V
2
V
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
1
-
Z
2
Z
g
V
Z
-
μ
Z
+
η
5
T
V
1
+
β
V
-
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
+
T
2
V
2
1
+
β
V
2
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
6
T
∗
-
T
∗
t
-
τ
2
+
T
2
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
7
C
∗
-
C
∗
t
-
τ
3
+
C
2
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the steady state conditions for S
2,(54)
π
=
d
T
2
+
1
-
ε
k
T
2
V
2
1
+
β
V
2
,
1
-
ξ
1
-
ε
k
T
2
V
2
1
+
β
V
2
=
e
m
1
τ
1
b
T
2
∗
,
ξ
1
-
ε
k
T
2
V
2
1
+
β
V
2
=
e
m
1
τ
1
a
C
2
∗
,
c
V
2
=
e
-
m
2
τ
2
N
T
b
T
2
∗
+
e
-
m
3
τ
3
N
C
a
C
2
∗
-
r
V
2
Z
2
,we get (
1
-
ε
)
k
T
2
V
2
/
1
+
β
V
2
=
η
1
b
T
2
∗
+
η
2
a
C
2
∗, and(55)
d
U
2
d
t
=
-
d
T
-
T
2
2
T
+
1
-
ε
k
T
2
V
2
1
+
β
V
2
1
+
β
V
2
V
1
+
β
V
V
2
+
η
1
b
T
2
∗
1
-
T
2
T
+
η
2
a
C
2
∗
1
-
T
2
T
-
η
1
b
T
2
∗
T
t
-
τ
1
V
t
-
τ
1
T
2
∗
1
+
β
V
2
T
2
V
2
T
∗
1
+
β
V
t
-
τ
1
+
η
1
b
T
2
∗
-
η
2
a
C
2
∗
T
t
-
τ
1
V
t
-
τ
1
C
2
∗
1
+
β
V
2
T
2
V
2
C
∗
1
+
β
V
t
-
τ
1
+
η
2
a
C
2
∗
-
η
1
b
T
2
∗
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
η
2
a
C
2
∗
C
∗
t
-
τ
3
V
2
V
C
2
∗
+
η
1
b
T
2
∗
+
η
2
a
C
2
∗
-
1
-
ε
k
T
2
V
2
1
+
β
V
2
V
V
2
+
η
1
b
T
2
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
2
a
C
2
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
1
b
T
2
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
2
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using (50), when j
=
2, we get(56)
d
U
2
d
t
=
-
d
T
-
T
2
2
T
+
1
-
ε
k
T
2
V
2
1
+
β
V
2
1
+
β
V
2
V
1
+
β
V
V
2
-
V
V
2
-
1
+
1
+
β
V
1
+
β
V
2
-
η
1
b
T
2
∗
T
2
T
-
1
-
ln
⁡
T
2
T
-
η
2
a
C
2
∗
T
2
T
-
1
-
ln
⁡
T
2
T
-
η
1
b
T
2
∗
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
2
∗
T
2
V
2
T
∗
1
+
β
V
t
-
τ
1
-
1
-
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
2
∗
T
2
V
2
T
∗
1
+
β
V
t
-
τ
1
-
η
2
a
C
2
∗
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
C
2
∗
T
2
V
2
C
∗
1
+
β
V
t
-
τ
1
-
1
-
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
C
2
∗
T
2
V
2
C
∗
1
+
β
V
t
-
τ
1
-
η
1
b
T
2
∗
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
1
-
ln
⁡
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
η
2
a
C
2
∗
C
∗
t
-
τ
3
V
2
V
C
2
∗
-
1
-
ln
⁡
C
∗
t
-
τ
3
V
2
V
C
2
∗
-
η
1
b
T
2
∗
1
+
β
V
1
+
β
V
2
-
1
-
ln
⁡
1
+
β
V
1
+
β
V
2
-
η
2
a
C
2
∗
1
+
β
V
1
+
β
V
2
-
1
-
ln
⁡
1
+
β
V
1
+
β
V
2
=
-
d
T
-
T
2
2
T
-
1
-
ε
k
T
2
V
2
1
+
β
V
2
β
V
-
V
2
2
1
+
β
V
2
1
+
β
V
V
2
-
η
1
b
T
2
∗
H
T
2
T
+
H
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
2
∗
T
2
V
2
T
∗
1
+
β
V
t
-
τ
1
+
H
T
∗
t
-
τ
2
V
2
V
T
2
∗
+
H
1
+
β
V
1
+
β
V
2
-
η
2
a
C
2
∗
H
T
2
T
+
H
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
C
2
∗
T
2
V
2
C
∗
1
+
β
V
t
-
τ
1
+
H
C
∗
t
-
τ
3
V
2
V
C
2
∗
+
H
1
+
β
V
1
+
β
V
2
.Similar to the proof of Theorem 5, one can show that S
2 is GAS.


## 3.1. Steady States

System (37) always admits an uninfected steady state S
0
=
(
T
0
,
0,0
,
0,0
), where T
0
=
π
/
d. Now, we define the basic reproduction number for system (37) as(38)
R
0
=
1
-
ε
k
T
0
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
c
. The system has two other steady states S
1
(
T
1
,
T
1
∗
,
C
1
∗
,
V
1
,
0
) and S
2
(
T
2
,
T
2
∗
,
C
2
∗
,
V
2
,
Z
2
), where(39)
T
1
=
β
π
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
+
c
1
-
ε
k
+
d
β
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
,
T
1
∗
=
e
-
m
1
τ
1
1
-
ξ
c
d
b
1
-
ε
k
+
d
β
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
R
0
-
1
,
C
1
∗
=
e
-
m
1
τ
1
ξ
c
d
a
1
-
ε
k
+
d
β
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
R
0
-
1
,
V
1
=
d
1
-
ε
k
+
d
β
R
0
-
1
,
T
2
=
π
g
+
β
μ
g
d
+
1
-
ε
k
μ
+
d
β
μ
,
T
2
∗
=
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
π
μ
b
g
d
+
1
-
ε
k
μ
+
d
β
μ
,
C
2
∗
=
e
-
m
1
τ
1
ξ
1
-
ε
k
π
μ
a
g
d
+
1
-
ε
k
μ
+
d
β
μ
,
V
2
=
μ
g
,
Z
2
=
c
r
d
g
R
0
d
g
+
1
-
ε
k
μ
+
d
β
μ
-
1
.We note that T
1, T
1
∗, C
1
∗, and V
1 are positive when R
0
>
1. And T
2
,
T
2
∗
,
C
2
∗
,
V
2
>
0 and Z
2
>
0 when d
g
R
0
/
d
g
+
(
1
-
ε
)
k
μ
>
1. Now, we define another threshold parameter R
1 as(40)
R
1
=
R
0
1
+
1
-
ε
k
μ
+
d
β
μ
/
d
g
.Clearly, R
0
>
R
1. From above, we have the following result.Lemma 6.
For system (37), one has the following:(i)
IfR
0
≤
1, then the system has only one positive steady state S
0.
(ii)
IfR
1
≤
1
<
R
0, then the system has two positive steady states S
0 and S
1.
(iii)
IfR
1
>
1, then the system has three positive steady states S
0, S
1, and S
2.


## 3.2. Global Stability Analysis

In this subsection, we investigate the global stability of system (37) by constructing suitable Lyapunov functionals and applying LaSalle invariance principle.Theorem 7.
For system (37), if R
0
≤
1, then S
0 is GAS.Proof.
Define(41)
U
0
=
T
0
H
T
T
0
+
η
1
T
∗
+
η
2
C
∗
+
η
3
V
+
η
4
Z
+
η
5
∫
0
τ
1
T
t
-
θ
V
t
-
θ
1
+
β
V
t
-
θ
d
θ
+
η
6
∫
0
τ
2
T
∗
t
-
θ
d
θ
+
η
7
∫
0
τ
3
C
∗
t
-
θ
d
θ
.Calculating the derivative of U
0 along the solutions of system (37), we obtain(42)
d
U
0
d
t
=
1
-
T
0
T
π
-
1
-
ε
k
T
V
1
+
β
V
-
d
T
+
η
1
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
b
T
∗
+
η
2
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
a
C
∗
+
η
3
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
g
V
Z
-
μ
Z
+
η
5
T
V
1
+
β
V
-
η
5
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
+
η
6
T
∗
-
η
6
T
∗
t
-
τ
2
+
η
7
C
∗
-
η
7
C
∗
t
-
τ
3
.Collecting terms of (42), we obtain(43)
d
U
0
d
t
=
1
-
T
0
T
π
-
d
T
+
1
-
ε
k
T
0
V
1
+
β
V
-
η
3
c
V
-
η
4
μ
Z
=
-
d
T
-
T
0
2
T
+
η
3
c
β
R
0
V
2
1
+
β
V
+
η
4
μ
Z
+
c
η
3
R
0
-
1
V
.It follows that S
0 is GAS when R
0
≤
1.Theorem 8.
For system (37), assume that R
1
≤
1
<
R
0; then, S
1 is GAS.Proof.
Consider(44)
U
1
=
T
1
H
T
T
1
+
η
1
T
1
∗
H
T
∗
T
1
∗
+
η
2
C
1
∗
H
C
∗
C
1
∗
+
η
3
V
1
H
V
V
1
+
η
4
Z
+
η
5
T
1
V
1
1
+
β
V
1
∫
0
τ
1
H
T
t
-
θ
V
t
-
θ
1
+
β
V
1
T
1
V
1
1
+
β
V
t
-
θ
d
θ
+
η
6
T
1
∗
∫
0
τ
2
H
T
∗
t
-
θ
T
1
∗
d
θ
+
η
7
C
1
∗
∫
0
τ
3
H
C
∗
t
-
θ
C
1
∗
d
θ
.Then, d
U
1
/
d
t is given by(45)
d
U
1
d
t
=
1
-
T
1
T
π
-
d
T
-
1
-
ε
k
T
V
1
+
β
V
+
η
1
1
-
T
1
∗
T
∗
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
b
T
∗
+
η
2
1
-
C
1
∗
C
∗
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
a
C
∗
+
η
3
1
-
V
1
V
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
g
V
Z
-
μ
Z
+
η
5
T
V
1
+
β
V
-
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
+
T
1
V
1
1
+
β
V
1
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
6
T
∗
-
T
∗
t
-
τ
2
+
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
7
C
∗
-
C
∗
t
-
τ
3
+
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Simplifying (45), we get(46)
d
U
1
d
t
=
1
-
T
1
T
π
-
d
T
+
1
-
ε
k
T
1
V
1
+
β
V
-
e
m
1
τ
1
η
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
T
1
∗
T
∗
+
η
1
b
T
1
∗
-
e
-
m
1
τ
1
η
2
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
C
1
∗
C
∗
+
η
2
a
C
1
∗
-
b
η
1
T
∗
t
-
τ
2
V
1
V
-
a
η
2
C
∗
t
-
τ
3
V
1
V
-
c
η
3
V
+
c
η
3
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
+
η
5
T
1
V
1
1
+
β
V
1
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Applying π
=
d
T
1
+
(
1
-
ε
)
k
T
1
V
1
/
1
+
β
V
1, we get(47)
d
U
1
d
t
=
1
-
T
1
T
d
T
1
-
d
T
+
1
-
ε
k
T
1
V
1
1
+
β
V
1
1
-
T
1
T
+
1
-
ε
k
T
1
V
1
+
β
V
-
e
-
m
1
τ
1
η
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
T
1
∗
T
∗
+
η
1
b
T
1
∗
-
e
-
m
1
τ
1
η
2
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
C
1
∗
C
∗
+
η
2
a
C
1
∗
-
b
η
1
T
∗
t
-
τ
2
V
1
V
-
a
η
2
C
∗
t
-
τ
3
V
1
V
-
c
η
3
V
+
c
η
3
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
+
T
1
V
1
e
-
m
1
τ
1
η
1
1
-
ξ
1
-
ε
k
1
+
β
V
1
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
T
1
V
1
e
-
m
1
τ
1
η
2
ξ
1
-
ε
k
1
+
β
V
1
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the steady state conditions for S
1,(48)
1
-
ξ
1
-
ε
k
T
1
V
1
1
+
β
V
1
=
e
m
1
τ
1
b
T
1
∗
,
ξ
1
-
ε
k
T
1
V
1
1
+
β
V
1
=
e
m
1
τ
1
a
C
1
∗
,
c
V
1
=
e
-
m
2
τ
2
N
T
b
T
1
∗
+
e
-
m
3
τ
3
N
C
a
C
1
∗
,we obtain (
1
-
ε
)
k
T
1
V
1
/
1
+
β
V
1
=
η
1
b
T
1
∗
+
η
2
a
C
1
∗ and(49)
d
U
1
d
t
=
-
d
T
-
T
1
2
T
+
1
-
ε
k
T
1
V
1
1
+
β
V
1
1
+
β
V
1
V
1
+
β
V
V
1
+
η
1
b
T
1
∗
1
-
T
1
T
+
η
2
a
C
1
∗
1
-
T
1
T
-
η
1
b
T
1
∗
T
t
-
τ
1
V
t
-
τ
1
T
1
∗
1
+
β
V
1
T
1
V
1
T
∗
1
+
β
V
t
-
τ
1
+
η
1
b
T
1
∗
-
η
2
a
C
1
∗
T
t
-
τ
1
V
t
-
τ
1
C
1
∗
1
+
β
V
1
T
1
V
1
C
∗
1
+
β
V
t
-
τ
1
+
η
2
a
C
1
∗
-
η
1
b
T
1
∗
V
1
T
∗
t
-
τ
2
V
T
1
∗
-
η
2
a
C
1
∗
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
η
1
b
T
1
∗
+
η
2
a
C
1
∗
-
1
-
ε
k
T
1
V
1
1
+
β
V
1
V
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
+
η
1
b
T
1
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
2
a
C
1
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Consider the following equalities:(50)
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
=
ln
⁡
T
j
T
+
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
T
j
∗
T
j
V
j
1
+
β
V
t
-
τ
1
T
∗
+
ln
⁡
V
j
T
∗
V
T
j
∗
+
ln
⁡
1
+
β
V
1
+
β
V
j
,
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
=
ln
⁡
T
j
T
+
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
j
C
j
∗
T
j
V
j
1
+
β
V
t
-
τ
1
C
∗
+
ln
⁡
V
j
C
∗
V
C
j
∗
+
ln
⁡
1
+
β
V
1
+
β
V
j
,
ln
⁡
T
∗
t
-
τ
2
T
∗
=
ln
⁡
T
∗
t
-
τ
2
V
j
V
T
j
∗
+
ln
⁡
V
T
j
∗
V
j
T
∗
,
ln
⁡
C
∗
t
-
τ
3
C
∗
=
ln
⁡
C
∗
t
-
τ
3
V
j
V
C
j
∗
+
ln
⁡
V
C
j
∗
V
j
C
∗
,
j
=
1,2
.Using (50), when j
=
1, we get(51)
d
U
1
d
t
=
-
d
T
-
T
1
2
T
+
1
-
ε
k
T
1
V
1
1
+
β
V
1
1
+
β
V
1
V
1
+
β
V
V
1
-
V
V
1
-
1
+
1
+
β
V
1
+
β
V
1
-
η
1
b
T
1
∗
T
1
T
-
1
-
ln
⁡
T
1
T
-
η
2
a
C
1
∗
T
1
T
-
1
-
ln
⁡
T
1
T
-
η
1
b
T
1
∗
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
T
1
∗
T
1
V
1
1
+
β
V
t
-
τ
1
T
∗
-
1
-
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
T
1
∗
T
1
V
1
1
+
β
V
t
-
τ
1
T
∗
-
η
2
a
C
1
∗
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
C
1
∗
T
1
V
1
1
+
β
V
t
-
τ
1
C
∗
-
1
-
l
n
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
C
1
∗
T
1
V
1
1
+
β
V
t
-
τ
1
C
∗
-
η
1
b
T
1
∗
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
1
-
ln
⁡
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
η
2
a
C
1
∗
C
∗
t
-
τ
3
V
1
V
C
1
∗
-
1
-
ln
⁡
C
∗
t
-
τ
3
V
1
V
C
1
∗
-
η
1
b
T
1
∗
1
+
β
V
1
+
β
V
1
-
1
-
ln
⁡
1
+
β
V
1
+
β
V
1
-
η
2
a
C
1
∗
1
+
β
V
1
+
β
V
1
-
1
-
l
n
1
+
β
V
1
+
β
V
1
+
r
η
3
V
1
-
μ
g
Z
=
-
d
T
-
T
1
2
T
-
1
-
ε
k
T
1
V
1
1
+
β
V
1
β
V
-
V
1
2
1
+
β
V
1
1
+
β
V
V
1
-
η
1
b
T
1
∗
H
T
1
T
+
H
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
T
1
∗
T
1
V
1
1
+
β
V
t
-
τ
1
T
∗
+
H
T
∗
t
-
τ
2
V
1
V
T
1
∗
+
H
1
+
β
V
1
+
β
V
1
-
η
2
a
C
1
∗
H
T
1
T
+
H
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
1
C
1
∗
T
1
V
1
1
+
β
V
t
-
τ
1
C
∗
+
H
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
H
1
+
β
V
1
+
β
V
1
+
r
η
3
d
g
+
μ
1
-
ε
k
+
d
β
μ
g
1
-
ε
k
+
d
β
g
R
1
-
1
Z
.We have if R
1
≤
1
<
R
0, then d
U
1
/
d
t
≤
0, where equality occurs at S
1. LIP implies global stability of S
1.Theorem 9.
For system (37), assume that R
1
>
1; then, S
2 is GAS.Proof.
We consider(52)
U
2
=
T
2
H
T
T
2
+
η
1
T
2
∗
H
T
∗
T
2
∗
+
η
2
C
2
∗
H
C
∗
C
2
∗
+
η
3
V
2
H
V
V
2
+
η
4
Z
2
H
Z
Z
2
+
η
5
T
2
V
2
1
+
β
V
2
∫
0
τ
1
H
T
t
-
θ
V
t
-
θ
1
+
β
V
2
T
2
V
2
1
+
β
V
t
-
θ
d
θ
+
η
6
T
2
∗
∫
0
τ
2
H
T
∗
t
-
θ
T
2
∗
d
θ
+
η
7
C
2
∗
∫
0
τ
3
H
C
∗
t
-
θ
C
2
∗
d
θ
.Calculate d
U
2
/
d
t along the trajectories of system (37) as(53)
d
U
2
d
t
=
1
-
T
2
T
π
-
d
T
-
1
-
ε
k
T
V
1
+
β
V
+
η
1
1
-
T
2
∗
T
∗
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
b
T
∗
+
η
2
1
-
C
2
∗
C
∗
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
-
a
C
∗
+
η
3
1
-
V
2
V
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
1
-
Z
2
Z
g
V
Z
-
μ
Z
+
η
5
T
V
1
+
β
V
-
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
t
-
τ
1
+
T
2
V
2
1
+
β
V
2
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
6
T
∗
-
T
∗
t
-
τ
2
+
T
2
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
7
C
∗
-
C
∗
t
-
τ
3
+
C
2
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the steady state conditions for S
2,(54)
π
=
d
T
2
+
1
-
ε
k
T
2
V
2
1
+
β
V
2
,
1
-
ξ
1
-
ε
k
T
2
V
2
1
+
β
V
2
=
e
m
1
τ
1
b
T
2
∗
,
ξ
1
-
ε
k
T
2
V
2
1
+
β
V
2
=
e
m
1
τ
1
a
C
2
∗
,
c
V
2
=
e
-
m
2
τ
2
N
T
b
T
2
∗
+
e
-
m
3
τ
3
N
C
a
C
2
∗
-
r
V
2
Z
2
,we get (
1
-
ε
)
k
T
2
V
2
/
1
+
β
V
2
=
η
1
b
T
2
∗
+
η
2
a
C
2
∗, and(55)
d
U
2
d
t
=
-
d
T
-
T
2
2
T
+
1
-
ε
k
T
2
V
2
1
+
β
V
2
1
+
β
V
2
V
1
+
β
V
V
2
+
η
1
b
T
2
∗
1
-
T
2
T
+
η
2
a
C
2
∗
1
-
T
2
T
-
η
1
b
T
2
∗
T
t
-
τ
1
V
t
-
τ
1
T
2
∗
1
+
β
V
2
T
2
V
2
T
∗
1
+
β
V
t
-
τ
1
+
η
1
b
T
2
∗
-
η
2
a
C
2
∗
T
t
-
τ
1
V
t
-
τ
1
C
2
∗
1
+
β
V
2
T
2
V
2
C
∗
1
+
β
V
t
-
τ
1
+
η
2
a
C
2
∗
-
η
1
b
T
2
∗
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
η
2
a
C
2
∗
C
∗
t
-
τ
3
V
2
V
C
2
∗
+
η
1
b
T
2
∗
+
η
2
a
C
2
∗
-
1
-
ε
k
T
2
V
2
1
+
β
V
2
V
V
2
+
η
1
b
T
2
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
2
a
C
2
∗
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
V
1
+
β
V
t
-
τ
1
+
η
1
b
T
2
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
2
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using (50), when j
=
2, we get(56)
d
U
2
d
t
=
-
d
T
-
T
2
2
T
+
1
-
ε
k
T
2
V
2
1
+
β
V
2
1
+
β
V
2
V
1
+
β
V
V
2
-
V
V
2
-
1
+
1
+
β
V
1
+
β
V
2
-
η
1
b
T
2
∗
T
2
T
-
1
-
ln
⁡
T
2
T
-
η
2
a
C
2
∗
T
2
T
-
1
-
ln
⁡
T
2
T
-
η
1
b
T
2
∗
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
2
∗
T
2
V
2
T
∗
1
+
β
V
t
-
τ
1
-
1
-
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
2
∗
T
2
V
2
T
∗
1
+
β
V
t
-
τ
1
-
η
2
a
C
2
∗
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
C
2
∗
T
2
V
2
C
∗
1
+
β
V
t
-
τ
1
-
1
-
ln
⁡
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
C
2
∗
T
2
V
2
C
∗
1
+
β
V
t
-
τ
1
-
η
1
b
T
2
∗
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
1
-
ln
⁡
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
η
2
a
C
2
∗
C
∗
t
-
τ
3
V
2
V
C
2
∗
-
1
-
ln
⁡
C
∗
t
-
τ
3
V
2
V
C
2
∗
-
η
1
b
T
2
∗
1
+
β
V
1
+
β
V
2
-
1
-
ln
⁡
1
+
β
V
1
+
β
V
2
-
η
2
a
C
2
∗
1
+
β
V
1
+
β
V
2
-
1
-
ln
⁡
1
+
β
V
1
+
β
V
2
=
-
d
T
-
T
2
2
T
-
1
-
ε
k
T
2
V
2
1
+
β
V
2
β
V
-
V
2
2
1
+
β
V
2
1
+
β
V
V
2
-
η
1
b
T
2
∗
H
T
2
T
+
H
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
T
2
∗
T
2
V
2
T
∗
1
+
β
V
t
-
τ
1
+
H
T
∗
t
-
τ
2
V
2
V
T
2
∗
+
H
1
+
β
V
1
+
β
V
2
-
η
2
a
C
2
∗
H
T
2
T
+
H
T
t
-
τ
1
V
t
-
τ
1
1
+
β
V
C
2
∗
T
2
V
2
C
∗
1
+
β
V
t
-
τ
1
+
H
C
∗
t
-
τ
3
V
2
V
C
2
∗
+
H
1
+
β
V
1
+
β
V
2
.Similar to the proof of Theorem 5, one can show that S
2 is GAS.


## 4. Model with General Incidence

In the following model, we assume the incidence rate is given by a general function of the concentration of CD4+ T cells and viruses:(57)
T
˙
t
=
π
-
Φ
T
t
,
V
t
-
d
T
,
(58)
T
˙
∗
t
=
e
-
m
1
τ
1
1
-
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
b
T
∗
t
,
(59)
C
˙
∗
t
=
e
-
m
1
τ
1
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
a
C
∗
t
,
(60)
V
˙
t
=
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
t
Z
t
-
c
V
t
,
(61)
Z
˙
t
=
g
V
t
Z
t
-
μ
Z
t
.The incidence rate is given by Φ
(
T
(
t
)
,
V
(
t
)
) which is assumed to be continuously differentiable; moreover, it satisfies the following conditions.Condition C1. (i) Consider Φ
(
T
,
V
)
>
0 and Φ
(
0
,
V
)
=
Φ
(
T
,
0
)
=
0, for all T
,
V
∈
(
0
,
∞
). (ii) Consider ∂
Φ
(
T
,
V
)
/
∂
T
>
0
,
∂
Φ
(
T
,
V
)
/
∂
V
>
0, and ∂
Φ
(
T
,
0
)
/
∂
V
>
0, for all T
,
V
∈
(
0
,
∞
).Condition C2. (i) Consider Φ
(
T
,
V
)
⩽
V
∂
Φ
(
T
,
0
)
/
∂
V, for all T
,
V
∈
(
0
,
∞
). (ii) Consider d
/
d
T
∂
Φ
(
T
,
0
)
/
∂
V
>
0
.The nonnegativity and boundedness of the solutions of model (57)–(61) can be shown as given in Section 2.


### 4.1. Steady States

Lemma 10.
Suppose that Conditions C1 and C2 are satisfied; then, there exist two bifurcation parametersR
0 and R
1 with R
0
>
R
1
>
0. Moreover,(i)
ifR
0
≤
1, then the system has only one positive steady state S
0,
(ii)
ifR
1
≤
1
<
R
0, then the system has two positive steady states S
0 and S
1,
(iii)
ifR
1
>
1, then the system has three positive steady states S
0, S
1, and S
2.Proof.
Let(62)
π
-
Φ
T
,
V
-
d
T
=
0
,
(63)
e
-
m
1
τ
1
1
-
ξ
Φ
T
,
V
-
b
T
∗
=
0
,
(64)
e
-
m
1
τ
1
ξ
Φ
T
,
V
-
a
C
∗
=
0
,
(65)
e
-
m
2
τ
2
N
T
b
T
∗
+
e
-
m
3
τ
3
N
C
a
C
∗
-
r
V
Z
-
c
V
=
0
,
(66)
g
V
Z
-
μ
Z
=
0
.Equation (66) admits two solutions, Z
=
0 and V
=
μ
/
g. Substituting Z
=
0 into (63) and (64), we obtain T
∗ and C
∗ as(67)
T
∗
=
e
-
m
1
τ
1
1
-
ξ
Φ
T
,
V
b
,
C
∗
=
e
-
m
1
τ
1
ξ
Φ
T
,
V
a
.Substituting (67) into (65), we get(68)
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
Φ
T
,
V
-
c
V
=
0
.Using Condition C1, we have V
=
0 is a solution of (68). Therefore, T
∗
=
C
∗
=
0 and T
=
T
0 which leads to S
0
=
(
T
0
,
0,0
,
0,0
). If V
≠
0, then, from (62) and (68), we obtain(69)
V
=
γ
Φ
T
,
V
c
=
γ
π
-
d
T
c
⟹
T
=
T
0
-
c
V
d
γ
,where γ
=
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C. Then, (68) becomes(70)
γ
Φ
T
0
-
c
d
γ
V
,
V
-
c
V
=
0
. Let us define a function Ψ
1 as(71)
Ψ
1
V
=
γ
Φ
T
0
-
c
d
γ
V
,
V
-
c
V
=
0
.Condition C1 implies that Ψ
1
(
0
)
=
0, and when V
=
V
∗
=
π
γ
/
c
>
0, then Ψ
1
(
V
∗
)
=
-
c
V
∗
<
0. We have(72)
Ψ
1
′
0
=
γ
-
c
d
γ
∂
Φ
T
0
,
0
∂
T
+
∂
Φ
T
0
,
0
∂
V
-
c
. From Condition C1, we have ∂
Φ
T
0
,
0
/
∂
T
=
0; then,(73)
Ψ
1
′
0
=
γ
∂
Φ
T
0
,
0
∂
V
-
c
=
c
γ
c
∂
Φ
T
0
,
0
∂
V
-
1
. Therefore, if Ψ
1
′
(
0
)
>
0, that is, γ
/
c
∂
Φ
T
0
,
0
/
∂
V
>
1, then there exists V
1
∈
(
0
,
V
∗
) such that Ψ
1
(
V
1
)
=
0. From (62), we define a function Ψ
2 as(74)
Ψ
2
T
=
π
-
d
T
-
Φ
T
,
V
1
=
0
. Using Condition C1, we have Ψ
2
(
0
)
=
π
>
0 and Ψ
2
(
T
0
)
=
-
Φ
(
T
0
,
V
1
)
<
0. Since Ψ
2 is a strictly decreasing function, then there exists a unique T
1
∈
(
0
,
T
0
) such that Ψ
2
(
T
1
)
=
0. It follows that T
1
∗
=
e
-
m
1
τ
1
(
1
-
ξ
)
Φ
T
1
,
V
1
/
b
>
0 and C
1
∗
=
e
-
m
1
τ
1
ξ
Φ
T
1
,
V
1
/
a
>
0. It means that an infected steady state without humoral immune response S
1
(
T
1
,
T
1
∗
,
C
1
∗
,
V
1
,
0
) exists when γ
/
c
∂
Φ
T
0
,
0
/
∂
V
>
1. Now, we define the parameter R
0 as(75)
R
0
=
γ
c
∂
Φ
T
0
,
0
∂
V
.The other possibility of (66) is V
2
=
μ
/
g. Insert V
2 in (62) and define Ψ
3 as(76)
Ψ
3
T
=
π
-
d
T
-
Φ
T
,
V
2
=
0
. Applying Condition C1, then Ψ
3 is a strictly decreasing function. Moreover, Ψ
3
(
0
)
=
π
>
0 and Ψ
3
(
T
0
)
=
-
Φ
(
T
0
,
V
2
)
<
0. Thus, there exists a unique T
2
∈
(
0
,
T
0
) such that Ψ
3
(
T
2
)
=
0. It follows from (63)–(65) that(77)
T
2
∗
=
e
-
m
1
τ
1
1
-
ξ
Φ
T
2
,
V
2
b
,
C
2
∗
=
e
-
m
1
τ
1
ξ
Φ
T
2
,
V
2
a
,
Z
2
=
c
r
γ
Φ
T
2
,
V
2
c
V
2
-
1
.Thus, T
2
∗
>
0 and C
2
∗
>
0; moreover, Z
2
>
0 when γ
Φ
T
2
,
V
2
/
c
V
2
>
1. Now, we define parameter R
1 as(78)
R
1
=
γ
Φ
T
2
,
V
2
c
V
2
. Hence, Z
2 can be rewritten as Z
2
=
c
/
r
R
1
-
1. It follows that S
2
(
T
2
,
T
2
∗
,
C
2
∗
,
V
2
,
Z
2
) exists when R
1
>
1.
Conditions C1 and C2 imply that(79)
R
1
=
γ
Φ
T
2
,
V
2
c
V
2
⩽
γ
c
∂
Φ
T
2
,
0
∂
V
<
γ
c
∂
Φ
T
0
,
0
∂
V
=
R
0
.


### 4.2. Global Stability Analysis

Theorem 11.
Let Conditions C1 and C2 be satisfied andR
0
≤
1; then, S
0 for system (57)–(61) is GAS.Proof.
Define(80)
U
0
=
T
-
T
0
-
∫
T
0
T
lim
V
→
0
+
⁡
Φ
T
0
,
V
Φ
ρ
,
V
d
ρ
+
η
1
T
∗
+
η
2
C
∗
+
η
3
V
+
η
4
Z
+
∫
0
τ
1
Φ
T
t
-
θ
,
V
t
-
θ
d
θ
+
η
6
∫
0
τ
2
T
∗
t
-
θ
d
θ
+
η
7
∫
0
τ
3
C
∗
t
-
θ
d
θ
.Then,(81)
d
U
0
d
t
=
1
-
lim
V
→
0
+
⁡
Φ
T
0
,
V
Φ
T
,
V
π
-
d
T
-
Φ
T
,
V
+
η
1
e
-
m
1
τ
1
1
-
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
b
T
∗
+
η
2
e
-
m
1
τ
1
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
a
C
∗
+
η
3
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
g
V
Z
-
μ
Z
+
Φ
T
,
V
-
Φ
T
t
-
τ
1
,
V
t
-
τ
1
+
η
6
T
∗
-
η
6
T
∗
t
-
τ
2
+
η
7
C
∗
-
η
7
C
∗
t
-
τ
3
.Equation (81) can be simplified as(82)
d
U
0
d
t
=
1
-
lim
V
→
0
+
⁡
Φ
T
0
,
V
Φ
T
,
V
π
-
d
T
+
Φ
T
,
V
lim
V
→
0
+
⁡
Φ
T
0
,
V
Φ
T
,
V
-
c
η
3
V
-
η
4
μ
Z
=
π
1
-
∂
Φ
T
0
,
0
/
∂
V
∂
Φ
T
,
0
/
∂
V
1
-
T
T
0
+
Φ
T
,
V
∂
Φ
T
0
,
0
/
∂
V
∂
Φ
T
,
0
/
∂
V
-
c
η
3
V
-
η
4
μ
Z
⩽
π
1
-
∂
Φ
T
0
,
0
/
∂
V
∂
Φ
T
,
0
/
∂
V
1
-
T
T
0
+
V
∂
Φ
T
0
,
0
∂
V
-
c
η
3
V
-
η
4
μ
Z
=
π
1
-
∂
Φ
T
0
,
0
/
∂
V
∂
Φ
T
,
0
/
∂
V
1
-
T
T
0
+
c
η
3
1
c
η
3
∂
Φ
T
0
,
0
∂
V
-
1
V
-
η
4
μ
Z
=
π
1
-
∂
Φ
T
0
,
0
/
∂
V
∂
Φ
T
,
0
/
∂
V
1
-
T
T
0
+
c
η
3
R
0
-
1
V
-
η
4
μ
Z
.Based on Condition C2, the first term of (82) is less than or equal to zero. Therefore, if R
0
≤
1, then d
U
0
/
d
t
≤
0, for all T
,
V
,
Z
>
0. Similar to the previous sections, one can show that S
0 is GAS.Condition C3. Consider(83)
Φ
T
,
V
Φ
T
,
V
i
-
V
V
i
1
-
Φ
T
,
V
i
Φ
T
,
V
⩽
0
,
∀
T
,
V
∈
0
,
∞
,
i
=
1,2
.Lemma 12.
Suppose thatR
0
>
1 and Conditions C1–C3 hold. Then, T
1
,
T
2
,
T
1
∗
,
T
2
∗
,
C
1
∗
,
C
2
∗
,
V
1
,
V
2 exist satisfying(84)
sgn
⁡
R
1
-
1
=
sgn
⁡
T
2
-
T
1
=
sgn
⁡
V
1
-
V
2
.Proof.
From Condition C1, forT
1
,
T
2
,
V
1
,
V
2
>
0, we have(85)
Φ
T
2
,
V
2
-
Φ
T
1
,
V
2
T
2
-
T
1
>
0
.
(86)
Φ
T
1
,
V
2
-
Φ
T
1
,
V
1
V
2
-
V
1
>
0
,Applying Condition C3 when i
=
1,  T
=
T
1, and V
=
V
2, we get(87)
Φ
T
1
,
V
2
V
1
-
Φ
T
1
,
V
1
V
2
Φ
T
1
,
V
2
-
Φ
T
1
,
V
1
⩽
0
.It follows from inequality (86) that (88)
Φ
T
1
,
V
2
V
1
-
Φ
T
1
,
V
1
V
2
V
1
-
V
2
>
0
.We assume that sgn
⁡
V
2
-
V
1
=
sgn
⁡
T
2
-
T
1. Then, we have(89)
π
-
d
T
2
-
π
-
d
T
1
=
Φ
T
2
,
V
2
-
Φ
T
1
,
V
1
=
Φ
T
2
,
V
2
-
Φ
T
1
,
V
2
+
Φ
T
1
,
V
2
-
Φ
T
1
,
V
1
,and using inequalities (85) and (86), we obtain sgn
⁡
T
1
-
T
2
=
sgn
⁡
T
2
-
T
1, and this leads to contradiction. Therefore, sgn
⁡
V
1
-
V
2
=
sgn
⁡
T
2
-
T
1. From the conditions of S
1, we have γ
Φ
T
1
,
V
1
/
c
V
1
=
1; then,(90)
R
1
-
1
=
γ
c
1
V
2
Φ
T
2
,
V
2
-
Φ
T
1
,
V
2
+
1
V
1
V
2
Φ
T
1
,
V
2
V
1
-
Φ
T
1
,
V
1
V
2
.Using inequalities (85) and (88), we obtain sgn
⁡
R
1
-
1
=
sgn
⁡
V
1
-
V
2.Theorem 13.
Assume thatR
1
≤
1
<
R
0 and Conditions C1–C3 are satisfied; then, S
1 for system (57)–(61) is GAS.Proof.
Let us define(91)
U
1
=
T
-
T
1
-
∫
T
1
T
Φ
T
1
,
V
1
Φ
ρ
,
V
1
d
ρ
+
η
1
T
1
∗
H
T
∗
T
1
∗
+
η
2
C
1
∗
H
C
∗
C
1
∗
+
η
3
V
1
H
V
V
1
+
η
4
Z
+
Φ
T
1
,
V
1
∫
0
τ
1
H
Φ
T
t
-
θ
,
V
t
-
θ
Φ
T
1
,
V
1
d
θ
+
η
6
T
1
∗
∫
0
τ
2
H
T
∗
t
-
θ
T
1
∗
d
θ
+
η
7
C
1
∗
∫
0
τ
3
H
C
∗
t
-
θ
C
1
∗
d
θ
.Then,(92)
d
U
1
d
t
=
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
π
-
d
T
-
Φ
T
,
V
+
η
1
1
-
T
1
∗
T
∗
e
-
m
1
τ
1
1
-
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
b
T
∗
+
η
2
1
-
C
1
∗
C
∗
e
-
m
1
τ
1
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
a
C
∗
+
η
3
1
-
V
1
V
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
g
V
Z
-
μ
Z
+
Φ
T
,
V
-
Φ
T
t
-
τ
1
,
V
t
-
τ
1
+
Φ
T
1
,
V
1
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
6
T
∗
-
T
∗
t
-
τ
2
+
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
7
C
∗
-
C
∗
t
-
τ
3
+
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Collecting terms of (92), we obtain (93)
d
U
1
d
t
=
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
π
-
d
T
+
Φ
T
,
V
Φ
T
1
,
V
1
Φ
T
,
V
1
-
e
-
m
1
τ
1
η
1
1
-
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
1
∗
T
∗
+
η
1
b
T
1
∗
-
e
-
m
1
τ
1
η
2
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
1
∗
C
∗
+
η
2
a
C
1
∗
-
b
η
1
T
∗
t
-
τ
2
V
1
V
-
a
η
2
C
∗
t
-
τ
3
V
1
V
-
c
η
3
V
+
c
η
3
V
1
+
η
3
r
V
1
Z
-
η
4
μ
Z
+
Φ
T
1
,
V
1
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Applying π
=
d
T
1
+
Φ
(
T
1
,
V
1
), we get(94)
d
U
1
d
t
=
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
d
T
1
-
d
T
+
Φ
T
1
,
V
1
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
+
Φ
T
,
V
Φ
T
1
,
V
1
Φ
T
,
V
1
-
e
-
m
1
τ
1
η
1
1
-
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
1
∗
T
∗
+
η
1
b
T
1
∗
-
e
-
m
1
τ
1
η
2
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
1
∗
C
∗
+
η
2
a
C
1
∗
-
η
1
b
T
∗
t
-
τ
2
V
1
V
-
η
2
a
C
∗
t
-
τ
3
V
1
V
-
c
η
3
V
+
c
η
3
V
1
+
η
3
r
V
1
Z
-
η
4
μ
Z
+
Φ
T
1
,
V
1
e
-
m
1
τ
1
η
1
1
-
ξ
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
Φ
T
1
,
V
1
e
-
m
1
τ
1
η
2
ξ
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the steady state conditions for S
1,(95)
1
-
ξ
Φ
T
1
,
V
1
=
e
m
1
τ
1
b
T
1
∗
,
ξ
Φ
T
1
,
V
1
=
e
m
1
τ
1
a
C
1
∗
,
c
V
1
=
e
-
m
2
τ
2
N
T
b
T
1
∗
+
e
-
m
3
τ
3
N
C
a
C
1
∗
,then we have Φ
(
T
1
,
V
1
)
=
η
1
b
T
1
∗
+
η
2
a
C
1
∗ and (96)
d
U
1
d
t
=
d
T
1
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
1
-
T
T
1
+
Φ
T
1
,
V
1
Φ
T
,
V
Φ
T
,
V
1
+
η
1
b
T
1
∗
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
+
η
2
a
C
1
∗
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
-
η
1
b
T
1
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
1
∗
Φ
T
1
,
V
1
T
∗
+
η
1
b
T
1
∗
-
η
2
a
C
1
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
1
∗
Φ
T
1
,
V
1
C
∗
+
η
2
a
C
1
∗
-
η
1
b
T
1
∗
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
η
2
a
C
1
∗
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
η
1
b
T
1
∗
+
η
2
a
C
1
∗
-
Φ
T
1
,
V
1
V
V
1
+
η
3
r
V
1
-
μ
g
Z
+
η
1
b
T
1
∗
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
2
a
C
1
∗
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the following equalities,(97)
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
=
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
j
∗
Φ
T
j
,
V
j
T
∗
+
ln
⁡
Φ
T
j
,
V
j
Φ
T
,
V
j
+
ln
⁡
Φ
T
,
V
j
V
Φ
T
,
V
V
j
+
ln
⁡
V
j
T
∗
V
T
j
∗
,
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
=
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
j
∗
Φ
T
j
,
V
j
C
∗
+
ln
⁡
Φ
T
j
,
V
j
Φ
T
,
V
j
ln
⁡
Φ
T
,
V
j
V
Φ
T
,
V
V
j
+
ln
⁡
V
j
C
∗
V
C
j
∗
,
ln
⁡
T
∗
t
-
τ
2
T
∗
=
ln
⁡
T
∗
t
-
τ
2
V
j
V
T
j
∗
+
ln
⁡
V
T
j
∗
V
j
T
∗
,
ln
⁡
C
∗
t
-
τ
3
C
∗
=
ln
⁡
C
∗
t
-
τ
3
V
j
V
C
j
∗
+
ln
⁡
V
C
j
∗
V
j
C
∗
,
j
=
1,2
,in case of j
=
1, we obtain(98)
d
U
1
d
t
=
d
T
1
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
1
-
T
T
1
+
Φ
T
1
,
V
1
Φ
T
,
V
Φ
T
,
V
1
-
V
V
1
-
1
+
Φ
T
,
V
1
V
Φ
T
,
V
V
1
-
η
1
b
T
1
∗
Φ
T
1
,
V
1
Φ
T
,
V
1
-
1
-
ln
⁡
Φ
T
1
,
V
1
Φ
T
,
V
1
-
η
2
a
C
1
∗
Φ
T
1
,
V
1
Φ
T
,
V
1
-
1
-
ln
⁡
Φ
T
1
,
V
1
Φ
T
,
V
1
-
η
1
b
T
1
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
1
∗
Φ
T
1
,
V
1
T
∗
-
1
-
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
1
∗
Φ
T
1
,
V
1
T
∗
-
η
2
a
C
1
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
1
∗
Φ
T
1
,
V
1
C
∗
-
1
-
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
1
∗
Φ
T
1
,
V
1
C
∗
-
η
1
b
T
1
∗
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
1
-
ln
⁡
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
η
2
a
C
1
∗
C
∗
t
-
τ
3
V
1
V
C
1
∗
-
1
-
ln
⁡
C
∗
t
-
τ
3
V
1
V
C
1
∗
-
η
1
b
T
1
∗
Φ
T
,
V
1
V
Φ
T
,
V
V
1
-
1
-
ln
⁡
Φ
T
,
V
1
V
Φ
T
,
V
V
1
-
η
2
a
C
1
∗
Φ
T
,
V
1
V
Φ
T
,
V
V
1
-
1
-
ln
⁡
Φ
T
,
V
1
V
Φ
T
,
V
V
1
+
r
η
3
V
1
-
μ
g
Z
=
d
T
1
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
1
-
T
T
1
+
Φ
T
1
,
V
1
Φ
T
,
V
Φ
T
,
V
1
-
V
V
1
1
-
Φ
T
,
V
1
Φ
T
,
V
-
η
1
b
T
1
∗
H
Φ
T
1
,
V
1
Φ
T
,
V
1
+
H
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
1
∗
Φ
T
1
,
V
1
T
∗
+
H
T
∗
t
-
τ
2
V
1
V
T
1
∗
+
H
Φ
T
,
V
1
V
Φ
T
,
V
V
1
-
η
2
a
C
1
∗
H
Φ
T
1
,
V
1
Φ
T
,
V
1
+
H
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
1
∗
Φ
T
1
,
V
1
C
∗
+
H
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
H
Φ
T
,
V
1
V
Φ
T
,
V
V
1
+
r
η
3
V
1
-
μ
g
Z
.Since R
1
≤
1, then Conditions C1 and C3 and Lemma 12 imply that d
U
1
/
d
t
≤
0, for all T
,
T
∗
,
C
∗
,
V
,
Z
>
0, where the equality occurs at S
1. LIP implies that S
1 is GAS.Theorem 14.
LetR
1
>
1 and Conditions C1–C3 be satisfied; then, S
2 for system (57)–(61) is GAS.Proof.
We consider a Lyapunov functional:(99)
U
2
=
T
-
T
2
-
∫
T
2
T
Φ
T
2
,
V
2
Φ
ρ
,
V
2
d
ρ
+
η
1
T
2
∗
H
T
∗
T
2
∗
+
η
2
C
2
∗
H
C
∗
C
2
∗
+
η
3
V
2
H
V
V
2
+
η
4
Z
2
H
Z
Z
2
+
Φ
T
2
,
V
2
∫
0
τ
1
H
Φ
T
t
-
θ
,
V
t
-
θ
Φ
T
2
,
V
2
d
θ
+
η
6
T
2
∗
∫
0
τ
2
H
T
∗
t
-
θ
T
2
∗
d
θ
+
η
7
C
2
∗
∫
0
τ
3
H
C
∗
t
-
θ
C
2
∗
d
θ
.Then,(100)
d
U
2
d
t
=
1
-
Φ
T
2
,
V
2
Φ
T
,
V
2
π
-
d
T
-
Φ
T
,
V
+
η
1
1
-
T
2
∗
T
∗
e
-
m
1
τ
1
1
-
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
b
T
∗
+
η
2
1
-
C
2
∗
C
∗
e
-
m
1
τ
1
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
a
C
∗
+
η
3
1
-
V
2
V
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
c
V
-
r
V
Z
+
η
4
1
-
Z
2
Z
g
V
Z
-
μ
Z
+
Φ
T
,
V
-
Φ
T
t
-
τ
1
,
V
t
-
τ
1
+
Φ
T
2
,
V
2
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
6
T
∗
-
T
∗
t
-
τ
2
+
T
2
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
7
C
∗
-
C
∗
t
-
τ
3
+
C
2
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the steady state conditions for S
2,(101)
π
=
d
T
2
+
Φ
T
2
,
V
2
,
1
-
ξ
Φ
T
2
,
V
2
=
e
m
1
τ
1
b
T
2
∗
,
ξ
Φ
T
2
,
V
2
=
e
m
1
τ
1
a
C
2
∗
,
c
V
2
=
e
-
m
2
τ
2
N
T
b
T
2
∗
+
e
-
m
3
τ
3
N
C
a
C
2
∗
-
r
V
2
Z
2
, we get Φ
(
T
2
,
V
2
)
=
η
1
b
T
2
∗
+
η
2
a
C
2
∗, and (102)
d
U
2
d
t
=
d
T
2
1
-
Φ
T
2
,
V
2
Φ
T
,
V
2
1
-
T
T
2
+
Φ
T
2
,
V
2
Φ
T
,
V
Φ
T
,
V
2
+
η
1
b
T
2
∗
1
-
Φ
T
2
,
V
2
Φ
T
,
V
2
+
η
2
a
C
2
∗
1
-
Φ
T
2
,
V
2
Φ
T
,
V
2
-
η
1
b
T
2
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
2
∗
Φ
T
2
,
V
2
T
∗
+
η
1
b
T
2
∗
-
η
2
a
C
2
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
2
∗
Φ
T
2
,
V
2
C
∗
+
η
2
a
C
2
∗
-
η
1
b
T
2
∗
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
η
2
a
C
2
∗
C
∗
t
-
τ
3
V
2
V
C
2
∗
+
η
1
b
T
2
∗
+
η
2
a
C
2
∗
-
Φ
T
2
,
V
2
V
V
2
+
η
1
b
T
2
∗
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
2
a
C
2
∗
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
1
b
T
2
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
2
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using (97), when j
=
2, we get(103)
d
U
2
d
t
=
d
T
2
1
-
Φ
T
2
,
V
2
Φ
T
,
V
2
1
-
T
T
2
+
Φ
T
2
,
V
2
Φ
T
,
V
Φ
T
,
V
2
-
V
V
2
-
1
+
Φ
T
,
V
2
V
Φ
T
,
V
V
2
-
η
1
b
T
2
∗
Φ
T
2
,
V
2
Φ
T
,
V
2
-
1
-
ln
⁡
Φ
T
2
,
V
2
Φ
T
,
V
2
-
η
2
a
C
2
∗
Φ
T
2
,
V
2
Φ
T
,
V
2
-
1
-
ln
⁡
Φ
T
2
,
V
2
Φ
T
,
V
2
-
η
1
b
T
2
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
2
∗
Φ
T
2
,
V
2
T
∗
-
1
-
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
2
∗
Φ
T
2
,
V
2
T
∗
-
η
2
a
C
2
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
2
∗
Φ
T
2
,
V
2
C
∗
-
1
-
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
2
∗
Φ
T
2
,
V
2
C
∗
-
η
1
b
T
2
∗
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
1
-
ln
⁡
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
η
2
a
C
2
∗
C
∗
t
-
τ
3
V
2
V
C
2
∗
-
1
-
ln
⁡
C
∗
t
-
τ
3
V
2
V
C
2
∗
-
η
1
b
T
2
∗
Φ
T
,
V
2
V
Φ
T
,
V
V
2
-
1
-
ln
⁡
Φ
T
,
V
2
V
Φ
T
,
V
V
2
-
η
2
a
C
2
∗
Φ
T
,
V
2
V
Φ
T
,
V
V
2
-
1
-
ln
⁡
Φ
T
,
V
2
V
Φ
T
,
V
V
2
=
d
T
2
1
-
Φ
T
2
,
V
2
Φ
T
,
V
2
1
-
T
T
2
+
Φ
T
2
,
V
2
Φ
T
,
V
Φ
T
,
V
2
-
V
V
2
1
-
Φ
T
,
V
2
Φ
T
,
V
-
η
1
b
T
2
∗
H
Φ
T
2
,
V
2
Φ
T
,
V
2
+
H
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
2
∗
Φ
T
2
,
V
2
T
∗
+
H
T
∗
t
-
τ
2
V
2
V
T
2
∗
+
H
Φ
T
,
V
2
V
Φ
T
,
V
V
2
-
η
2
a
C
2
∗
H
Φ
T
1
,
V
1
Φ
T
,
V
1
+
H
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
2
∗
Φ
T
2
,
V
2
C
∗
+
H
C
∗
t
-
τ
3
V
2
V
C
2
∗
+
H
Φ
T
,
V
2
V
Φ
T
,
V
V
2
.It is easy to show that S
2 is GAS.


## 4.1. Steady States

Lemma 10.
Suppose that Conditions C1 and C2 are satisfied; then, there exist two bifurcation parametersR
0 and R
1 with R
0
>
R
1
>
0. Moreover,(i)
ifR
0
≤
1, then the system has only one positive steady state S
0,
(ii)
ifR
1
≤
1
<
R
0, then the system has two positive steady states S
0 and S
1,
(iii)
ifR
1
>
1, then the system has three positive steady states S
0, S
1, and S
2.Proof.
Let(62)
π
-
Φ
T
,
V
-
d
T
=
0
,
(63)
e
-
m
1
τ
1
1
-
ξ
Φ
T
,
V
-
b
T
∗
=
0
,
(64)
e
-
m
1
τ
1
ξ
Φ
T
,
V
-
a
C
∗
=
0
,
(65)
e
-
m
2
τ
2
N
T
b
T
∗
+
e
-
m
3
τ
3
N
C
a
C
∗
-
r
V
Z
-
c
V
=
0
,
(66)
g
V
Z
-
μ
Z
=
0
.Equation (66) admits two solutions, Z
=
0 and V
=
μ
/
g. Substituting Z
=
0 into (63) and (64), we obtain T
∗ and C
∗ as(67)
T
∗
=
e
-
m
1
τ
1
1
-
ξ
Φ
T
,
V
b
,
C
∗
=
e
-
m
1
τ
1
ξ
Φ
T
,
V
a
.Substituting (67) into (65), we get(68)
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
Φ
T
,
V
-
c
V
=
0
.Using Condition C1, we have V
=
0 is a solution of (68). Therefore, T
∗
=
C
∗
=
0 and T
=
T
0 which leads to S
0
=
(
T
0
,
0,0
,
0,0
). If V
≠
0, then, from (62) and (68), we obtain(69)
V
=
γ
Φ
T
,
V
c
=
γ
π
-
d
T
c
⟹
T
=
T
0
-
c
V
d
γ
,where γ
=
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C. Then, (68) becomes(70)
γ
Φ
T
0
-
c
d
γ
V
,
V
-
c
V
=
0
. Let us define a function Ψ
1 as(71)
Ψ
1
V
=
γ
Φ
T
0
-
c
d
γ
V
,
V
-
c
V
=
0
.Condition C1 implies that Ψ
1
(
0
)
=
0, and when V
=
V
∗
=
π
γ
/
c
>
0, then Ψ
1
(
V
∗
)
=
-
c
V
∗
<
0. We have(72)
Ψ
1
′
0
=
γ
-
c
d
γ
∂
Φ
T
0
,
0
∂
T
+
∂
Φ
T
0
,
0
∂
V
-
c
. From Condition C1, we have ∂
Φ
T
0
,
0
/
∂
T
=
0; then,(73)
Ψ
1
′
0
=
γ
∂
Φ
T
0
,
0
∂
V
-
c
=
c
γ
c
∂
Φ
T
0
,
0
∂
V
-
1
. Therefore, if Ψ
1
′
(
0
)
>
0, that is, γ
/
c
∂
Φ
T
0
,
0
/
∂
V
>
1, then there exists V
1
∈
(
0
,
V
∗
) such that Ψ
1
(
V
1
)
=
0. From (62), we define a function Ψ
2 as(74)
Ψ
2
T
=
π
-
d
T
-
Φ
T
,
V
1
=
0
. Using Condition C1, we have Ψ
2
(
0
)
=
π
>
0 and Ψ
2
(
T
0
)
=
-
Φ
(
T
0
,
V
1
)
<
0. Since Ψ
2 is a strictly decreasing function, then there exists a unique T
1
∈
(
0
,
T
0
) such that Ψ
2
(
T
1
)
=
0. It follows that T
1
∗
=
e
-
m
1
τ
1
(
1
-
ξ
)
Φ
T
1
,
V
1
/
b
>
0 and C
1
∗
=
e
-
m
1
τ
1
ξ
Φ
T
1
,
V
1
/
a
>
0. It means that an infected steady state without humoral immune response S
1
(
T
1
,
T
1
∗
,
C
1
∗
,
V
1
,
0
) exists when γ
/
c
∂
Φ
T
0
,
0
/
∂
V
>
1. Now, we define the parameter R
0 as(75)
R
0
=
γ
c
∂
Φ
T
0
,
0
∂
V
.The other possibility of (66) is V
2
=
μ
/
g. Insert V
2 in (62) and define Ψ
3 as(76)
Ψ
3
T
=
π
-
d
T
-
Φ
T
,
V
2
=
0
. Applying Condition C1, then Ψ
3 is a strictly decreasing function. Moreover, Ψ
3
(
0
)
=
π
>
0 and Ψ
3
(
T
0
)
=
-
Φ
(
T
0
,
V
2
)
<
0. Thus, there exists a unique T
2
∈
(
0
,
T
0
) such that Ψ
3
(
T
2
)
=
0. It follows from (63)–(65) that(77)
T
2
∗
=
e
-
m
1
τ
1
1
-
ξ
Φ
T
2
,
V
2
b
,
C
2
∗
=
e
-
m
1
τ
1
ξ
Φ
T
2
,
V
2
a
,
Z
2
=
c
r
γ
Φ
T
2
,
V
2
c
V
2
-
1
.Thus, T
2
∗
>
0 and C
2
∗
>
0; moreover, Z
2
>
0 when γ
Φ
T
2
,
V
2
/
c
V
2
>
1. Now, we define parameter R
1 as(78)
R
1
=
γ
Φ
T
2
,
V
2
c
V
2
. Hence, Z
2 can be rewritten as Z
2
=
c
/
r
R
1
-
1. It follows that S
2
(
T
2
,
T
2
∗
,
C
2
∗
,
V
2
,
Z
2
) exists when R
1
>
1.
Conditions C1 and C2 imply that(79)
R
1
=
γ
Φ
T
2
,
V
2
c
V
2
⩽
γ
c
∂
Φ
T
2
,
0
∂
V
<
γ
c
∂
Φ
T
0
,
0
∂
V
=
R
0
.


## 4.2. Global Stability Analysis

Theorem 11.
Let Conditions C1 and C2 be satisfied andR
0
≤
1; then, S
0 for system (57)–(61) is GAS.Proof.
Define(80)
U
0
=
T
-
T
0
-
∫
T
0
T
lim
V
→
0
+
⁡
Φ
T
0
,
V
Φ
ρ
,
V
d
ρ
+
η
1
T
∗
+
η
2
C
∗
+
η
3
V
+
η
4
Z
+
∫
0
τ
1
Φ
T
t
-
θ
,
V
t
-
θ
d
θ
+
η
6
∫
0
τ
2
T
∗
t
-
θ
d
θ
+
η
7
∫
0
τ
3
C
∗
t
-
θ
d
θ
.Then,(81)
d
U
0
d
t
=
1
-
lim
V
→
0
+
⁡
Φ
T
0
,
V
Φ
T
,
V
π
-
d
T
-
Φ
T
,
V
+
η
1
e
-
m
1
τ
1
1
-
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
b
T
∗
+
η
2
e
-
m
1
τ
1
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
a
C
∗
+
η
3
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
g
V
Z
-
μ
Z
+
Φ
T
,
V
-
Φ
T
t
-
τ
1
,
V
t
-
τ
1
+
η
6
T
∗
-
η
6
T
∗
t
-
τ
2
+
η
7
C
∗
-
η
7
C
∗
t
-
τ
3
.Equation (81) can be simplified as(82)
d
U
0
d
t
=
1
-
lim
V
→
0
+
⁡
Φ
T
0
,
V
Φ
T
,
V
π
-
d
T
+
Φ
T
,
V
lim
V
→
0
+
⁡
Φ
T
0
,
V
Φ
T
,
V
-
c
η
3
V
-
η
4
μ
Z
=
π
1
-
∂
Φ
T
0
,
0
/
∂
V
∂
Φ
T
,
0
/
∂
V
1
-
T
T
0
+
Φ
T
,
V
∂
Φ
T
0
,
0
/
∂
V
∂
Φ
T
,
0
/
∂
V
-
c
η
3
V
-
η
4
μ
Z
⩽
π
1
-
∂
Φ
T
0
,
0
/
∂
V
∂
Φ
T
,
0
/
∂
V
1
-
T
T
0
+
V
∂
Φ
T
0
,
0
∂
V
-
c
η
3
V
-
η
4
μ
Z
=
π
1
-
∂
Φ
T
0
,
0
/
∂
V
∂
Φ
T
,
0
/
∂
V
1
-
T
T
0
+
c
η
3
1
c
η
3
∂
Φ
T
0
,
0
∂
V
-
1
V
-
η
4
μ
Z
=
π
1
-
∂
Φ
T
0
,
0
/
∂
V
∂
Φ
T
,
0
/
∂
V
1
-
T
T
0
+
c
η
3
R
0
-
1
V
-
η
4
μ
Z
.Based on Condition C2, the first term of (82) is less than or equal to zero. Therefore, if R
0
≤
1, then d
U
0
/
d
t
≤
0, for all T
,
V
,
Z
>
0. Similar to the previous sections, one can show that S
0 is GAS.Condition C3. Consider(83)
Φ
T
,
V
Φ
T
,
V
i
-
V
V
i
1
-
Φ
T
,
V
i
Φ
T
,
V
⩽
0
,
∀
T
,
V
∈
0
,
∞
,
i
=
1,2
.Lemma 12.
Suppose thatR
0
>
1 and Conditions C1–C3 hold. Then, T
1
,
T
2
,
T
1
∗
,
T
2
∗
,
C
1
∗
,
C
2
∗
,
V
1
,
V
2 exist satisfying(84)
sgn
⁡
R
1
-
1
=
sgn
⁡
T
2
-
T
1
=
sgn
⁡
V
1
-
V
2
.Proof.
From Condition C1, forT
1
,
T
2
,
V
1
,
V
2
>
0, we have(85)
Φ
T
2
,
V
2
-
Φ
T
1
,
V
2
T
2
-
T
1
>
0
.
(86)
Φ
T
1
,
V
2
-
Φ
T
1
,
V
1
V
2
-
V
1
>
0
,Applying Condition C3 when i
=
1,  T
=
T
1, and V
=
V
2, we get(87)
Φ
T
1
,
V
2
V
1
-
Φ
T
1
,
V
1
V
2
Φ
T
1
,
V
2
-
Φ
T
1
,
V
1
⩽
0
.It follows from inequality (86) that (88)
Φ
T
1
,
V
2
V
1
-
Φ
T
1
,
V
1
V
2
V
1
-
V
2
>
0
.We assume that sgn
⁡
V
2
-
V
1
=
sgn
⁡
T
2
-
T
1. Then, we have(89)
π
-
d
T
2
-
π
-
d
T
1
=
Φ
T
2
,
V
2
-
Φ
T
1
,
V
1
=
Φ
T
2
,
V
2
-
Φ
T
1
,
V
2
+
Φ
T
1
,
V
2
-
Φ
T
1
,
V
1
,and using inequalities (85) and (86), we obtain sgn
⁡
T
1
-
T
2
=
sgn
⁡
T
2
-
T
1, and this leads to contradiction. Therefore, sgn
⁡
V
1
-
V
2
=
sgn
⁡
T
2
-
T
1. From the conditions of S
1, we have γ
Φ
T
1
,
V
1
/
c
V
1
=
1; then,(90)
R
1
-
1
=
γ
c
1
V
2
Φ
T
2
,
V
2
-
Φ
T
1
,
V
2
+
1
V
1
V
2
Φ
T
1
,
V
2
V
1
-
Φ
T
1
,
V
1
V
2
.Using inequalities (85) and (88), we obtain sgn
⁡
R
1
-
1
=
sgn
⁡
V
1
-
V
2.Theorem 13.
Assume thatR
1
≤
1
<
R
0 and Conditions C1–C3 are satisfied; then, S
1 for system (57)–(61) is GAS.Proof.
Let us define(91)
U
1
=
T
-
T
1
-
∫
T
1
T
Φ
T
1
,
V
1
Φ
ρ
,
V
1
d
ρ
+
η
1
T
1
∗
H
T
∗
T
1
∗
+
η
2
C
1
∗
H
C
∗
C
1
∗
+
η
3
V
1
H
V
V
1
+
η
4
Z
+
Φ
T
1
,
V
1
∫
0
τ
1
H
Φ
T
t
-
θ
,
V
t
-
θ
Φ
T
1
,
V
1
d
θ
+
η
6
T
1
∗
∫
0
τ
2
H
T
∗
t
-
θ
T
1
∗
d
θ
+
η
7
C
1
∗
∫
0
τ
3
H
C
∗
t
-
θ
C
1
∗
d
θ
.Then,(92)
d
U
1
d
t
=
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
π
-
d
T
-
Φ
T
,
V
+
η
1
1
-
T
1
∗
T
∗
e
-
m
1
τ
1
1
-
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
b
T
∗
+
η
2
1
-
C
1
∗
C
∗
e
-
m
1
τ
1
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
a
C
∗
+
η
3
1
-
V
1
V
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
Z
-
c
V
+
η
4
g
V
Z
-
μ
Z
+
Φ
T
,
V
-
Φ
T
t
-
τ
1
,
V
t
-
τ
1
+
Φ
T
1
,
V
1
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
6
T
∗
-
T
∗
t
-
τ
2
+
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
7
C
∗
-
C
∗
t
-
τ
3
+
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Collecting terms of (92), we obtain (93)
d
U
1
d
t
=
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
π
-
d
T
+
Φ
T
,
V
Φ
T
1
,
V
1
Φ
T
,
V
1
-
e
-
m
1
τ
1
η
1
1
-
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
1
∗
T
∗
+
η
1
b
T
1
∗
-
e
-
m
1
τ
1
η
2
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
1
∗
C
∗
+
η
2
a
C
1
∗
-
b
η
1
T
∗
t
-
τ
2
V
1
V
-
a
η
2
C
∗
t
-
τ
3
V
1
V
-
c
η
3
V
+
c
η
3
V
1
+
η
3
r
V
1
Z
-
η
4
μ
Z
+
Φ
T
1
,
V
1
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Applying π
=
d
T
1
+
Φ
(
T
1
,
V
1
), we get(94)
d
U
1
d
t
=
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
d
T
1
-
d
T
+
Φ
T
1
,
V
1
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
+
Φ
T
,
V
Φ
T
1
,
V
1
Φ
T
,
V
1
-
e
-
m
1
τ
1
η
1
1
-
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
1
∗
T
∗
+
η
1
b
T
1
∗
-
e
-
m
1
τ
1
η
2
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
1
∗
C
∗
+
η
2
a
C
1
∗
-
η
1
b
T
∗
t
-
τ
2
V
1
V
-
η
2
a
C
∗
t
-
τ
3
V
1
V
-
c
η
3
V
+
c
η
3
V
1
+
η
3
r
V
1
Z
-
η
4
μ
Z
+
Φ
T
1
,
V
1
e
-
m
1
τ
1
η
1
1
-
ξ
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
Φ
T
1
,
V
1
e
-
m
1
τ
1
η
2
ξ
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the steady state conditions for S
1,(95)
1
-
ξ
Φ
T
1
,
V
1
=
e
m
1
τ
1
b
T
1
∗
,
ξ
Φ
T
1
,
V
1
=
e
m
1
τ
1
a
C
1
∗
,
c
V
1
=
e
-
m
2
τ
2
N
T
b
T
1
∗
+
e
-
m
3
τ
3
N
C
a
C
1
∗
,then we have Φ
(
T
1
,
V
1
)
=
η
1
b
T
1
∗
+
η
2
a
C
1
∗ and (96)
d
U
1
d
t
=
d
T
1
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
1
-
T
T
1
+
Φ
T
1
,
V
1
Φ
T
,
V
Φ
T
,
V
1
+
η
1
b
T
1
∗
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
+
η
2
a
C
1
∗
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
-
η
1
b
T
1
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
1
∗
Φ
T
1
,
V
1
T
∗
+
η
1
b
T
1
∗
-
η
2
a
C
1
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
1
∗
Φ
T
1
,
V
1
C
∗
+
η
2
a
C
1
∗
-
η
1
b
T
1
∗
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
η
2
a
C
1
∗
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
η
1
b
T
1
∗
+
η
2
a
C
1
∗
-
Φ
T
1
,
V
1
V
V
1
+
η
3
r
V
1
-
μ
g
Z
+
η
1
b
T
1
∗
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
2
a
C
1
∗
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
1
b
T
1
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
1
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the following equalities,(97)
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
=
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
j
∗
Φ
T
j
,
V
j
T
∗
+
ln
⁡
Φ
T
j
,
V
j
Φ
T
,
V
j
+
ln
⁡
Φ
T
,
V
j
V
Φ
T
,
V
V
j
+
ln
⁡
V
j
T
∗
V
T
j
∗
,
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
=
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
j
∗
Φ
T
j
,
V
j
C
∗
+
ln
⁡
Φ
T
j
,
V
j
Φ
T
,
V
j
ln
⁡
Φ
T
,
V
j
V
Φ
T
,
V
V
j
+
ln
⁡
V
j
C
∗
V
C
j
∗
,
ln
⁡
T
∗
t
-
τ
2
T
∗
=
ln
⁡
T
∗
t
-
τ
2
V
j
V
T
j
∗
+
ln
⁡
V
T
j
∗
V
j
T
∗
,
ln
⁡
C
∗
t
-
τ
3
C
∗
=
ln
⁡
C
∗
t
-
τ
3
V
j
V
C
j
∗
+
ln
⁡
V
C
j
∗
V
j
C
∗
,
j
=
1,2
,in case of j
=
1, we obtain(98)
d
U
1
d
t
=
d
T
1
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
1
-
T
T
1
+
Φ
T
1
,
V
1
Φ
T
,
V
Φ
T
,
V
1
-
V
V
1
-
1
+
Φ
T
,
V
1
V
Φ
T
,
V
V
1
-
η
1
b
T
1
∗
Φ
T
1
,
V
1
Φ
T
,
V
1
-
1
-
ln
⁡
Φ
T
1
,
V
1
Φ
T
,
V
1
-
η
2
a
C
1
∗
Φ
T
1
,
V
1
Φ
T
,
V
1
-
1
-
ln
⁡
Φ
T
1
,
V
1
Φ
T
,
V
1
-
η
1
b
T
1
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
1
∗
Φ
T
1
,
V
1
T
∗
-
1
-
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
1
∗
Φ
T
1
,
V
1
T
∗
-
η
2
a
C
1
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
1
∗
Φ
T
1
,
V
1
C
∗
-
1
-
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
1
∗
Φ
T
1
,
V
1
C
∗
-
η
1
b
T
1
∗
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
1
-
ln
⁡
T
∗
t
-
τ
2
V
1
V
T
1
∗
-
η
2
a
C
1
∗
C
∗
t
-
τ
3
V
1
V
C
1
∗
-
1
-
ln
⁡
C
∗
t
-
τ
3
V
1
V
C
1
∗
-
η
1
b
T
1
∗
Φ
T
,
V
1
V
Φ
T
,
V
V
1
-
1
-
ln
⁡
Φ
T
,
V
1
V
Φ
T
,
V
V
1
-
η
2
a
C
1
∗
Φ
T
,
V
1
V
Φ
T
,
V
V
1
-
1
-
ln
⁡
Φ
T
,
V
1
V
Φ
T
,
V
V
1
+
r
η
3
V
1
-
μ
g
Z
=
d
T
1
1
-
Φ
T
1
,
V
1
Φ
T
,
V
1
1
-
T
T
1
+
Φ
T
1
,
V
1
Φ
T
,
V
Φ
T
,
V
1
-
V
V
1
1
-
Φ
T
,
V
1
Φ
T
,
V
-
η
1
b
T
1
∗
H
Φ
T
1
,
V
1
Φ
T
,
V
1
+
H
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
1
∗
Φ
T
1
,
V
1
T
∗
+
H
T
∗
t
-
τ
2
V
1
V
T
1
∗
+
H
Φ
T
,
V
1
V
Φ
T
,
V
V
1
-
η
2
a
C
1
∗
H
Φ
T
1
,
V
1
Φ
T
,
V
1
+
H
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
1
∗
Φ
T
1
,
V
1
C
∗
+
H
C
∗
t
-
τ
3
V
1
V
C
1
∗
+
H
Φ
T
,
V
1
V
Φ
T
,
V
V
1
+
r
η
3
V
1
-
μ
g
Z
.Since R
1
≤
1, then Conditions C1 and C3 and Lemma 12 imply that d
U
1
/
d
t
≤
0, for all T
,
T
∗
,
C
∗
,
V
,
Z
>
0, where the equality occurs at S
1. LIP implies that S
1 is GAS.Theorem 14.
LetR
1
>
1 and Conditions C1–C3 be satisfied; then, S
2 for system (57)–(61) is GAS.Proof.
We consider a Lyapunov functional:(99)
U
2
=
T
-
T
2
-
∫
T
2
T
Φ
T
2
,
V
2
Φ
ρ
,
V
2
d
ρ
+
η
1
T
2
∗
H
T
∗
T
2
∗
+
η
2
C
2
∗
H
C
∗
C
2
∗
+
η
3
V
2
H
V
V
2
+
η
4
Z
2
H
Z
Z
2
+
Φ
T
2
,
V
2
∫
0
τ
1
H
Φ
T
t
-
θ
,
V
t
-
θ
Φ
T
2
,
V
2
d
θ
+
η
6
T
2
∗
∫
0
τ
2
H
T
∗
t
-
θ
T
2
∗
d
θ
+
η
7
C
2
∗
∫
0
τ
3
H
C
∗
t
-
θ
C
2
∗
d
θ
.Then,(100)
d
U
2
d
t
=
1
-
Φ
T
2
,
V
2
Φ
T
,
V
2
π
-
d
T
-
Φ
T
,
V
+
η
1
1
-
T
2
∗
T
∗
e
-
m
1
τ
1
1
-
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
b
T
∗
+
η
2
1
-
C
2
∗
C
∗
e
-
m
1
τ
1
ξ
Φ
T
t
-
τ
1
,
V
t
-
τ
1
-
a
C
∗
+
η
3
1
-
V
2
V
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
c
V
-
r
V
Z
+
η
4
1
-
Z
2
Z
g
V
Z
-
μ
Z
+
Φ
T
,
V
-
Φ
T
t
-
τ
1
,
V
t
-
τ
1
+
Φ
T
2
,
V
2
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
6
T
∗
-
T
∗
t
-
τ
2
+
T
2
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
7
C
∗
-
C
∗
t
-
τ
3
+
C
2
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using the steady state conditions for S
2,(101)
π
=
d
T
2
+
Φ
T
2
,
V
2
,
1
-
ξ
Φ
T
2
,
V
2
=
e
m
1
τ
1
b
T
2
∗
,
ξ
Φ
T
2
,
V
2
=
e
m
1
τ
1
a
C
2
∗
,
c
V
2
=
e
-
m
2
τ
2
N
T
b
T
2
∗
+
e
-
m
3
τ
3
N
C
a
C
2
∗
-
r
V
2
Z
2
, we get Φ
(
T
2
,
V
2
)
=
η
1
b
T
2
∗
+
η
2
a
C
2
∗, and (102)
d
U
2
d
t
=
d
T
2
1
-
Φ
T
2
,
V
2
Φ
T
,
V
2
1
-
T
T
2
+
Φ
T
2
,
V
2
Φ
T
,
V
Φ
T
,
V
2
+
η
1
b
T
2
∗
1
-
Φ
T
2
,
V
2
Φ
T
,
V
2
+
η
2
a
C
2
∗
1
-
Φ
T
2
,
V
2
Φ
T
,
V
2
-
η
1
b
T
2
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
2
∗
Φ
T
2
,
V
2
T
∗
+
η
1
b
T
2
∗
-
η
2
a
C
2
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
2
∗
Φ
T
2
,
V
2
C
∗
+
η
2
a
C
2
∗
-
η
1
b
T
2
∗
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
η
2
a
C
2
∗
C
∗
t
-
τ
3
V
2
V
C
2
∗
+
η
1
b
T
2
∗
+
η
2
a
C
2
∗
-
Φ
T
2
,
V
2
V
V
2
+
η
1
b
T
2
∗
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
2
a
C
2
∗
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
Φ
T
,
V
+
η
1
b
T
2
∗
ln
⁡
T
∗
t
-
τ
2
T
∗
+
η
2
a
C
2
∗
ln
⁡
C
∗
t
-
τ
3
C
∗
.Using (97), when j
=
2, we get(103)
d
U
2
d
t
=
d
T
2
1
-
Φ
T
2
,
V
2
Φ
T
,
V
2
1
-
T
T
2
+
Φ
T
2
,
V
2
Φ
T
,
V
Φ
T
,
V
2
-
V
V
2
-
1
+
Φ
T
,
V
2
V
Φ
T
,
V
V
2
-
η
1
b
T
2
∗
Φ
T
2
,
V
2
Φ
T
,
V
2
-
1
-
ln
⁡
Φ
T
2
,
V
2
Φ
T
,
V
2
-
η
2
a
C
2
∗
Φ
T
2
,
V
2
Φ
T
,
V
2
-
1
-
ln
⁡
Φ
T
2
,
V
2
Φ
T
,
V
2
-
η
1
b
T
2
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
2
∗
Φ
T
2
,
V
2
T
∗
-
1
-
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
2
∗
Φ
T
2
,
V
2
T
∗
-
η
2
a
C
2
∗
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
2
∗
Φ
T
2
,
V
2
C
∗
-
1
-
ln
⁡
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
2
∗
Φ
T
2
,
V
2
C
∗
-
η
1
b
T
2
∗
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
1
-
ln
⁡
T
∗
t
-
τ
2
V
2
V
T
2
∗
-
η
2
a
C
2
∗
C
∗
t
-
τ
3
V
2
V
C
2
∗
-
1
-
ln
⁡
C
∗
t
-
τ
3
V
2
V
C
2
∗
-
η
1
b
T
2
∗
Φ
T
,
V
2
V
Φ
T
,
V
V
2
-
1
-
ln
⁡
Φ
T
,
V
2
V
Φ
T
,
V
V
2
-
η
2
a
C
2
∗
Φ
T
,
V
2
V
Φ
T
,
V
V
2
-
1
-
ln
⁡
Φ
T
,
V
2
V
Φ
T
,
V
V
2
=
d
T
2
1
-
Φ
T
2
,
V
2
Φ
T
,
V
2
1
-
T
T
2
+
Φ
T
2
,
V
2
Φ
T
,
V
Φ
T
,
V
2
-
V
V
2
1
-
Φ
T
,
V
2
Φ
T
,
V
-
η
1
b
T
2
∗
H
Φ
T
2
,
V
2
Φ
T
,
V
2
+
H
Φ
T
t
-
τ
1
,
V
t
-
τ
1
T
2
∗
Φ
T
2
,
V
2
T
∗
+
H
T
∗
t
-
τ
2
V
2
V
T
2
∗
+
H
Φ
T
,
V
2
V
Φ
T
,
V
V
2
-
η
2
a
C
2
∗
H
Φ
T
1
,
V
1
Φ
T
,
V
1
+
H
Φ
T
t
-
τ
1
,
V
t
-
τ
1
C
2
∗
Φ
T
2
,
V
2
C
∗
+
H
C
∗
t
-
τ
3
V
2
V
C
2
∗
+
H
Φ
T
,
V
2
V
Φ
T
,
V
V
2
.It is easy to show that S
2 is GAS.


## 5. Numerical Simulations

In this section, we show an example of the general model (50)-(57) where Conditions C1–C3 can be satisfied:(104)
T
˙
t
=
π
-
d
T
t
-
1
-
ε
k
T
t
V
t
1
+
ω
T
t
+
β
V
t
,
T
˙
∗
t
=
e
-
m
1
τ
1
1
-
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
ω
T
t
-
τ
1
+
β
V
t
-
τ
1
-
b
T
∗
t
,
C
˙
∗
t
=
e
-
m
1
τ
1
ξ
1
-
ε
k
T
t
-
τ
1
V
t
-
τ
1
1
+
ω
T
t
-
τ
1
+
β
V
t
-
τ
1
-
a
C
∗
t
,
V
˙
t
=
e
-
m
2
τ
2
N
T
b
T
∗
t
-
τ
2
+
e
-
m
3
τ
3
N
C
a
C
∗
t
-
τ
3
-
r
V
t
Z
t
-
c
V
t
,
Z
˙
t
=
g
V
t
Z
t
-
μ
Z
t
,where ω,  β
⩾
0. In this example, the incidence rate is given by(105)
Φ
T
,
V
=
1
-
ε
k
T
V
1
+
ω
T
+
β
V
. We have the following:(106)
Φ
T
,
V
>
0
,
Φ
T
,
0
=
Φ
0
,
V
=
0
,
∀
T
>
0
,
V
>
0
,
∂
Φ
T
,
V
∂
T
=
1
-
ε
k
V
1
+
β
V
1
+
ω
T
+
β
V
2
>
0
,
∀
T
>
0
,
V
>
0
,
∂
Φ
T
,
V
∂
V
=
1
-
ε
k
T
1
+
ω
T
1
+
ω
T
+
β
V
2
>
0
,
∀
T
>
0
,
V
>
0
,
∂
Φ
T
,
0
∂
V
=
1
-
ε
k
T
1
+
ω
T
>
0
,
∀
T
>
0
,
Φ
T
,
V
=
1
-
ε
k
T
V
1
+
ω
T
+
β
V
⩽
1
-
ε
k
T
V
1
+
ω
T
=
V
∂
Φ
T
,
0
∂
V
,
∀
T
>
0
,
V
>
0
,
Φ
T
,
V
Φ
T
,
V
i
-
V
V
i
1
-
Φ
T
,
V
i
Φ
T
,
V
=
-
β
1
+
ω
T
V
-
V
i
2
V
i
1
+
ω
T
+
β
V
i
1
+
ω
T
+
β
V
⩽
0
,
∀
T
,
V
>
0
,
i
=
1,2
.Then, Conditions C1–C3 are satisfied. The parameters R
0 and R
1 will be(107)
R
0
=
1
-
ε
k
T
0
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
c
1
+
ω
T
0
,
R
1
=
1
-
ε
k
T
2
e
-
m
1
τ
1
+
m
2
τ
2
1
-
ξ
N
T
+
e
-
m
1
τ
1
+
m
3
τ
3
ξ
N
C
c
1
+
ω
T
2
+
β
V
2
,where(108)
T
2
=
1
2
ω
ω
T
0
-
1
+
ζ
V
2
+
1
+
ζ
V
2
-
ω
T
0
2
+
4
ω
T
0
1
+
β
V
2
,
V
2
=
μ
g
,
ζ
=
β
+
1
-
ε
k
d
.The values of some parameters of model (104) are listed in Table 1. The remaining parameters of the model will be chosen varied. Without loss of generality, we let τ
e
=
τ
1
=
τ
2
=
τ
3. We have used MATLAB for the numerical computations. Now, we study the following cases.Table 1
The parameters of model (104).
Parameter
Value
π
10 cells mm−3 day−1
d
0.01 day−1
b
0.3 day−1
ω
0.05 cells−1 mm3
β
0.05 virus−1 mm3
a
0.25 day−1
N
T
100 virus cells−1
ξ
0.5
N
C
50 virus cells−1
r
0.01 cells−1 mm3 day−1
μ
0.05 day−1
m
1
1 day−1
m
2
1 day−1
m
3
1 day−1
c
3 day−1Case 1 (effect of the parametersk and g on the stability of the system).
Now, we confirm the results of Theorems11–14. The evolution of the dynamics of model (104) was observed over a time interval [
0,1000
]. We have chosen the initial conditions as follows: T
(
0
)
=
600, T
∗
(
0
)
=
10, C
∗
(
0
)
=
5
,
V
(
0
)
=
10, and Z
(
0
)
=
80. We use three sets for the parameters k and g to get the following three subcases. We fix the values τ
e
=
0.5 and ε
=
0:(i)
R
0
≤
1. The values of k and g are chosen as k
=
0.002 virus−1 mm3 day−1 and g
=
0.0005 virus−1 mm3 day−1. The values of the two bifurcation parameters are computed as R
0
=
0.3607
<
1 and R
1
=
0.3114
<
1. This means that S
0 is GAS. One can see from Figures 1–5 that the numerical results confirm the results of Theorem 11. We observe that the states of the system eventually approach the steady state S
0
(
1000,0
,
0,0
,
0
). In this case, the HIV will be removed from the blood.
(ii)
R
1
≤
1
<
R
0. We take the values k
=
0.01 virus−1 mm3 day−1 and g
=
0.0005 virus−1 mm3 day−1. In this case, R
1
=
0.8291
<
1
<
R
0
=
1.8033. According to Theorem 13, S
1 is GAS. From Figures 1–5, we can see that there is a consistency between the numerical and theoretical results of Theorem 13. Moreover, the states of the system converge to the steady S
1
120.22,8.89,10.67,80.91,0.
(iii)
R
1
>
1. We choose k
=
0.01 virus−1 mm3 day−1 and g
=
0.002 virus−1 mm3 day−1. Then, we compute R
0
=
1.8033
>
1 and R
1
=
1.6976
>
1. From Figures 1–5, we can see that the states of the system approach the infected steady state with humoral immune response S
2
538.56,4.66,5.60,25.00,209.27. This supports the results of Theorem 14 that the infected steady state without humoral immune response S
2 is GAS.Figure 1
The concentration of uninfected CD4+ T cells for model (104).Figure 2
The concentration of short-lived infected cells for model (104).Figure 3
The concentration of long-lived chronically infected cells for model (104).Figure 4
The concentration of free virus particles for model (104).Figure 5
The concentration of B cells for model (104).Case 2 (effect of the parameterε on the stability of the system).
For this case, we takeτ
e
=
0.5 day−1, k
=
0.01 virus−1 mm3 day−1, and g
=
0.001 virus−1 mm3 day−1. In Figures 6–10, we show the effect of the drug efficacy ε on the stability of the steady states of the system. We can see from Figures 6–10 that as the drug efficacy ε is increased, the concentration of uninfected CD4+ T cells is increased, while the concentrations of free virus particles, B cells, short-lived infected cells, and long-lived chronically infected cells are decreased. From these figures and Table 2, we can see that the values of R
0 and R
1 are decreased as ε is increased. Using the values of the parameters given in Table 1, we obtain the following:(i)
If0
≤
ε
<
0.3730, then S
2 exists and it is GAS.
(ii)
If0.3730
≤
ε
<
0.4455, then S
1 exists and it is GAS.
(iii)
If0.4455
≤
ε
<
1, then S
0 is GAS.
It means that the numerical results and the results of Theorems11–14 are compatible. In this case, the treatment with sufficient drug efficacy can succeed to clear the HIV from the plasma.Table 2
The values of equilibria,R
0, and R
1 for model (104) with different values of ε.
Equilibria
R
0
R
1
ε
=
0.0
E
2  (231.88, 7.77, 9.32, 50, 123.86)
1.8033
1.4129
ε
=
0.2
E
2 (337.44, 6.70, 8.04, 50, 65.61)
1.4427
1.2187
ε
=
0.3730
E
1  (456.34, 5.50, 6.60, 50, 0)
1.1308
1
ε
=
0.4
E
1 (572.42, 4.32, 5.19, 39.32, 0)
1.0820
0.9624
ε
=
0.4455
E
0 (1000, 0, 0, 0, 0)
1
0.8974
ε
=
0.6
E
0 (1000, 0, 0, 0, 0)
0.7213
0.6632
ε
=
0.8
E
0 (1000, 0, 0, 0, 0)
0.3607
0.3388Figure 6
The concentration of uninfected CD4+ T cells for model (104).Figure 7
The concentration of short-lived infected cells for model (104).Figure 8
The concentration of long-lived chronically infected cells for model (104).Figure 9
The concentration of free virus particles for model (104).Figure 10
The concentration of B cells for model (104).Case 3 (effect of the parameterτ
e on the stability of the system).
We fix the parametersε
=
0.4, k
=
0.01 virus−1 mm3 day−1, and g
=
0.001 virus−1 mm m3 day−1. In Figures 11–15, we show the effect of the time delay parameter τ
e on the stability of the steady states of the system. We can see from Figures 11–15 that the time delay parameter τ
e plays a similar role as the drug efficacy parameter ε. From these figures and Table 3, we can see that the values of R
0 and R
1 are decreased when τ
e is increased. Using the values of the parameters given in Table 1, we obtain the following:(i)
If0
≤
τ
e
<
0.4808, then S
2 exists and it is GAS.
(ii)
If0.4808
≤
τ
e
<
0.5394, then S
1 exists and it is GAS.
(iii)
If0.5394
≤
τ
e, then S
0 is GAS.
Figures11–15 show that the numerical results are also compatible with the results of Theorems 11–14.Table 3
The values of equilibria,R
0, and R
1 for model (104) with different values of τ
e.
Delay
Equilibria
R
0
R
1
τ
e
=
0.0
E
2 (476.81, 8.72, 10.46, 50., 484.77)
2.9412
2.6160
τ
e
=
0.1
E
2 (476.81, 7.89, 9.47, 50, 342.53)
2.4080
2.1418
τ
e
=
0.3
E
2 (476.81, 6.46, 7.75, 50, 130.70)
1.6142
1.4357
τ
e
=
0.4808
E
1 (476.81, 5.39, 6.47, 50, 0)
1.1243
1.0
τ
e
=
0.5
E
1 (572.42, 4.32, 5.18, 39.32, 0)
1.0820
0.9624
τ
e
=
0.5394
E
0 (1000, 0, 0, 0, 0)
1.0
0.8894
τ
e
=
0.7
E
0 (1000, 0, 0, 0, 0)
0.7253
0.6451
τ
e
=
1.0
E
0 (1000, 0, 0, 0, 0)
0.3981
0.3540
τ
e
=
2.0
E
0 (1000, 0, 0, 0, 0)
0.0539
0.0479Figure 11
The concentration of uninfected CD4+ T cells for model (104).Figure 12
The concentration of short-lived infected cells for model (104).Figure 13
The concentration of long-lived chronically infected cells for model (104).Figure 14
The concentration of free virus particles for model (104).Figure 15
The concentration of B cells for model (104).


## 6. Conclusions

In this paper, we have proposed three HIV infection models with humoral immune response and two types of infected cells. The models incorporate two types of discrete delays. The incidence rate of infection is given by bilinear, saturated functional response and general nonlinear function in the first, second, and third model, respectively. We have derived two bifurcation parameters,R
0 and R
1. The global stability of all steady states of the models has been established using Lyapunov method. We have presented an example for the general incidence rate and performed numerical simulations to support our theoretical results.



---

*Source: 370968-2015-12-06.xml*