# Stability of Virus Infection Models with Antibodies and Chronically Infected Cells

Authors: Mustafa A. Obaid; A. M. Elaiw
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Two virus infection models with antibody immune response and chronically infected cells are proposed and analyzed. Bilinear incidence rate is considered in the first model, while the incidence rate is given by a saturated functional response in the second one. One main feature of these models is that it includes both short-lived infected cells and chronically infected cells. The chronically infected cells produce much smaller amounts of virus than the short-lived infected cells and die at a much slower rate. Our mathematical analysis establishes that the global dynamics of the two models are determined by two threshold parametersR
0 and R
1. By constructing Lyapunov functions and using LaSalle's invariance principle, we have established the global asymptotic stability of all steady states of the models. We have proven that, the uninfected steady state is globally asymptotically
stable (GAS) if R
0
<
1, the infected steady state without antibody immune response exists and it
is GAS if R
1
<
1
<
R
0, and the infected steady state with antibody immune response exists and it is GAS if R
1
>
1. We check our theorems with numerical simulation in the end.


---

## Body

## 1. Introduction

In recent years, many mathematical models have been proposed to study the dynamics of viral infections such as the human immunodeficiency virus (HIV), the hepatitis C virus (HCV), and the hepatitis B virus (HBV) (see, e.g., [1–17]). Such virus infection models can be very useful in the control of epidemic diseases and provide insights into the dynamics of viral load in vivo. Therefore, mathematical analysis of the virus infection models can play a significant role in the development of a better understanding of diseases and various drug therapy strategies. Most of the mathematical models of viral infection presented in the literature did not differentiate between the short-lived infected cells and chronically infected cells. The chronically infected cells produce much smaller amounts of virus than the short-lived infected cells and die at a much slower rate [18]. The virus dynamics model with chronically infected cells and under the effect of antiviral drug therapy was introduced in [18] as
(1)
T
˙
=
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
,
T
˙
*
=
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
-
δ
T
*
,
C
˙
*
=
α
(
1
-
ɛ
)
k
T
V
-
a
C
*
,
V
˙
=
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
,
where T, T
*, C
*, and V are the concentration of the uninfected cells, short-lived infected cells, chronically infected cells, and free virus particles, respectively. The constant λ is the rate at which new uninfected cells are generated and d is the natural death rate constant of uninfected cells. k is the infection rate constant. The fractions (
1
-
α
) and α with 0
<
α
<
1 are the probabilities that, upon infection, an uninfected cell will become either short-lived infected or chronically infected. δ and a are the death rate constants of the short-lived infected cells and chronically infected cells, respectively. N
T and N
C are the average number of virions produced in the lifetime of the short-lived infected and chronically infected cells, respectively. The chronically infected cells produce much smaller amounts of virus than the short-lived infected cells and die at a much slower rate (i.e., N
T
>
N
C and δ
>
a). The free viruses are cleared with rate constant c. The drug efficacy is denoted by ɛ and 0
≤
ɛ
≤
1.It is observed that the basic and global properties of model (1) are not studied in the literature. Moreover, model (1) did not take into consideration the immune response. During viral infections, the host immune system reacts with antigen-specific immune response. The immune system is described as having two “arms”: the cellular arm, which depends on T cells to mediate attacks on virally infected or cancerous cells, and the humoral arm, which depends on B cells. The B cell is a type of blood cell which belongs to a group of white blood cells (WBCs) called lymphocytes. WBCs protect the body from infection. The main job of B cells is to fight infection. B cells get activated when an infection occurs and they produce molecules called antibodies that attach to the surface of the infectious agent. These antibodies either kill the infection causing organism or make it prone to attack by other WBCs. They play a major role in the immune system, which guards the body against infection. Virus infection models with antibody immune response have been analyzed by many researchers (see [19–28]). However, in all of these works, the chronically infected cells have been neglected.In this paper, we propose two virus infection models with antibody immune response and chronically infected cells. In the first model, bilinear incidence rate which is based on the law of mass-action is considered. The second model generalizes the first one where the incidence rate is given by a saturation functional response. The global stability of all equilibria of the models is established using the method of Lyapunov function. We prove that the global dynamics of the models are determined by two threshold parametersR
0 and R
1. If R
0
≤
1, then the infection-free equilibrium is globally asymptotically stable (GAS), if R
1
≤
1
<
R
0, then the infected equilibrium without antibody immune response exists and it is GAS, and if R
1
>
1 then the infected equilibrium with antibody immune response exists and it is GAS.


## 2. Model with Bilinear Incidence Rate

In this section we propose a viral dynamics model with antibody immune response, taking into consideration the chronically infected cells. Based on the mass-action principle, we assume that the incidence rate of infection is bilinear; that is, the infection rate per virus and per uninfected cell is constant:(2)
T
˙
=
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
,
(3)
T
˙
*
=
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
-
δ
T
*
,
(4)
C
˙
*
=
α
(
1
-
ɛ
)
k
T
V
-
a
C
*
,
(5)
V
˙
=
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
,
(6)
Z
˙
=
g
V
Z
-
μ
Z
,
where Z is the concentration of antibody immune cells. The viruses are attacked by the antibodies with rate r
V
Z. The antibody immune cells are proliferated at rate g
V
Z and die at rate μ
Z. All the other variables and parameters of the model have the same meanings as given in (1).


### 2.1. Positive Invariance

We note that model (2)–(6) are biologically acceptable in the sense that no population goes negative. It is straightforward to check the positive invariance of the nonnegative orthant ℝ
+
5 by model (2)–(6) (see, e.g., [6]). In the following, we show the boundedness of the solution of model (2)–(6).Proposition 1.
There exist positive numbersL
i, i
=
1,2
,
3, such that the compact set
(7)
Ω
=
{
(
T
,
T
*
,
C
*
,
V
,
Z
)
∈
ℝ
+
4
:
0
≤
T
,
T
*
,
C
*
≤
L
1
,
r
0
≤
V
≤
L
2
,
0
≤
Z
≤
L
3
{
(
T
,
T
*
,
C
*
,
V
,
Z
)
∈
ℝ
+
4
:
0
≤
T
,
T
*
,
C
*
≤
L
1
,
}
is positively invariant.Proof.
To show the boundedness of the solutions we letG
1
(
t
)
=
T
(
t
)
+
T
*
(
t
)
+
C
*
(
t
); then
(8)
G
˙
1
(
t
)
=
λ
-
d
T
(
t
)
-
(
1
-
ɛ
)
k
T
(
t
)
V
(
t
)
+
(
1
-
α
)
(
1
-
ɛ
)
k
T
(
t
)
V
(
t
)
-
δ
T
*
+
α
(
1
-
ɛ
)
k
T
(
t
)
V
(
t
)
-
a
C
*
(
t
)
≤
λ
-
s
1
G
1
(
t
)
,
where s
1
=
min
⁡
⁡
{
d
,
a
,
δ
}. Hence G
1
(
t
)
≤
L
1, if G
1
(
0
)
≤
L
1 where L
1
=
λ
/
s
1. Since T
(
t
)
>
0, T
*
(
t
)
≥
0, and C
*
(
t
)
≥
0, then 0
≤
T
(
t
), T
*
(
t
), C
*
(
t
)
≤
L
1 if 0
≤
T
(
0
)
+
T
*
(
0
)
+
C
*
(
0
)
≤
L
1. Let G
2
(
t
)
=
V
(
t
)
+
(
r
/
g
)
Z
(
t
); then
(9)
G
˙
2
(
t
)
=
N
T
δ
T
*
(
t
)
+
N
C
a
C
*
(
t
)
-
c
V
(
t
)
-
r
μ
g
Z
(
t
)
≤
(
N
T
δ
+
N
C
a
)
L
1
-
s
2
(
V
(
t
)
+
r
g
Z
(
t
)
)
=
(
N
T
δ
+
N
C
a
)
L
1
-
s
2
G
2
(
t
)
,
where s
2
=
min
⁡
⁡
{
c
,
μ
}. Hence G
2
(
t
)
≤
L
2, if G
2
(
0
)
≤
L
2, where L
2
=
(
N
T
δ
+
N
C
a
)
L
1
/
s
2. Since V
(
t
)
≥
0 and Z
(
t
)
≥
0 then 0
≤
V
(
t
)
≤
L
2 and 0
≤
Z
(
t
)
≤
L
3 if 0
≤
V
(
0
)
+
(
r
/
g
)
Z
(
0
)
≤
L
2, where L
3
=
g
L
2
/
r.


### 2.2. Equilibria

System (2)–(6) always admits an infection-free equilibrium E
0
=
(
T
0
,
0,0
,
0,0
), where T
0
=
λ
/
d. In addition to E
0, the system can have an infected equilibrium without antibody immune response E
1
(
T
1
,
T
1
*
,
C
1
*
,
V
1
,
0
) and an infected equilibrium with antibody immune response E
2
(
T
2
,
T
2
*
,
C
2
*
,
V
2
,
Z
2
) where
(10)
T
1
=
c
(
1
-
ɛ
)
k
[
(
1
-
α
)
N
T
+
α
N
C
]
,
T
1
*
=
(
1
-
α
)
λ
{
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
-
c
}
δ
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
,
C
1
*
=
α
λ
{
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
-
c
}
a
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
,
V
1
=
d
{
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
-
c
}
(
1
-
ɛ
)
k
c
,
T
2
=
λ
g
g
d
+
(
1
-
ɛ
)
k
μ
,
T
2
*
=
(
1
-
α
)
(
1
-
ɛ
)
k
λ
μ
δ
(
d
g
+
(
1
-
ɛ
)
k
μ
)
,
C
2
*
=
α
(
1
-
ɛ
)
k
λ
μ
a
(
d
g
+
(
1
-
ɛ
)
k
μ
)
,
V
2
=
μ
g
,
Z
2
=
c
r
(
d
g
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
c
(
d
g
+
(
1
-
ɛ
)
k
μ
)
-
1
)
.
We discuss the local stability of the infection-free equilibrium E
0. At the infection-free equilibrium E
0
(
T
0
,
0,0
,
0,0
), the system has the Jacobian matrix given by
(11)
J
E
0
=
[
-
d
0
0
-
(
1
-
ɛ
)
k
T
0
0
0
-
δ
0
(
1
-
α
)
(
1
-
ɛ
)
k
T
0
0
0
0
-
a
α
(
1
-
ɛ
)
k
T
0
0
0
δ
N
T
a
N
C
-
c
0
0
0
0
0
-
μ
]
.
The characteristic equation of the Jacobian matrix evaluated at E
0 is
(12)
(
s
+
d
)
(
s
+
μ
)
(
s
3
+
a
1
s
2
+
a
2
s
+
a
3
)
=
0
,
where
(13)
a
1
=
a
+
c
+
δ
,
a
2
=
a
c
+
a
δ
+
c
δ
-
(
1
-
α
)
(
1
-
ɛ
)
k
T
0
N
T
δ
-
α
(
1
-
ɛ
)
k
T
0
N
C
a
,
a
3
=
a
c
δ
(
1
-
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
c
)
.
We observe that (12) has two negative eigenvalues s
1
=
-
d and s
2
=
-
μ. By the Routh-Hurwitz criterion, the remaining three eigenvalues of (12) have negative real parts if a
1
>
0, a
3
>
0, and a
1
a
2
-
a
3
>
0. We have a
1
>
0 and if (
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
/
c
<
1, then a
3
>
0 and
(14)
a
1
a
2
-
a
3
=
a
δ
2
+
a
2
δ
+
2
a
c
δ
+
a
(
a
+
c
)
[
c
-
α
(
1
-
ɛ
)
k
T
0
N
C
]
+
δ
(
δ
+
c
)
[
c
-
(
1
-
α
)
(
1
-
ɛ
)
k
T
0
N
T
]
>
0
.
Now we define the basic reproduction number for system (2)–(6) as
(15)
R
0
=
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
c
.
It follows that the equilibria E
1 and E
2 can be written as
(16)
T
1
=
T
0
R
0
,
T
1
*
=
(
1
-
α
)
λ
δ
(
R
0
-
1
)
R
0
,
C
1
*
=
α
λ
a
(
R
0
-
1
)
R
0
,
V
1
=
d
(
1
-
ɛ
)
k
(
R
0
-
1
)
,
T
2
=
λ
g
g
d
+
(
1
-
ɛ
)
k
μ
,
T
2
*
=
(
1
-
α
)
(
1
-
ɛ
)
k
λ
μ
δ
(
d
g
+
(
1
-
ɛ
)
k
μ
)
,
C
2
*
=
α
(
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 that Z
2
>
0 when d
g
R
0
/
(
d
g
+
(
1
-
ɛ
)
k
μ
)
>
1. Now we define another threshold parameter R
1 as
(17)
R
1
=
R
0
1
+
(
(
1
-
ɛ
)
k
μ
/
d
g
)
.
Clearly R
1
<
R
0.From (16) we have the following statements:(i)
ifR
0
≤
1, then there exists only positive equilibrium E
0;
(ii)
ifR
1
≤
1
<
R
0, then there exist two positive equilibria E
0 and E
1;
(iii)
ifR
1
>
1, then there exist three positive equilibria E
0, E
1, and E
2.


### 2.3. Global Stability Analysis

In this section, we study the global stability of all the equilibria of system (2)–(6) employing the method of Lyapunov function.Theorem 2.
For system (2)–(6), if R
0
≤
1, then E
0 is GAS.Proof.
Define a Lyapunov functionU
0 as follows:
(18)
U
0
=
T
0
(
T
T
0
-
1
-
ln
⁡
(
T
T
0
)
)
+
η
1
T
*
+
η
2
C
*
+
η
3
V
+
η
4
Z
,
where η
i,  i
=
1
,
…
,
4, are positive constants to be determined below. Calculating the derivative of U
0 along the solutions of the system (2)–(6) and applying λ
=
T
0
d, we obtain
(19)
d
U
0
d
t
=
(
1
-
T
0
T
)
(
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
)
+
η
1
(
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
-
δ
T
*
)
+
η
2
(
α
(
1
-
ɛ
)
k
T
V
-
a
C
*
)
+
η
3
(
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
)
+
η
4
(
g
V
Z
-
μ
Z
)
.
Let η
i,  i
=
1
,
…
,
4, be chosen such as
(20)
(
1
-
α
)
η
1
+
α
η
2
=
1
,
η
1
-
N
T
η
3
=
0
,
η
2
-
N
C
η
3
=
0
,
g
η
4
-
r
η
3
=
0
.
The solution of (20) is given by
(21)
η
1
=
N
T
(
1
-
α
)
N
T
+
α
N
C
,
η
2
=
N
C
(
1
-
α
)
N
T
+
α
N
C
,
η
3
=
1
(
1
-
α
)
N
T
+
α
N
C
,
η
4
=
r
g
[
(
1
-
α
)
N
T
+
α
N
C
]
.
The values of η
i,  i
=
1
,
…
,
4, given by (21) will be used throughout the paper. Then
(22)
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
+
η
3
c
(
R
0
-
1
)
V
-
η
4
μ
Z
.
If R
0
≤
1 then d
U
0
/
d
t
≤
0 for all T
,
V
,
Z
>
0. Thus the solutions of system (2)–(6) limit to M, the largest invariant subset of {
d
U
0
/
d
t
=
0
}. Clearly, it follows from (22) that d
U
0
/
d
t
=
0 if and only if T
=
T
0, V
=
0, and Z
=
0. Noting that M is invariant, for each element of M we have V
=
0 and Z
=
0, and then V
˙
=
0. From (5) we derive that
(23)
0
=
V
˙
=
N
T
δ
T
*
+
N
C
a
C
*
.
Since T
*, C
*
≥
0, then T
*
=
C
*
=
0. Hence d
U
0
/
d
t
=
0 if and only if T
=
T
0, T
*
=
0, C
*
=
0,  V
=
0, and Z
=
0. It follows from LaSalle's invariance principle that the infection-free equilibrium E
0 is GAS when R
0
≤
1.Theorem 3.
For system (2)–(6), if R
1
≤
1
<
R
0, then E
1 is GAS.Proof.
Define the following Lyapunov function:(24)
U
1
=
T
1
(
T
T
1
-
1
-
ln
⁡
(
T
T
1
)
)
+
η
1
T
1
*
(
T
*
T
1
*
-
1
-
ln
⁡
(
T
*
T
1
*
)
)
+
η
2
C
1
*
(
C
*
C
1
*
-
1
-
ln
⁡
(
C
*
C
1
*
)
)
+
η
3
V
1
(
V
V
1
-
1
-
ln
⁡
(
V
V
1
)
)
+
η
4
Z
.
The time derivative of U
1 along the trajectories of (2)–(6) is given by
(25)
d
U
1
d
t
=
(
1
-
T
1
T
)
(
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
)
+
η
1
(
1
-
T
1
*
T
*
)
(
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
-
δ
T
*
)
+
η
2
(
1
-
C
1
*
C
*
)
(
α
(
1
-
ɛ
)
k
T
V
-
a
C
*
)
+
η
3
(
1
-
V
1
V
)
(
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
)
+
η
4
(
g
V
Z
-
μ
Z
)
.
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
)
+
(
1
-
ɛ
)
k
T
1
V
-
η
1
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
T
1
*
T
*
+
η
1
δ
T
1
*
-
η
2
α
(
1
-
ɛ
)
k
T
V
C
1
*
C
*
+
η
2
a
C
1
*
-
δ
η
1
V
1
T
*
V
-
a
η
2
V
1
C
*
V
-
c
η
3
V
+
c
η
3
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
.
Using the following equilibrium conditions for E
1,
(27)
(
1
-
α
)
(
1
-
ɛ
)
k
T
1
V
1
=
δ
T
1
*
,
α
(
1
-
ɛ
)
k
T
1
V
1
=
a
C
1
*
,
c
V
1
=
N
T
δ
T
1
*
+
N
C
a
C
1
*
,
then we have (
1
-
ɛ
)
k
T
1
V
1
=
η
1
δ
T
1
*
+
η
2
a
C
1
* and
(28)
d
U
1
d
t
=
-
d
(
T
-
T
1
)
2
T
+
η
1
δ
T
1
*
(
1
-
T
1
T
)
+
η
2
a
C
1
*
(
1
-
T
1
T
)
-
η
1
δ
T
1
*
T
V
T
1
*
T
1
V
1
T
*
+
η
1
δ
T
1
*
-
η
2
a
C
1
*
T
V
C
1
*
T
1
V
1
C
*
+
η
2
a
C
1
*
-
η
1
δ
T
1
*
V
1
T
*
V
T
1
*
-
η
2
a
C
1
*
V
1
C
*
V
C
1
*
+
η
1
δ
T
1
*
+
η
2
a
C
1
*
+
r
η
3
(
V
1
-
μ
g
)
Z
=
-
d
(
T
-
T
1
)
2
T
+
η
1
δ
T
1
*
[
3
-
T
1
T
-
T
1
*
T
V
T
*
T
1
V
1
-
V
1
T
*
V
T
1
*
]
+
η
2
a
C
1
*
[
3
-
T
1
T
-
C
1
*
T
V
C
*
T
1
V
1
-
C
*
V
1
C
1
*
V
]
+
r
η
3
(
d
g
+
(
1
-
ɛ
)
k
μ
g
(
1
-
ɛ
)
k
)
(
R
1
-
1
)
Z
.
We have that if R
0
>
1, then T
1
,
T
1
*
,
C
1
*
,
V
1
>
0. Since the arithmetical mean is greater than or equal to the geometrical mean, then if R
1
≤
1 then d
U
1
/
d
t
≤
0 for all T
,
T
*
,
C
*
,
V
,
Z
>
0. It can be seen that d
U
1
/
d
t
=
0 if and only if T
=
T
1, T
*
=
T
1
*, C
*
=
C
1
*, V
=
V
1, and Z
=
0. LaSalle's invariance principle implies global stability of E
1.Theorem 4.
For system (2)–(6), if R
0
≤
1, then E
0 is GAS.Proof.
We consider a Lyapunov function(29)
U
2
=
T
2
(
T
T
2
-
1
-
ln
⁡
(
T
T
2
)
)
+
η
1
T
2
*
(
T
*
T
2
*
-
1
-
ln
⁡
(
T
*
T
2
*
)
)
+
η
2
C
2
*
(
C
*
C
2
*
-
1
-
ln
⁡
(
C
*
C
2
*
)
)
+
η
3
V
2
(
V
V
2
-
1
-
ln
⁡
(
V
V
2
)
)
+
η
4
Z
2
(
Z
Z
2
-
1
-
ln
⁡
(
Z
Z
2
)
)
.
Further, function U
2 along the trajectories of system (2)–(6) satisfies
(30)
d
U
2
d
t
=
(
1
-
T
2
T
)
(
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
)
+
η
1
(
1
-
T
2
*
T
*
)
(
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
-
δ
T
*
)
+
η
2
(
1
-
C
2
*
C
*
)
(
α
(
1
-
ɛ
)
k
T
V
-
a
C
*
)
+
η
3
(
1
-
V
2
V
)
(
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
)
+
η
4
(
1
-
Z
2
Z
)
(
g
V
Z
-
μ
Z
)
.
Using the following equilibrium conditions for E
2,
(31)
λ
=
d
T
2
+
(
1
-
ɛ
)
k
T
2
V
2
,
(
1
-
α
)
(
1
-
ɛ
)
k
T
2
V
2
=
δ
T
2
*
,
α
(
1
-
ɛ
)
k
T
2
V
2
=
a
C
2
*
,
c
V
2
+
r
V
2
Z
2
=
N
T
δ
T
2
*
+
N
C
a
C
2
*
,
we get
(32)
d
U
2
d
t
=
-
d
(
T
-
T
2
)
2
T
+
(
1
-
ɛ
)
k
T
2
V
2
(
1
-
T
2
T
)
+
(
1
-
ɛ
)
k
T
2
V
-
η
1
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
T
2
*
T
*
+
δ
η
1
T
2
*
-
η
2
α
(
1
-
ɛ
)
k
T
V
C
2
*
C
*
+
a
η
2
C
2
*
-
δ
η
1
V
2
T
*
V
-
a
η
2
V
2
C
*
V
-
c
η
3
V
+
c
η
3
V
2
+
r
η
4
V
2
Z
-
r
η
4
Z
2
V
+
μ
η
4
Z
2
-
μ
η
4
Z
=
-
d
(
T
-
T
2
)
2
T
+
η
1
δ
T
2
*
(
1
-
T
2
T
)
+
η
2
a
C
2
*
(
1
-
T
2
T
)
-
η
1
δ
T
2
*
T
V
T
2
*
T
2
V
2
T
*
+
η
1
δ
T
2
*
-
η
2
a
C
2
*
T
V
C
2
*
T
2
V
2
C
*
+
η
2
a
C
2
*
-
η
1
δ
T
2
*
V
2
T
*
V
T
2
*
-
η
2
a
C
2
*
V
2
C
*
V
C
2
*
+
η
1
δ
T
2
*
+
η
2
a
C
2
*
=
-
d
(
T
-
T
2
)
2
T
+
η
1
δ
T
2
*
[
3
-
T
2
T
-
T
2
*
T
V
T
*
T
2
V
2
-
V
2
T
*
V
T
2
*
]
+
η
2
a
C
2
*
[
3
-
T
2
T
-
C
2
*
T
V
C
*
T
2
V
2
-
C
*
V
2
C
2
*
V
]
.
Thus, if R
1
>
1, then T
2
,
T
2
*
,
C
2
*
,
V
2 and Z
2
>
0. Since the arithmetical mean is greater than or equal to the geometrical mean, then d
U
2
/
d
t
≤
0. It can be seen that d
U
2
/
d
t
=
0 if and only if T
=
T
2, T
*
=
T
2
*, C
*
=
C
2
*, and V
=
V
2. From (5), if V
=
V
2, then V
˙
=
0 and 0
=
N
T
δ
T
2
*
+
N
C
a
C
2
*
-
c
V
-
r
V
2
Z
=
0, so Z
=
Z
2 and hence d
U
2
/
d
t is equal to zero at E
2. So, the global stability of the equilibrium E
2 follows from LaSalle's invariance principle.


## 2.1. Positive Invariance

We note that model (2)–(6) are biologically acceptable in the sense that no population goes negative. It is straightforward to check the positive invariance of the nonnegative orthant ℝ
+
5 by model (2)–(6) (see, e.g., [6]). In the following, we show the boundedness of the solution of model (2)–(6).Proposition 1.
There exist positive numbersL
i, i
=
1,2
,
3, such that the compact set
(7)
Ω
=
{
(
T
,
T
*
,
C
*
,
V
,
Z
)
∈
ℝ
+
4
:
0
≤
T
,
T
*
,
C
*
≤
L
1
,
r
0
≤
V
≤
L
2
,
0
≤
Z
≤
L
3
{
(
T
,
T
*
,
C
*
,
V
,
Z
)
∈
ℝ
+
4
:
0
≤
T
,
T
*
,
C
*
≤
L
1
,
}
is positively invariant.Proof.
To show the boundedness of the solutions we letG
1
(
t
)
=
T
(
t
)
+
T
*
(
t
)
+
C
*
(
t
); then
(8)
G
˙
1
(
t
)
=
λ
-
d
T
(
t
)
-
(
1
-
ɛ
)
k
T
(
t
)
V
(
t
)
+
(
1
-
α
)
(
1
-
ɛ
)
k
T
(
t
)
V
(
t
)
-
δ
T
*
+
α
(
1
-
ɛ
)
k
T
(
t
)
V
(
t
)
-
a
C
*
(
t
)
≤
λ
-
s
1
G
1
(
t
)
,
where s
1
=
min
⁡
⁡
{
d
,
a
,
δ
}. Hence G
1
(
t
)
≤
L
1, if G
1
(
0
)
≤
L
1 where L
1
=
λ
/
s
1. Since T
(
t
)
>
0, T
*
(
t
)
≥
0, and C
*
(
t
)
≥
0, then 0
≤
T
(
t
), T
*
(
t
), C
*
(
t
)
≤
L
1 if 0
≤
T
(
0
)
+
T
*
(
0
)
+
C
*
(
0
)
≤
L
1. Let G
2
(
t
)
=
V
(
t
)
+
(
r
/
g
)
Z
(
t
); then
(9)
G
˙
2
(
t
)
=
N
T
δ
T
*
(
t
)
+
N
C
a
C
*
(
t
)
-
c
V
(
t
)
-
r
μ
g
Z
(
t
)
≤
(
N
T
δ
+
N
C
a
)
L
1
-
s
2
(
V
(
t
)
+
r
g
Z
(
t
)
)
=
(
N
T
δ
+
N
C
a
)
L
1
-
s
2
G
2
(
t
)
,
where s
2
=
min
⁡
⁡
{
c
,
μ
}. Hence G
2
(
t
)
≤
L
2, if G
2
(
0
)
≤
L
2, where L
2
=
(
N
T
δ
+
N
C
a
)
L
1
/
s
2. Since V
(
t
)
≥
0 and Z
(
t
)
≥
0 then 0
≤
V
(
t
)
≤
L
2 and 0
≤
Z
(
t
)
≤
L
3 if 0
≤
V
(
0
)
+
(
r
/
g
)
Z
(
0
)
≤
L
2, where L
3
=
g
L
2
/
r.


## 2.2. Equilibria

System (2)–(6) always admits an infection-free equilibrium E
0
=
(
T
0
,
0,0
,
0,0
), where T
0
=
λ
/
d. In addition to E
0, the system can have an infected equilibrium without antibody immune response E
1
(
T
1
,
T
1
*
,
C
1
*
,
V
1
,
0
) and an infected equilibrium with antibody immune response E
2
(
T
2
,
T
2
*
,
C
2
*
,
V
2
,
Z
2
) where
(10)
T
1
=
c
(
1
-
ɛ
)
k
[
(
1
-
α
)
N
T
+
α
N
C
]
,
T
1
*
=
(
1
-
α
)
λ
{
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
-
c
}
δ
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
,
C
1
*
=
α
λ
{
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
-
c
}
a
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
,
V
1
=
d
{
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
-
c
}
(
1
-
ɛ
)
k
c
,
T
2
=
λ
g
g
d
+
(
1
-
ɛ
)
k
μ
,
T
2
*
=
(
1
-
α
)
(
1
-
ɛ
)
k
λ
μ
δ
(
d
g
+
(
1
-
ɛ
)
k
μ
)
,
C
2
*
=
α
(
1
-
ɛ
)
k
λ
μ
a
(
d
g
+
(
1
-
ɛ
)
k
μ
)
,
V
2
=
μ
g
,
Z
2
=
c
r
(
d
g
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
c
(
d
g
+
(
1
-
ɛ
)
k
μ
)
-
1
)
.
We discuss the local stability of the infection-free equilibrium E
0. At the infection-free equilibrium E
0
(
T
0
,
0,0
,
0,0
), the system has the Jacobian matrix given by
(11)
J
E
0
=
[
-
d
0
0
-
(
1
-
ɛ
)
k
T
0
0
0
-
δ
0
(
1
-
α
)
(
1
-
ɛ
)
k
T
0
0
0
0
-
a
α
(
1
-
ɛ
)
k
T
0
0
0
δ
N
T
a
N
C
-
c
0
0
0
0
0
-
μ
]
.
The characteristic equation of the Jacobian matrix evaluated at E
0 is
(12)
(
s
+
d
)
(
s
+
μ
)
(
s
3
+
a
1
s
2
+
a
2
s
+
a
3
)
=
0
,
where
(13)
a
1
=
a
+
c
+
δ
,
a
2
=
a
c
+
a
δ
+
c
δ
-
(
1
-
α
)
(
1
-
ɛ
)
k
T
0
N
T
δ
-
α
(
1
-
ɛ
)
k
T
0
N
C
a
,
a
3
=
a
c
δ
(
1
-
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
c
)
.
We observe that (12) has two negative eigenvalues s
1
=
-
d and s
2
=
-
μ. By the Routh-Hurwitz criterion, the remaining three eigenvalues of (12) have negative real parts if a
1
>
0, a
3
>
0, and a
1
a
2
-
a
3
>
0. We have a
1
>
0 and if (
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
/
c
<
1, then a
3
>
0 and
(14)
a
1
a
2
-
a
3
=
a
δ
2
+
a
2
δ
+
2
a
c
δ
+
a
(
a
+
c
)
[
c
-
α
(
1
-
ɛ
)
k
T
0
N
C
]
+
δ
(
δ
+
c
)
[
c
-
(
1
-
α
)
(
1
-
ɛ
)
k
T
0
N
T
]
>
0
.
Now we define the basic reproduction number for system (2)–(6) as
(15)
R
0
=
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
c
.
It follows that the equilibria E
1 and E
2 can be written as
(16)
T
1
=
T
0
R
0
,
T
1
*
=
(
1
-
α
)
λ
δ
(
R
0
-
1
)
R
0
,
C
1
*
=
α
λ
a
(
R
0
-
1
)
R
0
,
V
1
=
d
(
1
-
ɛ
)
k
(
R
0
-
1
)
,
T
2
=
λ
g
g
d
+
(
1
-
ɛ
)
k
μ
,
T
2
*
=
(
1
-
α
)
(
1
-
ɛ
)
k
λ
μ
δ
(
d
g
+
(
1
-
ɛ
)
k
μ
)
,
C
2
*
=
α
(
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 that Z
2
>
0 when d
g
R
0
/
(
d
g
+
(
1
-
ɛ
)
k
μ
)
>
1. Now we define another threshold parameter R
1 as
(17)
R
1
=
R
0
1
+
(
(
1
-
ɛ
)
k
μ
/
d
g
)
.
Clearly R
1
<
R
0.From (16) we have the following statements:(i)
ifR
0
≤
1, then there exists only positive equilibrium E
0;
(ii)
ifR
1
≤
1
<
R
0, then there exist two positive equilibria E
0 and E
1;
(iii)
ifR
1
>
1, then there exist three positive equilibria E
0, E
1, and E
2.


## 2.3. Global Stability Analysis

In this section, we study the global stability of all the equilibria of system (2)–(6) employing the method of Lyapunov function.Theorem 2.
For system (2)–(6), if R
0
≤
1, then E
0 is GAS.Proof.
Define a Lyapunov functionU
0 as follows:
(18)
U
0
=
T
0
(
T
T
0
-
1
-
ln
⁡
(
T
T
0
)
)
+
η
1
T
*
+
η
2
C
*
+
η
3
V
+
η
4
Z
,
where η
i,  i
=
1
,
…
,
4, are positive constants to be determined below. Calculating the derivative of U
0 along the solutions of the system (2)–(6) and applying λ
=
T
0
d, we obtain
(19)
d
U
0
d
t
=
(
1
-
T
0
T
)
(
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
)
+
η
1
(
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
-
δ
T
*
)
+
η
2
(
α
(
1
-
ɛ
)
k
T
V
-
a
C
*
)
+
η
3
(
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
)
+
η
4
(
g
V
Z
-
μ
Z
)
.
Let η
i,  i
=
1
,
…
,
4, be chosen such as
(20)
(
1
-
α
)
η
1
+
α
η
2
=
1
,
η
1
-
N
T
η
3
=
0
,
η
2
-
N
C
η
3
=
0
,
g
η
4
-
r
η
3
=
0
.
The solution of (20) is given by
(21)
η
1
=
N
T
(
1
-
α
)
N
T
+
α
N
C
,
η
2
=
N
C
(
1
-
α
)
N
T
+
α
N
C
,
η
3
=
1
(
1
-
α
)
N
T
+
α
N
C
,
η
4
=
r
g
[
(
1
-
α
)
N
T
+
α
N
C
]
.
The values of η
i,  i
=
1
,
…
,
4, given by (21) will be used throughout the paper. Then
(22)
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
+
η
3
c
(
R
0
-
1
)
V
-
η
4
μ
Z
.
If R
0
≤
1 then d
U
0
/
d
t
≤
0 for all T
,
V
,
Z
>
0. Thus the solutions of system (2)–(6) limit to M, the largest invariant subset of {
d
U
0
/
d
t
=
0
}. Clearly, it follows from (22) that d
U
0
/
d
t
=
0 if and only if T
=
T
0, V
=
0, and Z
=
0. Noting that M is invariant, for each element of M we have V
=
0 and Z
=
0, and then V
˙
=
0. From (5) we derive that
(23)
0
=
V
˙
=
N
T
δ
T
*
+
N
C
a
C
*
.
Since T
*, C
*
≥
0, then T
*
=
C
*
=
0. Hence d
U
0
/
d
t
=
0 if and only if T
=
T
0, T
*
=
0, C
*
=
0,  V
=
0, and Z
=
0. It follows from LaSalle's invariance principle that the infection-free equilibrium E
0 is GAS when R
0
≤
1.Theorem 3.
For system (2)–(6), if R
1
≤
1
<
R
0, then E
1 is GAS.Proof.
Define the following Lyapunov function:(24)
U
1
=
T
1
(
T
T
1
-
1
-
ln
⁡
(
T
T
1
)
)
+
η
1
T
1
*
(
T
*
T
1
*
-
1
-
ln
⁡
(
T
*
T
1
*
)
)
+
η
2
C
1
*
(
C
*
C
1
*
-
1
-
ln
⁡
(
C
*
C
1
*
)
)
+
η
3
V
1
(
V
V
1
-
1
-
ln
⁡
(
V
V
1
)
)
+
η
4
Z
.
The time derivative of U
1 along the trajectories of (2)–(6) is given by
(25)
d
U
1
d
t
=
(
1
-
T
1
T
)
(
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
)
+
η
1
(
1
-
T
1
*
T
*
)
(
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
-
δ
T
*
)
+
η
2
(
1
-
C
1
*
C
*
)
(
α
(
1
-
ɛ
)
k
T
V
-
a
C
*
)
+
η
3
(
1
-
V
1
V
)
(
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
)
+
η
4
(
g
V
Z
-
μ
Z
)
.
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
)
+
(
1
-
ɛ
)
k
T
1
V
-
η
1
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
T
1
*
T
*
+
η
1
δ
T
1
*
-
η
2
α
(
1
-
ɛ
)
k
T
V
C
1
*
C
*
+
η
2
a
C
1
*
-
δ
η
1
V
1
T
*
V
-
a
η
2
V
1
C
*
V
-
c
η
3
V
+
c
η
3
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
.
Using the following equilibrium conditions for E
1,
(27)
(
1
-
α
)
(
1
-
ɛ
)
k
T
1
V
1
=
δ
T
1
*
,
α
(
1
-
ɛ
)
k
T
1
V
1
=
a
C
1
*
,
c
V
1
=
N
T
δ
T
1
*
+
N
C
a
C
1
*
,
then we have (
1
-
ɛ
)
k
T
1
V
1
=
η
1
δ
T
1
*
+
η
2
a
C
1
* and
(28)
d
U
1
d
t
=
-
d
(
T
-
T
1
)
2
T
+
η
1
δ
T
1
*
(
1
-
T
1
T
)
+
η
2
a
C
1
*
(
1
-
T
1
T
)
-
η
1
δ
T
1
*
T
V
T
1
*
T
1
V
1
T
*
+
η
1
δ
T
1
*
-
η
2
a
C
1
*
T
V
C
1
*
T
1
V
1
C
*
+
η
2
a
C
1
*
-
η
1
δ
T
1
*
V
1
T
*
V
T
1
*
-
η
2
a
C
1
*
V
1
C
*
V
C
1
*
+
η
1
δ
T
1
*
+
η
2
a
C
1
*
+
r
η
3
(
V
1
-
μ
g
)
Z
=
-
d
(
T
-
T
1
)
2
T
+
η
1
δ
T
1
*
[
3
-
T
1
T
-
T
1
*
T
V
T
*
T
1
V
1
-
V
1
T
*
V
T
1
*
]
+
η
2
a
C
1
*
[
3
-
T
1
T
-
C
1
*
T
V
C
*
T
1
V
1
-
C
*
V
1
C
1
*
V
]
+
r
η
3
(
d
g
+
(
1
-
ɛ
)
k
μ
g
(
1
-
ɛ
)
k
)
(
R
1
-
1
)
Z
.
We have that if R
0
>
1, then T
1
,
T
1
*
,
C
1
*
,
V
1
>
0. Since the arithmetical mean is greater than or equal to the geometrical mean, then if R
1
≤
1 then d
U
1
/
d
t
≤
0 for all T
,
T
*
,
C
*
,
V
,
Z
>
0. It can be seen that d
U
1
/
d
t
=
0 if and only if T
=
T
1, T
*
=
T
1
*, C
*
=
C
1
*, V
=
V
1, and Z
=
0. LaSalle's invariance principle implies global stability of E
1.Theorem 4.
For system (2)–(6), if R
0
≤
1, then E
0 is GAS.Proof.
We consider a Lyapunov function(29)
U
2
=
T
2
(
T
T
2
-
1
-
ln
⁡
(
T
T
2
)
)
+
η
1
T
2
*
(
T
*
T
2
*
-
1
-
ln
⁡
(
T
*
T
2
*
)
)
+
η
2
C
2
*
(
C
*
C
2
*
-
1
-
ln
⁡
(
C
*
C
2
*
)
)
+
η
3
V
2
(
V
V
2
-
1
-
ln
⁡
(
V
V
2
)
)
+
η
4
Z
2
(
Z
Z
2
-
1
-
ln
⁡
(
Z
Z
2
)
)
.
Further, function U
2 along the trajectories of system (2)–(6) satisfies
(30)
d
U
2
d
t
=
(
1
-
T
2
T
)
(
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
)
+
η
1
(
1
-
T
2
*
T
*
)
(
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
-
δ
T
*
)
+
η
2
(
1
-
C
2
*
C
*
)
(
α
(
1
-
ɛ
)
k
T
V
-
a
C
*
)
+
η
3
(
1
-
V
2
V
)
(
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
)
+
η
4
(
1
-
Z
2
Z
)
(
g
V
Z
-
μ
Z
)
.
Using the following equilibrium conditions for E
2,
(31)
λ
=
d
T
2
+
(
1
-
ɛ
)
k
T
2
V
2
,
(
1
-
α
)
(
1
-
ɛ
)
k
T
2
V
2
=
δ
T
2
*
,
α
(
1
-
ɛ
)
k
T
2
V
2
=
a
C
2
*
,
c
V
2
+
r
V
2
Z
2
=
N
T
δ
T
2
*
+
N
C
a
C
2
*
,
we get
(32)
d
U
2
d
t
=
-
d
(
T
-
T
2
)
2
T
+
(
1
-
ɛ
)
k
T
2
V
2
(
1
-
T
2
T
)
+
(
1
-
ɛ
)
k
T
2
V
-
η
1
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
T
2
*
T
*
+
δ
η
1
T
2
*
-
η
2
α
(
1
-
ɛ
)
k
T
V
C
2
*
C
*
+
a
η
2
C
2
*
-
δ
η
1
V
2
T
*
V
-
a
η
2
V
2
C
*
V
-
c
η
3
V
+
c
η
3
V
2
+
r
η
4
V
2
Z
-
r
η
4
Z
2
V
+
μ
η
4
Z
2
-
μ
η
4
Z
=
-
d
(
T
-
T
2
)
2
T
+
η
1
δ
T
2
*
(
1
-
T
2
T
)
+
η
2
a
C
2
*
(
1
-
T
2
T
)
-
η
1
δ
T
2
*
T
V
T
2
*
T
2
V
2
T
*
+
η
1
δ
T
2
*
-
η
2
a
C
2
*
T
V
C
2
*
T
2
V
2
C
*
+
η
2
a
C
2
*
-
η
1
δ
T
2
*
V
2
T
*
V
T
2
*
-
η
2
a
C
2
*
V
2
C
*
V
C
2
*
+
η
1
δ
T
2
*
+
η
2
a
C
2
*
=
-
d
(
T
-
T
2
)
2
T
+
η
1
δ
T
2
*
[
3
-
T
2
T
-
T
2
*
T
V
T
*
T
2
V
2
-
V
2
T
*
V
T
2
*
]
+
η
2
a
C
2
*
[
3
-
T
2
T
-
C
2
*
T
V
C
*
T
2
V
2
-
C
*
V
2
C
2
*
V
]
.
Thus, if R
1
>
1, then T
2
,
T
2
*
,
C
2
*
,
V
2 and Z
2
>
0. Since the arithmetical mean is greater than or equal to the geometrical mean, then d
U
2
/
d
t
≤
0. It can be seen that d
U
2
/
d
t
=
0 if and only if T
=
T
2, T
*
=
T
2
*, C
*
=
C
2
*, and V
=
V
2. From (5), if V
=
V
2, then V
˙
=
0 and 0
=
N
T
δ
T
2
*
+
N
C
a
C
2
*
-
c
V
-
r
V
2
Z
=
0, so Z
=
Z
2 and hence d
U
2
/
d
t is equal to zero at E
2. So, the global stability of the equilibrium E
2 follows from LaSalle's invariance principle.


## 3. Model with Saturation Incidence Rate

In model (2)–(6), the infection process is characterized by bilinear incidence rate (
1
-
ɛ
)
k
x
v. However, there are a number of reasons why this bilinear incidence can be insufficient to describe infection process in detail (see, e.g., [29–31]). For example, a less than linear response in v could occur when the concentration of viruses becomes higher, where the infectious fraction is high so that exposure is very likely [29]. Experiments reported in [32] strongly suggested that the infection rate of microparasitic infections is an increasing function of the parasite dose and is usually sigmoidal in shape (see, e.g., [33]). In [33], to place the model on more sound biological grounds, Regoes et al. replaced the mass-action infection rate with a dose-dependent infection rates. In this section, the incidence rate is given by a saturation functional response:
(33)
T
˙
=
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
1
+
β
V
,
(34)
T
˙
*
=
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
δ
T
*
,
(35)
C
˙
*
=
α
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
a
C
*
,
(36)
V
˙
=
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
,
(37)
Z
˙
=
g
V
Z
-
μ
Z
,
where β
>
0 is a constant, which represents the saturation infection rate constant.All the variables and parameters have the same meanings as given in model (2)–(6).


### 3.1. Equilibria

Similar to the previous section, we can define two threshold parametersR
0 and R
1 for system (33)–(37) as
(38)
R
0
=
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
c
,
R
1
=
R
0
1
+
(
d
β
μ
+
(
1
-
ɛ
)
k
μ
/
d
g
)
.
Clearly R
1
<
R
0. It is clear that system (33)–(37) has an infection-free equilibrium E
0
=
(
T
0
,
0,0
,
0,0
), where T
0
=
λ
/
d. In addition to E
0, the system can have an infected equilibrium without antibody immune response E
1
(
T
1
,
T
1
*
,
C
1
*
,
V
1
,
0
), where
(39)
T
1
=
β
λ
[
(
1
-
α
)
N
T
+
α
N
C
]
+
c
(
(
1
-
ɛ
)
k
+
d
β
)
[
(
1
-
α
)
N
T
+
α
N
C
]
,
T
1
*
=
(
1
-
α
)
c
d
δ
(
(
1
-
ɛ
)
k
+
d
β
)
[
(
1
-
α
)
N
T
+
α
N
C
]
(
R
0
-
1
)
,
C
1
*
=
α
c
d
a
(
(
1
-
ɛ
)
k
+
d
β
)
[
(
1
-
α
)
N
T
+
α
N
C
]
(
R
0
-
1
)
,
V
1
=
d
(
1
-
ɛ
)
k
+
d
β
(
R
0
-
1
)
,
and infected equilibrium with antibody immune response E
2
(
T
2
,
T
2
*
,
C
2
*
,
V
2
,
Z
2
), where
(40)
T
2
=
λ
(
g
+
β
μ
)
g
d
+
(
1
-
ɛ
)
k
μ
+
d
β
μ
,
T
2
*
=
(
1
-
α
)
(
1
-
ɛ
)
k
λ
μ
δ
(
d
g
+
(
1
-
ɛ
)
k
μ
+
d
β
μ
)
,
C
2
*
=
α
(
1
-
ɛ
)
k
λ
μ
a
(
d
g
+
(
1
-
ɛ
)
k
μ
+
d
β
μ
)
,
V
2
=
μ
g
,
Z
2
=
c
r
(
R
1
-
1
)
.
It is clear from (39) and (40) that(i)
ifR
0
≤
1, then there exists only positive equilibrium E
0;
(ii)
ifR
1
≤
1
<
R
0, then there exist two positive equilibria E
0 and E
1;
(iii)
ifR
1
>
1, then there exist three positive equilibria E
0, E
1, and E
2.


### 3.2. Global Stability Analysis

In this section, we study the global stability of all the equilibria of system (33)–(37) employing the method of Lyapunov function and LaSalle's invariance principle.Theorem 5.
For system (33)–(37), if R
0
≤
1, then E
0 is GAS.Proof.
Define a Lyapunov functionU
0 as follows:
(41)
U
0
=
T
0
(
T
T
0
-
1
-
ln
⁡
(
T
T
0
)
)
+
η
1
T
*
+
η
2
C
*
+
η
3
V
+
η
4
Z
.
Calculating the derivative of U
0 along the solutions of system (33)–(37) and applying λ
=
T
0
d, we obtain
(42)
d
U
0
d
t
=
(
1
-
T
0
T
)
(
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
1
+
β
V
)
+
η
1
(
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
δ
T
*
)
+
η
2
(
α
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
a
C
*
)
+
η
3
(
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
)
+
η
4
(
g
V
Z
-
μ
Z
)
=
(
1
-
T
0
T
)
(
λ
-
d
T
)
+
(
1
-
ɛ
)
k
T
0
V
1
+
β
V
-
c
η
3
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
.
Similar to the proof of Theorem 2, one can easily show that E
0 is GAS when R
0
≤
1.Theorem 6.
For system (33)–(37), if R
1
≤
1
<
R
0, then E
1 is GAS.Proof.
Construct a Lyapunov function as follows:(43)
U
1
=
T
1
(
T
T
1
-
1
-
ln
⁡
(
T
T
1
)
)
+
η
1
T
1
*
(
T
*
T
1
*
-
1
-
ln
⁡
(
T
*
T
1
*
)
)
+
η
2
C
1
*
(
C
*
C
1
*
-
1
-
ln
⁡
(
C
*
C
1
*
)
)
+
η
3
V
1
(
V
V
1
-
1
-
ln
⁡
(
V
V
1
)
)
+
η
4
Z
.
The derivative of U
1 along the trajectories of system (33)–(37) is given by
(44)
d
U
1
d
t
=
(
1
-
T
1
T
)
(
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
1
+
β
V
)
+
η
1
(
1
-
T
1
*
T
*
)
(
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
δ
T
*
)
+
η
2
(
1
-
C
1
*
C
*
)
(
α
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
a
C
*
)
+
η
3
(
1
-
V
1
V
)
(
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
)
+
η
4
(
g
V
Z
-
μ
Z
)
.
Applying λ
=
d
T
1
+
(
(
1
-
ɛ
)
k
T
1
V
1
/
(
1
+
β
V
1
)
) we get
(45)
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
-
η
1
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
1
+
β
V
T
1
*
T
*
+
η
1
δ
T
1
*
-
η
2
α
(
1
-
ɛ
)
k
T
V
1
+
β
V
C
1
*
C
*
+
η
2
a
C
1
*
-
η
1
δ
V
1
T
*
V
-
η
2
a
V
1
C
*
V
-
c
η
3
V
+
c
η
3
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
.
Using the following equilibrium conditions for E
1,
(46)
(
1
-
α
)
(
1
-
ɛ
)
k
T
1
V
1
1
+
β
V
1
=
δ
T
1
*
,
α
(
1
-
ɛ
)
k
T
1
V
1
1
+
β
V
1
=
a
C
1
*
,
c
V
1
=
N
T
δ
T
1
*
+
N
C
a
C
1
*
,
we get
(47)
d
U
1
d
t
=
-
d
(
T
-
T
1
)
2
T
+
η
1
δ
T
1
*
(
1
-
T
1
T
)
+
η
2
a
C
1
*
(
1
-
T
1
T
)
+
(
1
-
ɛ
)
k
T
1
V
1
1
+
β
V
1
[
V
(
1
+
β
V
1
)
V
1
(
1
+
β
V
)
-
V
V
1
]
-
η
1
δ
T
1
*
T
V
T
1
*
(
1
+
β
V
1
)
T
1
V
1
T
*
(
1
+
β
V
)
+
η
1
δ
T
1
*
-
η
2
a
C
1
*
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
+
η
2
a
C
1
*
-
η
1
δ
T
1
*
V
1
T
*
V
T
1
*
-
η
2
a
C
1
*
V
1
C
*
V
C
1
*
+
η
1
δ
T
1
*
+
η
2
a
C
1
*
+
r
η
3
(
V
1
-
μ
g
)
Z
=
-
d
(
T
-
T
1
)
2
T
+
(
1
-
ɛ
)
k
T
1
V
1
1
+
β
V
1
[
-
1

+

V
(
1
+
β
V
1
)
V
1
(
1
+
β
V
)
-
V
V
1

+

1
+
β
V
1
+
β
V
1
]
+
η
1
δ
T
1
*
[
4
-
T
1
T
-
T
V
T
1
*
(
1
+
β
V
1
)
T
1
V
1
T
*
(
1
+
β
V
)
-
V
1
T
*
V
T
1
*
+
η
1
δ
T
1
*
w
-
1
+
β
V
1
+
β
V
1
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
]
+
η
2
a
C
1
*
[
4
-
T
1
T
-
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
-
C
*
V
1
C
1
*
V
+
η
2
a
C
1
*
w
-
1
+
β
V
1
+
β
V
1
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
]
+
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
V
1
(
1
+
β
V
)
(
1
+
β
V
1
)
]
+
η
1
δ
T
1
*
[
4
-
T
1
T
-
T
V
T
1
*
(
1
+
β
V
1
)
T
1
V
1
T
*
(
1
+
β
V
)
-
V
1
T
*
V
T
1
*
+
η
1
δ
T
1
*
w
-
1
+
β
V
1
+
β
V
1
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
]
+
η
2
a
C
1
*
[
4
-
T
1
T
-
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
-
C
*
V
1
C
1
*
V
+
η
2
a
C
1
*
w
-
1
+
β
V
1
+
β
V
1
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
]
+
r
η
3
(
d
g
+
(
1
-
ɛ
)
k
μ
+
d
β
μ
g
(
1
-
ɛ
)
k
+
d
g
β
)
(
R
1
-
1
)
Z
.
We have that if R
1
≤
1
<
R
0, then d
U
1
/
d
t
≤
0 where equality occurs at E
1. LaSalle's invariance principle implies global stability of E
1.Theorem 7.
For system (33)–(37), if R
1
>
1, then E
2 is GAS.Proof.
We consider a Lyapunov function as follows:(48)
U
2
=
T
2
(
T
T
2
-
1
-
ln
⁡
(
T
T
2
)
)
+
η
1
T
2
*
(
T
*
T
2
*
-
1
-
ln
⁡
(
T
*
T
2
*
)
)
+
η
2
C
2
*
(
C
*
C
2
*
-
1
-
ln
⁡
(
C
*
C
2
*
)
)
+
η
3
V
2
(
V
V
2
-
1
-
ln
⁡
(
V
V
2
)
)
+
η
4
Z
2
(
Z
Z
2
-
1
-
ln
⁡
(
Z
Z
2
)
)
.
Further, function U
2 along the trajectories of system (33)–(37) satisfies
(49)
d
U
2
d
t
=
(
1
-
T
2
T
)
(
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
1
+
β
V
)
+
η
1
(
1
-
T
2
*
T
*
)
(
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
δ
T
*
)
+
η
2
(
1
-
C
2
*
C
*
)
(
α
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
a
C
*
)
+
η
3
(
1
-
V
2
V
)
(
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
)
+
η
4
(
1
-
Z
2
Z
)
(
g
V
Z
-
μ
Z
)
.
Using the following equilibrium conditions for E
2,
(50)
λ
=
d
T
2
+
(
1
-
ɛ
)
k
T
2
V
2
1
+
β
V
2
,
δ
T
2
*
=
(
1
-
α
)
(
1
-
ɛ
)
k
T
2
V
2
1
+
β
V
2
,
a
C
2
*
=
α
(
1
-
ɛ
)
k
T
2
V
2
1
+
β
V
2
,
c
V
2
+
r
V
2
Z
2
=
N
T
δ
T
2
*
+
N
C
a
C
2
*
,
we get
(51)
d
U
2
d
t
=
-
d
(
T
-
T
2
)
2
T
+
(
1
-
ɛ
)
k
T
2
V
2
1
+
β
V
2
(
1
-
T
2
T
)
+
(
1
-
ɛ
)
k
T
2
V
1
+
β
V
-
η
1
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
1
+
β
V
T
2
*
T
*
+
η
1
δ
T
2
*
-
η
2
α
(
1
-
ɛ
)
k
T
V
1
+
β
V
C
2
*
C
*
+
η
2
a
C
2
*
-
η
1
δ
V
2
T
*
V
-
η
2
a
V
2
C
*
V
-
η
3
c
V
+
η
3
c
V
2
+
η
3
r
V
2
Z
-
η
4
g
Z
2
V
+
μ
η
4
Z
2
-
μ
η
4
Z
=
-
d
(
T
-
T
2
)
2
T
+
η
1
δ
T
2
*
(
1
-
T
2
T
)
+
η
2
a
C
2
*
(
1
-
T
2
T
)
+
(
1
-
ɛ
)
k
T
2
V
2
1
+
β
V
2
[
V
(
1
+
β
V
2
)
V
2
(
1
+
β
V
)
-
V
V
2
]
-
η
1
δ
T
2
*
T
V
T
2
*
(
1
+
β
V
2
)
T
2
V
2
T
*
(
1
+
β
V
)
+
η
1
δ
T
2
*
-
η
2
a
C
2
*
T
V
C
2
*
(
1
+
β
V
2
)
T
2
V
2
C
*
(
1
+
β
V
)
+
η
2
a
C
2
*
-
η
1
δ
T
2
*
V
2
T
*
V
T
2
*
-
η
2
a
C
2
*
V
2
C
*
V
C
2
*
+
η
1
δ
T
2
*
+
η
2
a
C
2
*
=
-
d
(
T
-
T
2
)
2
T
-
(
1
-
ɛ
)
k
T
2
V
2
1
+
β
V
2
[
β
(
V
-
V
2
)
2
V
2
(
1
+
β
V
)
(
1
+
β
V
2
)
]
+
η
1
δ
T
2
*
[
4
-
T
2
T
-
T
V
T
2
*
(
1
+
β
V
2
)
T
2
V
2
T
*
(
1
+
β
V
)
-
V
2
T
*
V
T
2
*
+
η
1
δ
T
2
*
w
-
1
+
β
V
1
+
β
V
2
T
V
T
2
*
(
1
+
β
V
2
)
T
2
V
2
T
*
(
1
+
β
V
)
]
+
η
2
a
C
2
*
[
4
-
T
2
T
-
T
V
C
2
*
(
1
+
β
V
2
)
T
2
V
2
C
*
(
1
+
β
V
)
-
C
*
V
2
C
2
*
V
+
η
1
δ
T
2
*
w
-
1
+
β
V
1
+
β
V
2
T
V
C
2
*
(
1
+
β
V
2
)
T
2
V
2
C
*
(
1
+
β
V
)
]
.
Similar to the proof of Theorem 4, one can show that E
2 is GAS.


## 3.1. Equilibria

Similar to the previous section, we can define two threshold parametersR
0 and R
1 for system (33)–(37) as
(38)
R
0
=
(
1
-
ɛ
)
k
T
0
[
(
1
-
α
)
N
T
+
α
N
C
]
c
,
R
1
=
R
0
1
+
(
d
β
μ
+
(
1
-
ɛ
)
k
μ
/
d
g
)
.
Clearly R
1
<
R
0. It is clear that system (33)–(37) has an infection-free equilibrium E
0
=
(
T
0
,
0,0
,
0,0
), where T
0
=
λ
/
d. In addition to E
0, the system can have an infected equilibrium without antibody immune response E
1
(
T
1
,
T
1
*
,
C
1
*
,
V
1
,
0
), where
(39)
T
1
=
β
λ
[
(
1
-
α
)
N
T
+
α
N
C
]
+
c
(
(
1
-
ɛ
)
k
+
d
β
)
[
(
1
-
α
)
N
T
+
α
N
C
]
,
T
1
*
=
(
1
-
α
)
c
d
δ
(
(
1
-
ɛ
)
k
+
d
β
)
[
(
1
-
α
)
N
T
+
α
N
C
]
(
R
0
-
1
)
,
C
1
*
=
α
c
d
a
(
(
1
-
ɛ
)
k
+
d
β
)
[
(
1
-
α
)
N
T
+
α
N
C
]
(
R
0
-
1
)
,
V
1
=
d
(
1
-
ɛ
)
k
+
d
β
(
R
0
-
1
)
,
and infected equilibrium with antibody immune response E
2
(
T
2
,
T
2
*
,
C
2
*
,
V
2
,
Z
2
), where
(40)
T
2
=
λ
(
g
+
β
μ
)
g
d
+
(
1
-
ɛ
)
k
μ
+
d
β
μ
,
T
2
*
=
(
1
-
α
)
(
1
-
ɛ
)
k
λ
μ
δ
(
d
g
+
(
1
-
ɛ
)
k
μ
+
d
β
μ
)
,
C
2
*
=
α
(
1
-
ɛ
)
k
λ
μ
a
(
d
g
+
(
1
-
ɛ
)
k
μ
+
d
β
μ
)
,
V
2
=
μ
g
,
Z
2
=
c
r
(
R
1
-
1
)
.
It is clear from (39) and (40) that(i)
ifR
0
≤
1, then there exists only positive equilibrium E
0;
(ii)
ifR
1
≤
1
<
R
0, then there exist two positive equilibria E
0 and E
1;
(iii)
ifR
1
>
1, then there exist three positive equilibria E
0, E
1, and E
2.


## 3.2. Global Stability Analysis

In this section, we study the global stability of all the equilibria of system (33)–(37) employing the method of Lyapunov function and LaSalle's invariance principle.Theorem 5.
For system (33)–(37), if R
0
≤
1, then E
0 is GAS.Proof.
Define a Lyapunov functionU
0 as follows:
(41)
U
0
=
T
0
(
T
T
0
-
1
-
ln
⁡
(
T
T
0
)
)
+
η
1
T
*
+
η
2
C
*
+
η
3
V
+
η
4
Z
.
Calculating the derivative of U
0 along the solutions of system (33)–(37) and applying λ
=
T
0
d, we obtain
(42)
d
U
0
d
t
=
(
1
-
T
0
T
)
(
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
1
+
β
V
)
+
η
1
(
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
δ
T
*
)
+
η
2
(
α
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
a
C
*
)
+
η
3
(
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
)
+
η
4
(
g
V
Z
-
μ
Z
)
=
(
1
-
T
0
T
)
(
λ
-
d
T
)
+
(
1
-
ɛ
)
k
T
0
V
1
+
β
V
-
c
η
3
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
.
Similar to the proof of Theorem 2, one can easily show that E
0 is GAS when R
0
≤
1.Theorem 6.
For system (33)–(37), if R
1
≤
1
<
R
0, then E
1 is GAS.Proof.
Construct a Lyapunov function as follows:(43)
U
1
=
T
1
(
T
T
1
-
1
-
ln
⁡
(
T
T
1
)
)
+
η
1
T
1
*
(
T
*
T
1
*
-
1
-
ln
⁡
(
T
*
T
1
*
)
)
+
η
2
C
1
*
(
C
*
C
1
*
-
1
-
ln
⁡
(
C
*
C
1
*
)
)
+
η
3
V
1
(
V
V
1
-
1
-
ln
⁡
(
V
V
1
)
)
+
η
4
Z
.
The derivative of U
1 along the trajectories of system (33)–(37) is given by
(44)
d
U
1
d
t
=
(
1
-
T
1
T
)
(
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
1
+
β
V
)
+
η
1
(
1
-
T
1
*
T
*
)
(
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
δ
T
*
)
+
η
2
(
1
-
C
1
*
C
*
)
(
α
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
a
C
*
)
+
η
3
(
1
-
V
1
V
)
(
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
)
+
η
4
(
g
V
Z
-
μ
Z
)
.
Applying λ
=
d
T
1
+
(
(
1
-
ɛ
)
k
T
1
V
1
/
(
1
+
β
V
1
)
) we get
(45)
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
-
η
1
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
1
+
β
V
T
1
*
T
*
+
η
1
δ
T
1
*
-
η
2
α
(
1
-
ɛ
)
k
T
V
1
+
β
V
C
1
*
C
*
+
η
2
a
C
1
*
-
η
1
δ
V
1
T
*
V
-
η
2
a
V
1
C
*
V
-
c
η
3
V
+
c
η
3
V
1
+
r
η
3
V
1
Z
-
μ
η
4
Z
.
Using the following equilibrium conditions for E
1,
(46)
(
1
-
α
)
(
1
-
ɛ
)
k
T
1
V
1
1
+
β
V
1
=
δ
T
1
*
,
α
(
1
-
ɛ
)
k
T
1
V
1
1
+
β
V
1
=
a
C
1
*
,
c
V
1
=
N
T
δ
T
1
*
+
N
C
a
C
1
*
,
we get
(47)
d
U
1
d
t
=
-
d
(
T
-
T
1
)
2
T
+
η
1
δ
T
1
*
(
1
-
T
1
T
)
+
η
2
a
C
1
*
(
1
-
T
1
T
)
+
(
1
-
ɛ
)
k
T
1
V
1
1
+
β
V
1
[
V
(
1
+
β
V
1
)
V
1
(
1
+
β
V
)
-
V
V
1
]
-
η
1
δ
T
1
*
T
V
T
1
*
(
1
+
β
V
1
)
T
1
V
1
T
*
(
1
+
β
V
)
+
η
1
δ
T
1
*
-
η
2
a
C
1
*
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
+
η
2
a
C
1
*
-
η
1
δ
T
1
*
V
1
T
*
V
T
1
*
-
η
2
a
C
1
*
V
1
C
*
V
C
1
*
+
η
1
δ
T
1
*
+
η
2
a
C
1
*
+
r
η
3
(
V
1
-
μ
g
)
Z
=
-
d
(
T
-
T
1
)
2
T
+
(
1
-
ɛ
)
k
T
1
V
1
1
+
β
V
1
[
-
1

+

V
(
1
+
β
V
1
)
V
1
(
1
+
β
V
)
-
V
V
1

+

1
+
β
V
1
+
β
V
1
]
+
η
1
δ
T
1
*
[
4
-
T
1
T
-
T
V
T
1
*
(
1
+
β
V
1
)
T
1
V
1
T
*
(
1
+
β
V
)
-
V
1
T
*
V
T
1
*
+
η
1
δ
T
1
*
w
-
1
+
β
V
1
+
β
V
1
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
]
+
η
2
a
C
1
*
[
4
-
T
1
T
-
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
-
C
*
V
1
C
1
*
V
+
η
2
a
C
1
*
w
-
1
+
β
V
1
+
β
V
1
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
]
+
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
V
1
(
1
+
β
V
)
(
1
+
β
V
1
)
]
+
η
1
δ
T
1
*
[
4
-
T
1
T
-
T
V
T
1
*
(
1
+
β
V
1
)
T
1
V
1
T
*
(
1
+
β
V
)
-
V
1
T
*
V
T
1
*
+
η
1
δ
T
1
*
w
-
1
+
β
V
1
+
β
V
1
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
]
+
η
2
a
C
1
*
[
4
-
T
1
T
-
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
-
C
*
V
1
C
1
*
V
+
η
2
a
C
1
*
w
-
1
+
β
V
1
+
β
V
1
T
V
C
1
*
(
1
+
β
V
1
)
T
1
V
1
C
*
(
1
+
β
V
)
]
+
r
η
3
(
d
g
+
(
1
-
ɛ
)
k
μ
+
d
β
μ
g
(
1
-
ɛ
)
k
+
d
g
β
)
(
R
1
-
1
)
Z
.
We have that if R
1
≤
1
<
R
0, then d
U
1
/
d
t
≤
0 where equality occurs at E
1. LaSalle's invariance principle implies global stability of E
1.Theorem 7.
For system (33)–(37), if R
1
>
1, then E
2 is GAS.Proof.
We consider a Lyapunov function as follows:(48)
U
2
=
T
2
(
T
T
2
-
1
-
ln
⁡
(
T
T
2
)
)
+
η
1
T
2
*
(
T
*
T
2
*
-
1
-
ln
⁡
(
T
*
T
2
*
)
)
+
η
2
C
2
*
(
C
*
C
2
*
-
1
-
ln
⁡
(
C
*
C
2
*
)
)
+
η
3
V
2
(
V
V
2
-
1
-
ln
⁡
(
V
V
2
)
)
+
η
4
Z
2
(
Z
Z
2
-
1
-
ln
⁡
(
Z
Z
2
)
)
.
Further, function U
2 along the trajectories of system (33)–(37) satisfies
(49)
d
U
2
d
t
=
(
1
-
T
2
T
)
(
λ
-
d
T
-
(
1
-
ɛ
)
k
T
V
1
+
β
V
)
+
η
1
(
1
-
T
2
*
T
*
)
(
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
δ
T
*
)
+
η
2
(
1
-
C
2
*
C
*
)
(
α
(
1
-
ɛ
)
k
T
V
1
+
β
V
-
a
C
*
)
+
η
3
(
1
-
V
2
V
)
(
N
T
δ
T
*
+
N
C
a
C
*
-
c
V
-
r
V
Z
)
+
η
4
(
1
-
Z
2
Z
)
(
g
V
Z
-
μ
Z
)
.
Using the following equilibrium conditions for E
2,
(50)
λ
=
d
T
2
+
(
1
-
ɛ
)
k
T
2
V
2
1
+
β
V
2
,
δ
T
2
*
=
(
1
-
α
)
(
1
-
ɛ
)
k
T
2
V
2
1
+
β
V
2
,
a
C
2
*
=
α
(
1
-
ɛ
)
k
T
2
V
2
1
+
β
V
2
,
c
V
2
+
r
V
2
Z
2
=
N
T
δ
T
2
*
+
N
C
a
C
2
*
,
we get
(51)
d
U
2
d
t
=
-
d
(
T
-
T
2
)
2
T
+
(
1
-
ɛ
)
k
T
2
V
2
1
+
β
V
2
(
1
-
T
2
T
)
+
(
1
-
ɛ
)
k
T
2
V
1
+
β
V
-
η
1
(
1
-
α
)
(
1
-
ɛ
)
k
T
V
1
+
β
V
T
2
*
T
*
+
η
1
δ
T
2
*
-
η
2
α
(
1
-
ɛ
)
k
T
V
1
+
β
V
C
2
*
C
*
+
η
2
a
C
2
*
-
η
1
δ
V
2
T
*
V
-
η
2
a
V
2
C
*
V
-
η
3
c
V
+
η
3
c
V
2
+
η
3
r
V
2
Z
-
η
4
g
Z
2
V
+
μ
η
4
Z
2
-
μ
η
4
Z
=
-
d
(
T
-
T
2
)
2
T
+
η
1
δ
T
2
*
(
1
-
T
2
T
)
+
η
2
a
C
2
*
(
1
-
T
2
T
)
+
(
1
-
ɛ
)
k
T
2
V
2
1
+
β
V
2
[
V
(
1
+
β
V
2
)
V
2
(
1
+
β
V
)
-
V
V
2
]
-
η
1
δ
T
2
*
T
V
T
2
*
(
1
+
β
V
2
)
T
2
V
2
T
*
(
1
+
β
V
)
+
η
1
δ
T
2
*
-
η
2
a
C
2
*
T
V
C
2
*
(
1
+
β
V
2
)
T
2
V
2
C
*
(
1
+
β
V
)
+
η
2
a
C
2
*
-
η
1
δ
T
2
*
V
2
T
*
V
T
2
*
-
η
2
a
C
2
*
V
2
C
*
V
C
2
*
+
η
1
δ
T
2
*
+
η
2
a
C
2
*
=
-
d
(
T
-
T
2
)
2
T
-
(
1
-
ɛ
)
k
T
2
V
2
1
+
β
V
2
[
β
(
V
-
V
2
)
2
V
2
(
1
+
β
V
)
(
1
+
β
V
2
)
]
+
η
1
δ
T
2
*
[
4
-
T
2
T
-
T
V
T
2
*
(
1
+
β
V
2
)
T
2
V
2
T
*
(
1
+
β
V
)
-
V
2
T
*
V
T
2
*
+
η
1
δ
T
2
*
w
-
1
+
β
V
1
+
β
V
2
T
V
T
2
*
(
1
+
β
V
2
)
T
2
V
2
T
*
(
1
+
β
V
)
]
+
η
2
a
C
2
*
[
4
-
T
2
T
-
T
V
C
2
*
(
1
+
β
V
2
)
T
2
V
2
C
*
(
1
+
β
V
)
-
C
*
V
2
C
2
*
V
+
η
1
δ
T
2
*
w
-
1
+
β
V
1
+
β
V
2
T
V
C
2
*
(
1
+
β
V
2
)
T
2
V
2
C
*
(
1
+
β
V
)
]
.
Similar to the proof of Theorem 4, one can show that E
2 is GAS.


## 4. Numerical Simulations

We now use simple numerical simulations to illustrate our theoretical results for the two models. In both models we will fix the following data:λ
=
10 mm−3 day−1,d
=
0.01 day−1, k
=
0.001 mm3 day−1,  δ
=
0.5 day−1,  α
=
0.5,  a
=
0.1 day-
1, c
=
3 day−1, N
T
=
10, N
C
=
5, r
=
0.01 mm3 day−1, and μ
=
0.1 day−1. The other parameters will be chosen below. All computations were carried out by MATLAB.


### 4.1. Model with Bilinear Incidence Rate

In this section, we perform simulation results for model (2)–(6) to check our theoretical results given in Theorems 2–4. We have the following cases.(i)
R
0
≤
1. We choose ɛ
=
0.63 and g
=
0.01 mm3 day−1. Using these data we compute R
0
=
0.92 and R
1
=
0.672. Figures 1, 2, 3, 4, and 5 show that the numerical results are consistent with Theorem 2. We can see that, the concentration of uninfected cells is increased and converges to its normal value λ
/
d
=
1000 mm−3, while the concentrations of short-lived infected cells, chronically infected cells, free viruses, and antibody immune cells are decaying and tend to zero.
(ii)
R
1
≤
1
<
R
0. We take ɛ
=
0 and g
=
0.005 mm3 day−1. In this case, R
0
=
2.5 and R
1
=
0.833. Figures 1–5 show that the numerical results are consistent with Theorem 3. We can see that the trajectory of the system will tend to the infected equilibrium without antibody immune response E
1
(
400,6
,
27.77,15,0
). In this case, the infection becomes chronic but with no persistent antibody immune response.
(iii)
R
1
>
1. We choose ɛ
=
0 and g
=
0.01 mm3 day−1. Then we compute R
0
=
2.5 and R
1
=
1.25. From Figures 1–5 we can see that our simulation results are consistent with the theoretical results of Theorem 4. We observe that the trajectory of the system will tend to the infected equilibrium with antibody immune response E
2
(
500.04,5
,
23.15,10,57.03
). In this case, the infection becomes chronic but with persistent antibody immune response.Figure 1
The evolution of uninfected cells for model (2)–(6).Figure 2
The evolution of short-lived infected cells for model (2)–(6).Figure 3
The evolution of chronically infected cells for model (2)–(6).Figure 4
The evolution of free viruses for model (2)–(6).Figure 5
The evolution of antibody immune cells for model (2)–(6).We note that the values of the parametersg, r, and μ have no impact on the value of R
0, since R
0 is independent of those parameters. This fact seems to suggest that antibodies do not play a role in eliminating the viruses. From the definition of R
1, we can see that R
1 can be increased by increasing g or decreasing μ.Figures1 and 4 show that the presence of antibody immune response (i.e., R
1
>
1) reduces the concentration of free viruses and increases the concentration of uninfected cells. This can be seen by comparing the virus and uninfected cell components in the equilibria E
1 and E
2 under the condition R
1
>
1. For model (2)–(6), simple calculation shows that
(52)
V
1
-
V
2
=
(
d
g
+
(
1
-
ɛ
)
k
μ
g
(
1
-
ɛ
)
k
)
(
R
1
-
1
)
.
It follows that if R
1
>
1, then V
2
<
V
1. From (2) and at any equilibrium point E
¯
(
T
¯
,
T
¯
*
,
C
¯
*
,
V
¯
,
Z
¯
) we have
(53)
T
¯
=
λ
d
+
(
1
-
ɛ
)
k
V
¯
.
Clearly, T
¯ is a decreasing function of V
¯. This yields that if R
1
>
1, then V
2
<
V
1 and T
2
>
T
1.


### 4.2. Model with Saturation Functional Response

In this section, we perform simulation results to check Theorems5–7. The parameter β is chosen as α
=
0.2 mm3. We have the following cases.(i)
R
0
≤
1. We take ɛ
=
0.63 and g
=
0.01 mm3 day−1. Using these data, we compute R
0
=
0.92 and R
1
=
0.273. The simulation results of this case are shown in Figures 6, 7, 8, 9, and 10. We can see that the numerical results are consistent with Theorem 5. It is observed that the viruses will be cleared and the uninfected cells will return to their normal value.
(ii)
R
1
≤
1
<
R
0. To satisfy this condition, we take ɛ
=
0 and g
=
0.005 mm3 day−1. This will give R
0
=
2.5 and R
1
=
0.833. Figures 6–10 show that the numerical results are consistent with Theorem 6. We see that the infected equilibrium E
1
(
800
,
2
,
9.25
,
5
,
0
) is GAS, and the infection becomes chronic but with no persistent antibody immune response.
(iii)
R
1
>
1. This condition is satisfied by choosing ɛ
=
0 and g
=
0.01 mm3 day−1. This yields R
0
=
2.5 and R
1
=
1.25. Figures 6–10 demonstrate the global stability of E
2
(
832.58
,
1.67
,
7.71
,
3.34
,
74.55
). Then, the infection becomes chronic but with persistent antibody immune response.Figure 6
The evolution of uninfected cells for model (33)–(37).Figure 7
The evolution of short-lived infected cells for model (33)–(37).Figure 8
The evolution of chronically infected cells for model (33)–(37).Figure 9
The evolution of free viruses for model (33)–(37).Figure 10
The evolution of antibody immune cells for model (33)–(37).From the definition of the parameterR
0, we can see that the value of the saturation infection rate constant β has no impact on the value of R
0. This means that saturation does not play a role in eliminating the virus. From the definition of R
1, we can see that R
1 can be increased by increasing g or decreasing μ and β.Figures6 and 9 show that if R
1
>
1 the antibody immune response reduces the concentration of free viruses and increases the concentration of uninfected cells. For model (33)–(37), simple calculation shows that
(54)
V
1
-
V
2
=
(
d
g
+
(
1
-
ɛ
)
k
μ
+
d
β
μ
g
(
1
-
ɛ
)
k
+
d
g
β
)
(
R
1
-
1
)
.
As a result, if R
1
>
1, then V
2
<
V
1. From (33) and at any equilibrium point E
¯
(
T
¯
,
T
¯
*
,
C
¯
*
,
V
¯
,
Z
¯
) we have
(55)
T
¯
=
(
1
+
β
V
¯
)
λ
d
+
(
1
-
ɛ
)
k
V
¯
+
d
V
¯
β
,
d
T
¯
d
V
¯
=
-
(
1
-
ɛ
)
k
λ
(
d
+
(
1
-
ɛ
)
k
V
¯
+
d
V
¯
β
)
2
.
Then, T
¯ is a decreasing function of V
¯. It follows that if R
1
>
1 then V
2
<
V
1 and T
2
>
T
1.


## 4.1. Model with Bilinear Incidence Rate

In this section, we perform simulation results for model (2)–(6) to check our theoretical results given in Theorems 2–4. We have the following cases.(i)
R
0
≤
1. We choose ɛ
=
0.63 and g
=
0.01 mm3 day−1. Using these data we compute R
0
=
0.92 and R
1
=
0.672. Figures 1, 2, 3, 4, and 5 show that the numerical results are consistent with Theorem 2. We can see that, the concentration of uninfected cells is increased and converges to its normal value λ
/
d
=
1000 mm−3, while the concentrations of short-lived infected cells, chronically infected cells, free viruses, and antibody immune cells are decaying and tend to zero.
(ii)
R
1
≤
1
<
R
0. We take ɛ
=
0 and g
=
0.005 mm3 day−1. In this case, R
0
=
2.5 and R
1
=
0.833. Figures 1–5 show that the numerical results are consistent with Theorem 3. We can see that the trajectory of the system will tend to the infected equilibrium without antibody immune response E
1
(
400,6
,
27.77,15,0
). In this case, the infection becomes chronic but with no persistent antibody immune response.
(iii)
R
1
>
1. We choose ɛ
=
0 and g
=
0.01 mm3 day−1. Then we compute R
0
=
2.5 and R
1
=
1.25. From Figures 1–5 we can see that our simulation results are consistent with the theoretical results of Theorem 4. We observe that the trajectory of the system will tend to the infected equilibrium with antibody immune response E
2
(
500.04,5
,
23.15,10,57.03
). In this case, the infection becomes chronic but with persistent antibody immune response.Figure 1
The evolution of uninfected cells for model (2)–(6).Figure 2
The evolution of short-lived infected cells for model (2)–(6).Figure 3
The evolution of chronically infected cells for model (2)–(6).Figure 4
The evolution of free viruses for model (2)–(6).Figure 5
The evolution of antibody immune cells for model (2)–(6).We note that the values of the parametersg, r, and μ have no impact on the value of R
0, since R
0 is independent of those parameters. This fact seems to suggest that antibodies do not play a role in eliminating the viruses. From the definition of R
1, we can see that R
1 can be increased by increasing g or decreasing μ.Figures1 and 4 show that the presence of antibody immune response (i.e., R
1
>
1) reduces the concentration of free viruses and increases the concentration of uninfected cells. This can be seen by comparing the virus and uninfected cell components in the equilibria E
1 and E
2 under the condition R
1
>
1. For model (2)–(6), simple calculation shows that
(52)
V
1
-
V
2
=
(
d
g
+
(
1
-
ɛ
)
k
μ
g
(
1
-
ɛ
)
k
)
(
R
1
-
1
)
.
It follows that if R
1
>
1, then V
2
<
V
1. From (2) and at any equilibrium point E
¯
(
T
¯
,
T
¯
*
,
C
¯
*
,
V
¯
,
Z
¯
) we have
(53)
T
¯
=
λ
d
+
(
1
-
ɛ
)
k
V
¯
.
Clearly, T
¯ is a decreasing function of V
¯. This yields that if R
1
>
1, then V
2
<
V
1 and T
2
>
T
1.


## 4.2. Model with Saturation Functional Response

In this section, we perform simulation results to check Theorems5–7. The parameter β is chosen as α
=
0.2 mm3. We have the following cases.(i)
R
0
≤
1. We take ɛ
=
0.63 and g
=
0.01 mm3 day−1. Using these data, we compute R
0
=
0.92 and R
1
=
0.273. The simulation results of this case are shown in Figures 6, 7, 8, 9, and 10. We can see that the numerical results are consistent with Theorem 5. It is observed that the viruses will be cleared and the uninfected cells will return to their normal value.
(ii)
R
1
≤
1
<
R
0. To satisfy this condition, we take ɛ
=
0 and g
=
0.005 mm3 day−1. This will give R
0
=
2.5 and R
1
=
0.833. Figures 6–10 show that the numerical results are consistent with Theorem 6. We see that the infected equilibrium E
1
(
800
,
2
,
9.25
,
5
,
0
) is GAS, and the infection becomes chronic but with no persistent antibody immune response.
(iii)
R
1
>
1. This condition is satisfied by choosing ɛ
=
0 and g
=
0.01 mm3 day−1. This yields R
0
=
2.5 and R
1
=
1.25. Figures 6–10 demonstrate the global stability of E
2
(
832.58
,
1.67
,
7.71
,
3.34
,
74.55
). Then, the infection becomes chronic but with persistent antibody immune response.Figure 6
The evolution of uninfected cells for model (33)–(37).Figure 7
The evolution of short-lived infected cells for model (33)–(37).Figure 8
The evolution of chronically infected cells for model (33)–(37).Figure 9
The evolution of free viruses for model (33)–(37).Figure 10
The evolution of antibody immune cells for model (33)–(37).From the definition of the parameterR
0, we can see that the value of the saturation infection rate constant β has no impact on the value of R
0. This means that saturation does not play a role in eliminating the virus. From the definition of R
1, we can see that R
1 can be increased by increasing g or decreasing μ and β.Figures6 and 9 show that if R
1
>
1 the antibody immune response reduces the concentration of free viruses and increases the concentration of uninfected cells. For model (33)–(37), simple calculation shows that
(54)
V
1
-
V
2
=
(
d
g
+
(
1
-
ɛ
)
k
μ
+
d
β
μ
g
(
1
-
ɛ
)
k
+
d
g
β
)
(
R
1
-
1
)
.
As a result, if R
1
>
1, then V
2
<
V
1. From (33) and at any equilibrium point E
¯
(
T
¯
,
T
¯
*
,
C
¯
*
,
V
¯
,
Z
¯
) we have
(55)
T
¯
=
(
1
+
β
V
¯
)
λ
d
+
(
1
-
ɛ
)
k
V
¯
+
d
V
¯
β
,
d
T
¯
d
V
¯
=
-
(
1
-
ɛ
)
k
λ
(
d
+
(
1
-
ɛ
)
k
V
¯
+
d
V
¯
β
)
2
.
Then, T
¯ is a decreasing function of V
¯. It follows that if R
1
>
1 then V
2
<
V
1 and T
2
>
T
1.


## 5. Conclusions

In this paper, we have proposed two virus infection models with antibody immune response taking into account the chronically infected cells. In the first model we have assumed that the incidence rate of infection is bilinear while in the second model the incidence rate is given by saturation functional response. We have shown that the dynamics of the models are fully determined by two threshold parametersR
0 and R
1. The parameter R
0 determines whether a chronic infection can be established while R
1 determines whether a persistent antibody response can be established. By constructing Lyapunov function and using LaSalle's invariance principle, we have investigated the global stability of all equilibria of the two models. We have proven that if R
0
≤
1 then the infection-free equilibrium E
0 is GAS, and the viruses are cleared. If R
1
≤
1
<
R
0, then the infected equilibrium without antibody immune response E
1 exists and it is GAS, and the infection becomes chronic but with no persistent antibody immune response. If R
1
>
1, then the infected equilibrium with antibody immune response E
2 exists and it is GAS, and the infection is chronic with persistent antibody immune response. Numerical simulations have been performed for the two models. Our simulation results confirm the analytic results given in Theorems 2–7.



---

*Source: 650371-2014-04-03.xml*