# Modeling Saturated Diagnosis and Vaccination in Reducing HIV/AIDS Infection

Authors: Can Chen; Yanni Xiao
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A mathematical model is proposed to consider the effects of saturated diagnosis and vaccination on HIV/AIDS infection. By employing center manifold theory, we prove that there exists a backward bifurcation which suggests that the disease cannot be eradicated even if the basic reproduction number is less than unity. Global stability of the disease-free equilibrium is investigated for appropriate conditions. When the basic reproduction number is greater than unity, the system is uniformly persistent. The proposed model is applied to describe HIV infection among injecting drug users (IDUs) in Yunnan province, China. Numerical studies indicate that new cases and prevalence are sensitive to transmission rate, vaccination rate, and vaccine efficacy. The findings suggest that increasing vaccination rate and vaccine efficacy and enhancing interventions like reducing share injectors can greatly reduce the transmission of HIV among IDUs in Yunnan province, China.


---

## Body

## 1. Introduction

Acquired immunodeficiency syndrome (AIDS) is spreading rapidly in the world ever since it was firstly detected in 1981 and continues to threaten the health of human seriously, especially among sex workers and injecting drug users. Furthermore, AIDS also influences the economy of many countries which has attracted great attention of governments. For such a severe scenario, the governments have taken intervention measures to reduce HIV transmission.Mathematical models play a vital role in gaining a quantitative insight into HIV transmission dynamics and suggesting the effective control strategies. In order to study the effect of various intervention strategies on HIV transmission, extensive mathematical models have been formulated. Traditionally, models of HIV/AIDS dynamics often incorporate staged progression (see, e.g., [1–3]), but these did not include any control measures. Hyman et al. [4] extended these models to consider screening and contact tracing and discussed which strategy would slow infectiousness. Compartmental models with staged progression that incorporate the imperfect vaccine were constructed in [5] to predict HIV epidemic, but they did not consider diagnosis. Elbasha and Gumel [6] considered that a proportion of new recruits are vaccinated and upon becoming infected with HIV, susceptible and vaccinated individuals enter the classes of infected and vaccine infected people, separately. They showed the existence of backward bifurcation via numerical simulations. Sharomi et al. [7] explored the role of the choice of incidence function in HIV models formulated in [6] and obtained that the phenomenon of backward bifurcation can be removed by substituting the standard incidence function with a mass action incidence. In South Africa, testing and screening campaign was launched for HIV; Nyabadza and Mukandavire [8] analyzed their effects by developing HIV models. More recently, a model of HIV/AIDS with diagnosis was presented in [9]. The authors estimated parameter values and predicted its transmission in China in the next few years.The majority of mathematical models consider only one control strategy, vaccination or diagnosis, for instance [5, 9]; however, curbing HIV/AIDS infection needs comprehensive strategies, since, under the serious threat of HIV, it may be more rational to adopt various measures for different high risk groups. These motivate us to consider two combined intervention measures, vaccination and diagnosis. Furthermore, due to the limited resources, we then choose a nonlinear function which can be used to describe saturation effect. We use a parameter h, representing the half saturation constant, in the diagnosis function to measure the effect of HIV individuals being late for diagnosis [10]. When the number of infected individuals I is low, the number of actual per capita diagnosed individuals is proportional to I, whereas when the number of infected individuals I is sufficiently large, there is a saturation effect which makes the number of diagnosed individuals approach to be constant due to the limitation of human and economic power. The number of new diagnosed cases per unit time is saturated with the total infected population.The paper is organized as follows. The model is formulated in Section2. The existence of backward bifurcation and the stability of the disease-free equilibrium are discussed in Section 3. In Section 4, uniform persistence of the model is investigated. Numerical simulation results are concluded in Section 5. In Section 6, we give a brief summary.


## 2. The Model

The model describes the spread of HIV/AIDS in a high risk population. The total high risk population size at timet, denoted by (N
(
t
)), is subdivided into susceptible individuals (S
(
t
)), vaccinated susceptible individuals (V
(
t
)), HIV infected but not yet diagnosed individuals (I
(
t
)), diagnosed HIV-positive individuals (D
(
t
)), and those who have developed AIDS (A
(
t
)), so that N
(
t
)
=
S
(
t
)
+
V
(
t
)
+
I
(
t
)
+
D
(
t
)
+
A
(
t
).The equations of the model are(1)
S
′
=
Π
+
ω
V
-
λ
(
t
)
S
-
ξ
S
-
μ
S
,
V
′
=
ξ
S
-
(
1
-
ϵ
)
λ
(
t
)
V
-
ω
V
-
μ
V
,
I
′
=
λ
(
t
)
(
S
+
(
1
-
ϵ
)
V
)
-
σ
1
I
-
q
I
1
+
h
I
-
μ
I
,
D
′
=
q
I
1
+
h
I
-
σ
2
D
-
μ
D
,
A
′
=
σ
1
I
+
σ
2
D
-
(
μ
+
ψ
)
A
,
where the incidence rate λ
(
t
)
=
β
(
(
I
(
t
)
+
η
1
D
(
t
)
+
η
2
A
(
t
)
)
/
N
(
t
)
), β denotes the transmission rate, and η
1 and η
2 illustrate the modification factors in the transmission coefficient of diagnosed HIV-positive individuals and AIDS patients, respectively. People enter into the susceptible class at a rate Π, become infected at a rate λ
(
t
)
S, and become vaccinated at a rate ξ. Also μ is the natural death rate; ω denotes the waning rate of vaccine; ϵ represents the vaccine efficacy; q is the diagnosis rate; σ
1 and σ
2 are the progression rate to diagnose HIV-positive individuals and AIDS patients, respectively; ψ is the disease-induced death rate.Since the model monitors change in the human population, the variables and parameters are assumed to be nonnegative for allt
≥
0. The system will be analyzed in a suitable feasible region Ω
⊆
R
+
5, where Ω
=
{
(
S
,
V
,
I
,
D
,
A
)
∈
R
+
5
∣
S
+
V
+
I
+
D
+
A
≤
Π
/
μ
}. We can easily prove that the solutions of system (1) with nonnegative initial conditions remain nonnegative, and the feasible region Ω is positively invariant and attracting with respect to system (1) for all t
>
0.


## 3. Model Analysis

### 3.1. Disease-Free Equilibrium and the Basic Reproduction Number

Model (1) has a disease-free equilibrium (DFE), obtained by setting the right-hand sides of system (1) to zero, represented as
(2)
E
0
:
(
S
¯
,
V
¯
,
I
¯
,
D
¯
,
A
¯
)

=

(

Π
(
μ
+
ω
)
μ
(
μ
+
ξ
+
ω
)
,
Π
ξ
μ
(
μ
+
ξ
+
ω
)
,
0,0
,
0

)
,
N
¯
=
S
¯
+
V
¯
+
I
¯
+
D
¯
+
A
¯
=
Π
μ
.
Following [11], the reproduction number can be established by using the next generation operator approach. The matrices for new infection and transition terms, respectively, given by F and V, are
(3)
F
=
(
β
m
η
1
β
m
η
2
β
m
0
0
0
0
0
0
)
,
V
=
(
σ
1
+
q
+
μ
0
0
-
q
σ
2
+
μ
0
-
σ
1
-
σ
2
μ
+
ψ
)
,
where m
=
(
S
¯
+
(
1
-
ϵ
)
V
¯
)
/
N
¯
=
(
μ
+
ω
+
(
1
-
ϵ
)
ξ
)
/
(
μ
+
ω
+
ξ
).Denote byR
0 the basic reproduction number as
(4)
R
0
=
ρ
(
F
V
-
1
)
=
β
m
σ
1
+
q
+
μ
(
1
+
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
)
,
that is, the spectral radius of the next generation matrix F
V
-
1. Biologically speaking, R
0 is the average number of new secondary infections generated by a single HIV infected individual, introduced into a susceptible population in which some individuals have been vaccinated.


### 3.2. Existence of Backward Bifurcation

Employing the center manifold theory as described in [12], we investigate the existence of backward bifurcation. In order to apply the center manifold theory, we make the following changes of variables:
(5)
S
=
x
1
,
V
=
x
2
,
I
=
x
3
,
D
=
x
4
,
A
=
x
5
,
so that N
=
∑
n
=
1
5
x
n. We now use the vector notation X
=
(
x
1
,
x
2
,
x
3
,
x
4
,
x
5
)
T, where (
·
)
T denotes a matrix transpose. System (1) can then be written as X
˙
=
F
=
(
f
1
,
f
2
,
f
3
,
f
4
,
f
5
)
T, so that
(6)
x
˙
1
(
t
)
=
f
1
=
Π
+
ω
x
2
-
λ
x
1
-
ξ
x
1
-
μ
x
1
,
x
˙
2
(
t
)
=
f
2
=
ξ
x
1
-
(
1
-
ϵ
)
λ
x
2
-
ω
x
2
-
μ
x
2
,
x
˙
3
(
t
)
=
f
3
=
λ
(
x
1
+
(
1
-
ϵ
)
x
2
)
-
(
σ
1
+
q
1
+
h
x
3
+
μ
)
x
3
,
x
˙
4
(
t
)
=
f
4
=
q
x
3
1
+
h
x
3
-
(
σ
2
+
μ
)
x
4
,
x
˙
5
(
t
)
=
f
5
=
σ
1
x
3
+
σ
2
x
4
-
(
ψ
+
μ
)
x
5
,
where λ
=
β
(
(
x
3
+
η
1
x
4
+
η
2
x
5
)
/
N
).Ifβ is taken as the bifurcation parameter and we consider the case R
0
=
1, solving for β gives β
=
β
*; that is,
(7)
β
*
m
(
1
+
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
)
=
σ
1
+
q
+
μ
.
First of all, observe that the eigenvalues of the Jacobian matrix J
(
E
0
) at β
=
β
* [13], that is, J
(
E
0
)
|
β
=
β
*,(8)
J
(
E
0
)
|
β
=
β
*
=
(
-
μ
-
ξ
ω
-
β
*
S
¯
N
¯
-
η
1
β
*
S
¯
N
¯
-
η
2
β
*
S
¯
N
¯
ξ
-
μ
-
ω
-
β
*
(
1
-
ϵ
)
V
¯
N
¯
-
η
1
β
*
(
1
-
ϵ
)
V
¯
N
¯
-
η
2
β
*
(
1
-
ϵ
)
V
¯
N
¯
0
0
β
*
m
-
σ
1
-
q
-
μ
η
1
β
*
m
η
2
β
*
m
0
0
q
-
σ
2
-
μ
0
0
0
σ
1
σ
2
-
μ
-
ψ
)
,are given by
(9)
λ
1
=
-
μ
,
λ
2
=
-
(
μ
+
ω
+
ξ
)
,
λ
3
=
0
.
The other two eigenvalues satisfy the following equation:
(10)
λ
2
+
b
λ
+
c
=
0
,
where
(11)
b
=
(
σ
2
+
μ
)
+
(
μ
+
ψ
)
+
(
σ
1
+
q
+
μ
)
-
β
*
m
,
c
=
(
σ
2
+
μ
)
(
μ
+
ψ
)
+
(
σ
2
+
μ
)
(
σ
1
+
q
+
μ
-
β
*
m
)
+
(
μ
+
ψ
)
(
σ
1
+
q
+
μ
-
β
*
m
)
-
η
1
q
β
*
m
-
σ
1
η
2
β
*
m
.
Substituting (7) into b and c, we find
(12)
b
=
(
σ
2
+
μ
)
+
(
μ
+
ψ
)
+
β
*
m
[
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
]
,
c
=
(
σ
2
+
μ
)
(
μ
+
ψ
)
×
{
1
+
β
*
m
μ
+
ψ
[
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
]
hih
+
β
*
m
σ
2
+
μ
[
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
]
hih
-
β
*
m
η
1
q
+
σ
1
η
2
(
σ
2
+
μ
)
(
μ
+
ψ
)
}
=
(
σ
2
+
μ
)
(
μ
+
ψ
)
×
{
1
+
β
*
m
[
η
2
(
μ
+
ψ
)
2
(
σ
2
q
σ
2
+
μ
+
σ
1
)
hhhhiihhhhh
+
η
1
q
(
σ
2
+
μ
)
2
+
η
2
σ
2
q
(
μ
+
ψ
)
(
σ
2
+
μ
)
2
]
}
.
Clearly, b and c are positive. Equation (10) has two roots with negative real parts. Hence, λ
3
=
0 is a simple zero eigenvalue and all other eigenvalues have negative real parts. The assumptions in [12] are satisfied. Therefore, the center manifold theory can be used to analyze the dynamics of system (1) near β
=
β
*(or, equivalently, R
0
=
1). The Jacobian matrix of system (1) at β
=
β
* has a right eigenvector w, given by w
=
(
w
1
,
w
2
,
w
3
,
w
4
,
w
5
)
T. And it can be computed from the system (
J
(
E
0
)
|
β
=
β
*
)
·
w
=
0; that is,
(13)
0
=
-
(
μ
+
ξ
)
w
1
+
ω
w
2
-
β
*
x
1
¯
N
¯
w
3
-
η
1
β
*
x
1
¯
N
¯
w
4
-
η
2
β
*
x
1
¯
N
¯
w
5
,
0
=
ξ
w
1
-
(
μ
+
ω
)
w
2
-
(
1
-
ϵ
)
β
*
x
2
¯
N
¯
w
3
-
(
1
-
ϵ
)
η
1
β
*
x
2
¯
N
¯
w
4
-
(
1
-
ϵ
)
η
2
β
*
x
2
¯
N
¯
w
5
,
0
=
(
β
*
m
-
σ
1
-
q
-
μ
)
w
3
+
η
1
β
*
m
w
4
+
η
2
β
*
m
w
5
,
0
=
q
w
3
-
(
σ
2
+
μ
)
w
4
,
0
=
σ
1
w
3
+
σ
2
w
4
-
(
μ
+
ψ
)
w
5
;
from (13), we derive the following solutions:
(14)
w
1
=
-
(
μ
+
w
)
2
+
(
1
-
ϵ
)
w
ξ
μ
(
μ
+
w
+
ξ
)
2
β
*
(
w
3
+
η
1
w
4
+
η
2
w
5
)
<
0
,
w
2
=
-
ξ
(
μ
+
w
+
(
1
-
ϵ
)
(
μ
+
ξ
)
)
μ
(
μ
+
w
+
ξ
)
2
β
*
(
w
3
+
η
1
w
4
+
η
2
w
5
)
<
0
,
w
3
=
w
3
>
0
,
w
4
=
q
σ
2
+
μ
w
3
>
0
,
w
5
=
σ
1
w
3
+
σ
2
w
4
μ
+
ψ
>
0
.
The left eigenvector of J
(
E
0
)
|
β
=
β
* is v, denoted by v
=
(
v
1
,
v
2
,
v
3
,
v
4
,
v
5
). And it can be computed from the system (
J
(
E
0
)
|
β
=
β
*
)
T
·
v
=
0; that is,
(15)
0
=
-
(
μ
+
ξ
)
v
1
+
ξ
v
2
,
0
=
ω
v
1
-
(
μ
+
ω
)
v
2
,
0
=
-
β
*
x
1
¯
N
¯
v
1
-
(
1
-
ϵ
)
β
*
x
2
¯
N
¯
v
2
+
(
β
*
m
-
σ
1
-
q
-
μ
)
v
3
+
q
v
4
+
σ
1
v
5
,
0
=
-
η
1
β
*
x
1
¯
N
¯
v
1
-
(
1
-
ϵ
)
η
1
β
*
x
2
¯
N
¯
v
2
+
η
1
β
*
m
v
3
-
(
σ
2
+
μ
)
v
4
+
σ
2
v
5
,
0
=
-
η
2
β
*
x
1
¯
N
¯
v
1
-
(
1
-
ϵ
)
η
2
β
*
x
2
¯
N
¯
v
2
+
η
2
β
*
m
v
3
-
(
μ
+
ψ
)
v
5
,
with the following solutions:
(16)
v
1
=
0
,
v
2
=
0
,
v
3
=
v
3
>
0
,
v
4
=
β
*
m
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
v
3
>
0
,
v
5
=
η
2
β
*
m
μ
+
ψ
v
3
>
0
.The local bifurcation analysis nearβ
=
β
*(R
0
=
1) is then determined by the signs of two associated constants, denoted by a and b, respectively, as
(17)
a
=
∑
k
,
i
,
j
=
1
5
‍
v
k
w
i
w
j
∂
2
f
k
∂
x
i
∂
x
j
(
0,0
)
,
b
=
∑
k
,
i
=
1
5
‍
v
k
w
i
∂
2
f
k
∂
x
i
∂
β
*
(
0,0
)
.
The computations of a and b are done as follows: for system (6) the associated nonzero partial derivatives of F at the disease-free equilibrium are
(18)
∂
2
f
3
∂
x
1
∂
x
3
=
β
*
N
¯
(
1
-
m
)
,
∂
2
f
3
∂
x
1
∂
x
4
=
η
1
β
*
N
¯
(
1
-
m
)
,
∂
2
f
3
∂
x
1
∂
x
5
=
η
2
β
*
N
¯
(
1
-
m
)
,
∂
2
f
3
∂
x
2
∂
x
3
=
β
*
N
¯
(
1
-
ϵ
-
m
)
,
∂
2
f
3
∂
x
2
∂
x
4
=
η
1
β
*
N
¯
(
1
-
ϵ
-
m
)
,
∂
2
f
3
∂
x
2
∂
x
5
=
η
2
β
*
N
¯
(
1
-
ϵ
-
m
)
,
∂
2
f
3
∂
x
3
∂
x
3
=
-
2
β
*
N
¯
m
+
2
q
h
,
∂
2
f
3
∂
x
3
∂
x
4
=
-
β
*
(
1
+
η
1
)
N
¯
m
,
∂
2
f
3
∂
x
3
∂
x
5
=
-
β
*
(
1
+
η
2
)
N
¯
m
,
∂
2
f
3
∂
x
4
∂
x
4
=
-
2
η
1
β
*
N
¯
m
,
∂
2
f
3
∂
x
4
∂
x
5
=
-
(
η
1
+
η
2
)
β
*
N
¯
m
,
∂
2
f
3
∂
x
5
∂
x
5
=
-
2
η
2
β
*
N
¯
m
,
∂
2
f
4
∂
x
3
∂
x
3
=
-
2
q
h
,
∂
2
f
3
∂
x
3
∂
β
*
=
m
,
∂
2
f
3
∂
x
4
∂
β
*
=
η
1
m
,
∂
2
f
3
∂
x
5
∂
β
*
=
η
2
m
.
Substituting (18) into (17), we get
(19)
a
=
2
β
*
N
¯
v
3
{
(
β
*
(
w
3
+
η
1
w
4
+
η
2
w
5
)
(
1
-
ϵ
)
ϵ
ξ
(
μ
+
w
+
ξ
)
2
-
m
(
w
3
+
w
4
+
w
5
)
)
(
w
3
+
η
1
w
4
+
η
2
w
5
)
hhhhhh
×
(
β
*
(
w
3
+
η
1
w
4
+
η
2
w
5
)
(
1
-
ϵ
)
ϵ
ξ
(
μ
+
w
+
ξ
)
2
hhhhhhhhh
-
m
(
w
3
+
w
4
+
w
5
)
(
1
-
ϵ
)
ϵ
ξ
(
μ
+
w
+
ξ
)
2
)
}
+
2
q
h
w
3
2
v
3
(
1
-
β
*
m
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
)
,
b
=
v
3
(
w
3
+
η
1
w
4
+
η
2
w
5
)
m
.From [14, 15], we know that if a
>
0, b
>
0, there exists a backward bifurcation. Since the bifurcation coefficient, b, is always positive, then we establish the following result.Theorem 1.
Ifa
>
0, system (1) exhibits a backward bifurcation when R
0
=
1.Due to existence of backward bifurcation we know that, for positivea, there exists another critical value R
c, which is less than unity, for model (1). Moreover, there is no endemic equilibrium for R
0
<
R
c; there are two distinct endemic equilibria for R
c
<
R
0
<
1, and a unique endemic equilibrium exists for R
0
=
R
c
<
1 or R
0
>
1. Numerical studies will confirm this in the end of this subsection.We now analyze the endemic equilibrium of model (1). The equilibrium of model (1) can be obtained as follows:
(20)
S
*
=
Π
(
(
1
-
ϵ
)
λ
*
+
μ
+
ω
)
(
(
1
-
ϵ
)
λ
*
+
μ
+
ω
)
(
λ
*
+
μ
+
ξ
)
-
ω
ξ
,
V
*
=
Π
ξ
(
(
1
-
ϵ
)
λ
*
+
μ
+
ω
)
(
λ
*
+
μ
+
ξ
)
-
ω
ξ
,
D
*
=
1
σ
2
+
μ
q
I
*
1
+
h
I
*
,
A
*
=
σ
1
μ
+
ψ
I
*
+
σ
2
μ
+
ψ
1
σ
2
+
μ
q
I
*
1
+
h
I
*
,
N
*
=
Π
-
ψ
A
*
μ
,
where
(21)
λ
*
=
λ
(
I
*
)
=
β
I
*
+
η
1
D
*
+
η
2
A
*
N
*
=
β
I
*
N
*
{
1
+
η
1
σ
2
+
μ
q
1
+
h
I
*
hhhhh
+
η
2
(
σ
1
μ
+
ψ
+
σ
2
μ
+
ψ
1
σ
2
+
μ
q
1
+
h
I
*
)
}
.
Substituting (20) into the third equation of system (1), it is easy to derive the following equation:
(22)
f
(
I
)
=
g
(
I
)
,
where
(23)
f
(
I
)
=
(
S
+
(
1
-
ϵ
)
V
)
β
N
×
{
1
+
η
1
σ
2
+
μ
q
1
+
h
I
hhh
+
η
2
(
σ
1
μ
+
ψ
+
σ
2
μ
+
ψ
1
σ
2
+
μ
q
1
+
h
I
)
}
,
(24)
g
(
I
)
=
σ
1
+
μ
+
q
1
+
h
I
.
Clearly, I
*
=
0 is a fixed point, which corresponds to the disease-free equilibrium E
0. For I
=
0, we can obtain
(25)
g
(
0
)
=
σ
1
+
q
+
μ
,
f
(
0
)
=
β
m
{
1
+
η
1
q
σ
2
+
μ
+
σ
2
μ
+
ψ
(
σ
1
+
σ
2
q
σ
2
+
μ
)
}
,
R
0
=
f
(
0
)
g
(
0
)
.
Define
(26)
F
(
I
)
=
f
(
I
)
g
(
I
)
.
From model (1), it can be shown that if I
* is a positive solution of F
(
I
)
=
1, then S
*, V
*, D
*, and A
* are positive. Thus, the equilibrium is biologically relevant. Unfortunately, it is hard to solve the equation F
(
I
)
=
1 analytically; in the following we numerically show that this equation can have two positive roots, which confirms the existence of backward bifurcation. In Figure 1(a), F
(
I
) is plotted versus I for different values of β and all other parameters are fixed. Figure 1(a) shows that an increase in β would lead to curve F
(
I
) becoming tangent to line 1 and defining a critical value F
(
I
*
)
|
β
=
β
c
=
1 and F
′
(
I
*
)
|
β
=
β
c
=
0 hold true. Figure 1(b) shows the occurrence of the backward bifurcation as parameter β varies. We write R
0
(
β
) as the threshold value to indicate β as the bifurcation parameter while all other parameters are fixed. Define [16]
(27)
R
c
=
R
0
(
β
c
)
,
below which the disease-free equilibrium is unique equilibrium.(a) Plot of the functionF
(
I
) with different values of the transmission coefficient β. (b) Backward bifurcation when the transmission coefficient β varies. Other parameters: Π
=
100, μ
=
1
/
45, w
=
1
/
20, ϵ
=
0.4, ξ
=
0.6, σ
1
=
0.7, σ
2
=
0.4, ψ
=
0.5, q
=
1, h
=
1, η
1
=
0.4, and η
2
=
0.7.
(a)
(b)


### 3.3. Stability Analysis of Equilibria

First, we have the following result on the local stability ofE
0.Theorem 2.
The disease-free equilibriumE
0 of system (1) is locally asymptotically stable if R
0
<
1 and unstable otherwise.Proof.
By checking the Jacobian matrix of system (1) evaluated at E
0, we know that the characteristic equation for J
(
E
0
) has two eigenvalues as
(28)
λ
1
=
-
μ
,
λ
2
=
-
(
μ
+
ω
+
ξ
)
,
and the others satisfy the following equation:
(29)
h
(
λ
)
=
β
m
λ
+
σ
1
+
q
+
μ
+
η
1
β
m
q
(
λ
+
σ
1
+
q
+
μ
)
(
λ
+
σ
2
+
μ
)
+
η
2
β
m
q
σ
2
(
λ
+
σ
1
+
q
+
μ
)
(
λ
+
σ
2
+
μ
)
(
λ
+
μ
+
ψ
)
+
η
2
β
m
σ
1
(
λ
+
σ
1
+
q
+
μ
)
(
λ
+
μ
+
ψ
)
=
1
.
If the real parts of the roots of the equationh
(
λ
)
=
1 are nonnegative, that is, ℜ
(
λ
)
≥
0, then [17]
(30)
|
h
(
λ
)
|
≤
h
(
0
)
=
R
0
.
Hence, if R
0
<
1, ∀
λ such that ℜ
(
λ
)
≥
0, then |
h
(
λ
)
|
≤
R
0
<
1, showing that there are no solutions to h
(
λ
)
=
1 with positive real part. Hence, E
0 is locally asymptotically stable if R
0
<
1. This proof is completed.Then, using Lyapunov function we can get global stability ofE
0.Theorem 3.
The disease-free equilibriumE
0 of system (1) is globally asymptotically stable if R
0
<
min
⁡
{
R
c
,
(
σ
1
+
μ
)
m
/
(
σ
1
+
q
+
μ
)
}.Proof.
We note that no endemic equilibrium exists forR
0
<
R
c. Then, E
0 is a unique equilibrium of system (1). We now consider a Lyapunov function:
(31)
V
=
I
+
β
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
D
+
η
2
β
μ
+
ψ
A
.
The time derivative of V is given by
(32)
V
˙
=
I
˙
+
β
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
D
˙
+
η
2
β
μ
+
ψ
A
˙
=
{
β
(
I
+
η
1
D
+
η
2
A
)
S
+
(
1
-
ϵ
)
V
N
-
σ
1
I
-
q
I
1
+
h
I
-
μ
I
}
+
β
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
(
q
I
1
+
h
I
-
σ
2
D
-
μ
D
)
+
η
2
β
μ
+
ψ
(
σ
1
I
+
σ
2
D
-
(
μ
+
ψ
)
A
)
≤
{
β
(
I
+
η
1
D
+
η
2
A
)
-
σ
1
I
-
q
I
1
+
h
I
-
μ
I
}
+
β
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
(
q
I
-
σ
2
D
-
μ
D
)
+
η
2
β
μ
+
ψ
(
σ
1
I
+
σ
2
D
-
(
μ
+
ψ
)
A
)
=
I
{
β
-
σ
1
-
q
1
+
h
I
-
μ
+
β
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
q
hii
+
η
2
β
μ
+
ψ
σ
1
}
≤
I
{
R
0
σ
1
+
q
+
μ
m
-
(
σ
1
+
μ
)
}
=
(
σ
1
+
q
+
μ
)
I
m
{
R
0
-
(
σ
1
+
μ
)
m
σ
1
+
q
+
μ
}
.
Note thatV
˙
≤
0 if R
0
<
(
σ
1
+
μ
)
m
/
(
σ
1
+
q
+
μ
). Furthermore, V
˙
=
0 if and only if I
=
0. Therefore, the largest compact invariant set in Ω: V
˙
=
0, when R
0
<
min
⁡
{
R
c
,
(
σ
1
+
μ
)
m
/
(
σ
1
+
q
+
μ
)
}, is the singleton E
0. Thus, E
0 is globally asymptotically stable if R
0
<
min
⁡
{
R
c
,
(
σ
1
+
μ
)
m
/
(
σ
1
+
q
+
μ
)
}. This completes the proof.


## 3.1. Disease-Free Equilibrium and the Basic Reproduction Number

Model (1) has a disease-free equilibrium (DFE), obtained by setting the right-hand sides of system (1) to zero, represented as
(2)
E
0
:
(
S
¯
,
V
¯
,
I
¯
,
D
¯
,
A
¯
)

=

(

Π
(
μ
+
ω
)
μ
(
μ
+
ξ
+
ω
)
,
Π
ξ
μ
(
μ
+
ξ
+
ω
)
,
0,0
,
0

)
,
N
¯
=
S
¯
+
V
¯
+
I
¯
+
D
¯
+
A
¯
=
Π
μ
.
Following [11], the reproduction number can be established by using the next generation operator approach. The matrices for new infection and transition terms, respectively, given by F and V, are
(3)
F
=
(
β
m
η
1
β
m
η
2
β
m
0
0
0
0
0
0
)
,
V
=
(
σ
1
+
q
+
μ
0
0
-
q
σ
2
+
μ
0
-
σ
1
-
σ
2
μ
+
ψ
)
,
where m
=
(
S
¯
+
(
1
-
ϵ
)
V
¯
)
/
N
¯
=
(
μ
+
ω
+
(
1
-
ϵ
)
ξ
)
/
(
μ
+
ω
+
ξ
).Denote byR
0 the basic reproduction number as
(4)
R
0
=
ρ
(
F
V
-
1
)
=
β
m
σ
1
+
q
+
μ
(
1
+
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
)
,
that is, the spectral radius of the next generation matrix F
V
-
1. Biologically speaking, R
0 is the average number of new secondary infections generated by a single HIV infected individual, introduced into a susceptible population in which some individuals have been vaccinated.


## 3.2. Existence of Backward Bifurcation

Employing the center manifold theory as described in [12], we investigate the existence of backward bifurcation. In order to apply the center manifold theory, we make the following changes of variables:
(5)
S
=
x
1
,
V
=
x
2
,
I
=
x
3
,
D
=
x
4
,
A
=
x
5
,
so that N
=
∑
n
=
1
5
x
n. We now use the vector notation X
=
(
x
1
,
x
2
,
x
3
,
x
4
,
x
5
)
T, where (
·
)
T denotes a matrix transpose. System (1) can then be written as X
˙
=
F
=
(
f
1
,
f
2
,
f
3
,
f
4
,
f
5
)
T, so that
(6)
x
˙
1
(
t
)
=
f
1
=
Π
+
ω
x
2
-
λ
x
1
-
ξ
x
1
-
μ
x
1
,
x
˙
2
(
t
)
=
f
2
=
ξ
x
1
-
(
1
-
ϵ
)
λ
x
2
-
ω
x
2
-
μ
x
2
,
x
˙
3
(
t
)
=
f
3
=
λ
(
x
1
+
(
1
-
ϵ
)
x
2
)
-
(
σ
1
+
q
1
+
h
x
3
+
μ
)
x
3
,
x
˙
4
(
t
)
=
f
4
=
q
x
3
1
+
h
x
3
-
(
σ
2
+
μ
)
x
4
,
x
˙
5
(
t
)
=
f
5
=
σ
1
x
3
+
σ
2
x
4
-
(
ψ
+
μ
)
x
5
,
where λ
=
β
(
(
x
3
+
η
1
x
4
+
η
2
x
5
)
/
N
).Ifβ is taken as the bifurcation parameter and we consider the case R
0
=
1, solving for β gives β
=
β
*; that is,
(7)
β
*
m
(
1
+
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
)
=
σ
1
+
q
+
μ
.
First of all, observe that the eigenvalues of the Jacobian matrix J
(
E
0
) at β
=
β
* [13], that is, J
(
E
0
)
|
β
=
β
*,(8)
J
(
E
0
)
|
β
=
β
*
=
(
-
μ
-
ξ
ω
-
β
*
S
¯
N
¯
-
η
1
β
*
S
¯
N
¯
-
η
2
β
*
S
¯
N
¯
ξ
-
μ
-
ω
-
β
*
(
1
-
ϵ
)
V
¯
N
¯
-
η
1
β
*
(
1
-
ϵ
)
V
¯
N
¯
-
η
2
β
*
(
1
-
ϵ
)
V
¯
N
¯
0
0
β
*
m
-
σ
1
-
q
-
μ
η
1
β
*
m
η
2
β
*
m
0
0
q
-
σ
2
-
μ
0
0
0
σ
1
σ
2
-
μ
-
ψ
)
,are given by
(9)
λ
1
=
-
μ
,
λ
2
=
-
(
μ
+
ω
+
ξ
)
,
λ
3
=
0
.
The other two eigenvalues satisfy the following equation:
(10)
λ
2
+
b
λ
+
c
=
0
,
where
(11)
b
=
(
σ
2
+
μ
)
+
(
μ
+
ψ
)
+
(
σ
1
+
q
+
μ
)
-
β
*
m
,
c
=
(
σ
2
+
μ
)
(
μ
+
ψ
)
+
(
σ
2
+
μ
)
(
σ
1
+
q
+
μ
-
β
*
m
)
+
(
μ
+
ψ
)
(
σ
1
+
q
+
μ
-
β
*
m
)
-
η
1
q
β
*
m
-
σ
1
η
2
β
*
m
.
Substituting (7) into b and c, we find
(12)
b
=
(
σ
2
+
μ
)
+
(
μ
+
ψ
)
+
β
*
m
[
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
]
,
c
=
(
σ
2
+
μ
)
(
μ
+
ψ
)
×
{
1
+
β
*
m
μ
+
ψ
[
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
]
hih
+
β
*
m
σ
2
+
μ
[
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
]
hih
-
β
*
m
η
1
q
+
σ
1
η
2
(
σ
2
+
μ
)
(
μ
+
ψ
)
}
=
(
σ
2
+
μ
)
(
μ
+
ψ
)
×
{
1
+
β
*
m
[
η
2
(
μ
+
ψ
)
2
(
σ
2
q
σ
2
+
μ
+
σ
1
)
hhhhiihhhhh
+
η
1
q
(
σ
2
+
μ
)
2
+
η
2
σ
2
q
(
μ
+
ψ
)
(
σ
2
+
μ
)
2
]
}
.
Clearly, b and c are positive. Equation (10) has two roots with negative real parts. Hence, λ
3
=
0 is a simple zero eigenvalue and all other eigenvalues have negative real parts. The assumptions in [12] are satisfied. Therefore, the center manifold theory can be used to analyze the dynamics of system (1) near β
=
β
*(or, equivalently, R
0
=
1). The Jacobian matrix of system (1) at β
=
β
* has a right eigenvector w, given by w
=
(
w
1
,
w
2
,
w
3
,
w
4
,
w
5
)
T. And it can be computed from the system (
J
(
E
0
)
|
β
=
β
*
)
·
w
=
0; that is,
(13)
0
=
-
(
μ
+
ξ
)
w
1
+
ω
w
2
-
β
*
x
1
¯
N
¯
w
3
-
η
1
β
*
x
1
¯
N
¯
w
4
-
η
2
β
*
x
1
¯
N
¯
w
5
,
0
=
ξ
w
1
-
(
μ
+
ω
)
w
2
-
(
1
-
ϵ
)
β
*
x
2
¯
N
¯
w
3
-
(
1
-
ϵ
)
η
1
β
*
x
2
¯
N
¯
w
4
-
(
1
-
ϵ
)
η
2
β
*
x
2
¯
N
¯
w
5
,
0
=
(
β
*
m
-
σ
1
-
q
-
μ
)
w
3
+
η
1
β
*
m
w
4
+
η
2
β
*
m
w
5
,
0
=
q
w
3
-
(
σ
2
+
μ
)
w
4
,
0
=
σ
1
w
3
+
σ
2
w
4
-
(
μ
+
ψ
)
w
5
;
from (13), we derive the following solutions:
(14)
w
1
=
-
(
μ
+
w
)
2
+
(
1
-
ϵ
)
w
ξ
μ
(
μ
+
w
+
ξ
)
2
β
*
(
w
3
+
η
1
w
4
+
η
2
w
5
)
<
0
,
w
2
=
-
ξ
(
μ
+
w
+
(
1
-
ϵ
)
(
μ
+
ξ
)
)
μ
(
μ
+
w
+
ξ
)
2
β
*
(
w
3
+
η
1
w
4
+
η
2
w
5
)
<
0
,
w
3
=
w
3
>
0
,
w
4
=
q
σ
2
+
μ
w
3
>
0
,
w
5
=
σ
1
w
3
+
σ
2
w
4
μ
+
ψ
>
0
.
The left eigenvector of J
(
E
0
)
|
β
=
β
* is v, denoted by v
=
(
v
1
,
v
2
,
v
3
,
v
4
,
v
5
). And it can be computed from the system (
J
(
E
0
)
|
β
=
β
*
)
T
·
v
=
0; that is,
(15)
0
=
-
(
μ
+
ξ
)
v
1
+
ξ
v
2
,
0
=
ω
v
1
-
(
μ
+
ω
)
v
2
,
0
=
-
β
*
x
1
¯
N
¯
v
1
-
(
1
-
ϵ
)
β
*
x
2
¯
N
¯
v
2
+
(
β
*
m
-
σ
1
-
q
-
μ
)
v
3
+
q
v
4
+
σ
1
v
5
,
0
=
-
η
1
β
*
x
1
¯
N
¯
v
1
-
(
1
-
ϵ
)
η
1
β
*
x
2
¯
N
¯
v
2
+
η
1
β
*
m
v
3
-
(
σ
2
+
μ
)
v
4
+
σ
2
v
5
,
0
=
-
η
2
β
*
x
1
¯
N
¯
v
1
-
(
1
-
ϵ
)
η
2
β
*
x
2
¯
N
¯
v
2
+
η
2
β
*
m
v
3
-
(
μ
+
ψ
)
v
5
,
with the following solutions:
(16)
v
1
=
0
,
v
2
=
0
,
v
3
=
v
3
>
0
,
v
4
=
β
*
m
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
v
3
>
0
,
v
5
=
η
2
β
*
m
μ
+
ψ
v
3
>
0
.The local bifurcation analysis nearβ
=
β
*(R
0
=
1) is then determined by the signs of two associated constants, denoted by a and b, respectively, as
(17)
a
=
∑
k
,
i
,
j
=
1
5
‍
v
k
w
i
w
j
∂
2
f
k
∂
x
i
∂
x
j
(
0,0
)
,
b
=
∑
k
,
i
=
1
5
‍
v
k
w
i
∂
2
f
k
∂
x
i
∂
β
*
(
0,0
)
.
The computations of a and b are done as follows: for system (6) the associated nonzero partial derivatives of F at the disease-free equilibrium are
(18)
∂
2
f
3
∂
x
1
∂
x
3
=
β
*
N
¯
(
1
-
m
)
,
∂
2
f
3
∂
x
1
∂
x
4
=
η
1
β
*
N
¯
(
1
-
m
)
,
∂
2
f
3
∂
x
1
∂
x
5
=
η
2
β
*
N
¯
(
1
-
m
)
,
∂
2
f
3
∂
x
2
∂
x
3
=
β
*
N
¯
(
1
-
ϵ
-
m
)
,
∂
2
f
3
∂
x
2
∂
x
4
=
η
1
β
*
N
¯
(
1
-
ϵ
-
m
)
,
∂
2
f
3
∂
x
2
∂
x
5
=
η
2
β
*
N
¯
(
1
-
ϵ
-
m
)
,
∂
2
f
3
∂
x
3
∂
x
3
=
-
2
β
*
N
¯
m
+
2
q
h
,
∂
2
f
3
∂
x
3
∂
x
4
=
-
β
*
(
1
+
η
1
)
N
¯
m
,
∂
2
f
3
∂
x
3
∂
x
5
=
-
β
*
(
1
+
η
2
)
N
¯
m
,
∂
2
f
3
∂
x
4
∂
x
4
=
-
2
η
1
β
*
N
¯
m
,
∂
2
f
3
∂
x
4
∂
x
5
=
-
(
η
1
+
η
2
)
β
*
N
¯
m
,
∂
2
f
3
∂
x
5
∂
x
5
=
-
2
η
2
β
*
N
¯
m
,
∂
2
f
4
∂
x
3
∂
x
3
=
-
2
q
h
,
∂
2
f
3
∂
x
3
∂
β
*
=
m
,
∂
2
f
3
∂
x
4
∂
β
*
=
η
1
m
,
∂
2
f
3
∂
x
5
∂
β
*
=
η
2
m
.
Substituting (18) into (17), we get
(19)
a
=
2
β
*
N
¯
v
3
{
(
β
*
(
w
3
+
η
1
w
4
+
η
2
w
5
)
(
1
-
ϵ
)
ϵ
ξ
(
μ
+
w
+
ξ
)
2
-
m
(
w
3
+
w
4
+
w
5
)
)
(
w
3
+
η
1
w
4
+
η
2
w
5
)
hhhhhh
×
(
β
*
(
w
3
+
η
1
w
4
+
η
2
w
5
)
(
1
-
ϵ
)
ϵ
ξ
(
μ
+
w
+
ξ
)
2
hhhhhhhhh
-
m
(
w
3
+
w
4
+
w
5
)
(
1
-
ϵ
)
ϵ
ξ
(
μ
+
w
+
ξ
)
2
)
}
+
2
q
h
w
3
2
v
3
(
1
-
β
*
m
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
)
,
b
=
v
3
(
w
3
+
η
1
w
4
+
η
2
w
5
)
m
.From [14, 15], we know that if a
>
0, b
>
0, there exists a backward bifurcation. Since the bifurcation coefficient, b, is always positive, then we establish the following result.Theorem 1.
Ifa
>
0, system (1) exhibits a backward bifurcation when R
0
=
1.Due to existence of backward bifurcation we know that, for positivea, there exists another critical value R
c, which is less than unity, for model (1). Moreover, there is no endemic equilibrium for R
0
<
R
c; there are two distinct endemic equilibria for R
c
<
R
0
<
1, and a unique endemic equilibrium exists for R
0
=
R
c
<
1 or R
0
>
1. Numerical studies will confirm this in the end of this subsection.We now analyze the endemic equilibrium of model (1). The equilibrium of model (1) can be obtained as follows:
(20)
S
*
=
Π
(
(
1
-
ϵ
)
λ
*
+
μ
+
ω
)
(
(
1
-
ϵ
)
λ
*
+
μ
+
ω
)
(
λ
*
+
μ
+
ξ
)
-
ω
ξ
,
V
*
=
Π
ξ
(
(
1
-
ϵ
)
λ
*
+
μ
+
ω
)
(
λ
*
+
μ
+
ξ
)
-
ω
ξ
,
D
*
=
1
σ
2
+
μ
q
I
*
1
+
h
I
*
,
A
*
=
σ
1
μ
+
ψ
I
*
+
σ
2
μ
+
ψ
1
σ
2
+
μ
q
I
*
1
+
h
I
*
,
N
*
=
Π
-
ψ
A
*
μ
,
where
(21)
λ
*
=
λ
(
I
*
)
=
β
I
*
+
η
1
D
*
+
η
2
A
*
N
*
=
β
I
*
N
*
{
1
+
η
1
σ
2
+
μ
q
1
+
h
I
*
hhhhh
+
η
2
(
σ
1
μ
+
ψ
+
σ
2
μ
+
ψ
1
σ
2
+
μ
q
1
+
h
I
*
)
}
.
Substituting (20) into the third equation of system (1), it is easy to derive the following equation:
(22)
f
(
I
)
=
g
(
I
)
,
where
(23)
f
(
I
)
=
(
S
+
(
1
-
ϵ
)
V
)
β
N
×
{
1
+
η
1
σ
2
+
μ
q
1
+
h
I
hhh
+
η
2
(
σ
1
μ
+
ψ
+
σ
2
μ
+
ψ
1
σ
2
+
μ
q
1
+
h
I
)
}
,
(24)
g
(
I
)
=
σ
1
+
μ
+
q
1
+
h
I
.
Clearly, I
*
=
0 is a fixed point, which corresponds to the disease-free equilibrium E
0. For I
=
0, we can obtain
(25)
g
(
0
)
=
σ
1
+
q
+
μ
,
f
(
0
)
=
β
m
{
1
+
η
1
q
σ
2
+
μ
+
σ
2
μ
+
ψ
(
σ
1
+
σ
2
q
σ
2
+
μ
)
}
,
R
0
=
f
(
0
)
g
(
0
)
.
Define
(26)
F
(
I
)
=
f
(
I
)
g
(
I
)
.
From model (1), it can be shown that if I
* is a positive solution of F
(
I
)
=
1, then S
*, V
*, D
*, and A
* are positive. Thus, the equilibrium is biologically relevant. Unfortunately, it is hard to solve the equation F
(
I
)
=
1 analytically; in the following we numerically show that this equation can have two positive roots, which confirms the existence of backward bifurcation. In Figure 1(a), F
(
I
) is plotted versus I for different values of β and all other parameters are fixed. Figure 1(a) shows that an increase in β would lead to curve F
(
I
) becoming tangent to line 1 and defining a critical value F
(
I
*
)
|
β
=
β
c
=
1 and F
′
(
I
*
)
|
β
=
β
c
=
0 hold true. Figure 1(b) shows the occurrence of the backward bifurcation as parameter β varies. We write R
0
(
β
) as the threshold value to indicate β as the bifurcation parameter while all other parameters are fixed. Define [16]
(27)
R
c
=
R
0
(
β
c
)
,
below which the disease-free equilibrium is unique equilibrium.(a) Plot of the functionF
(
I
) with different values of the transmission coefficient β. (b) Backward bifurcation when the transmission coefficient β varies. Other parameters: Π
=
100, μ
=
1
/
45, w
=
1
/
20, ϵ
=
0.4, ξ
=
0.6, σ
1
=
0.7, σ
2
=
0.4, ψ
=
0.5, q
=
1, h
=
1, η
1
=
0.4, and η
2
=
0.7.
(a)
(b)


## 3.3. Stability Analysis of Equilibria

First, we have the following result on the local stability ofE
0.Theorem 2.
The disease-free equilibriumE
0 of system (1) is locally asymptotically stable if R
0
<
1 and unstable otherwise.Proof.
By checking the Jacobian matrix of system (1) evaluated at E
0, we know that the characteristic equation for J
(
E
0
) has two eigenvalues as
(28)
λ
1
=
-
μ
,
λ
2
=
-
(
μ
+
ω
+
ξ
)
,
and the others satisfy the following equation:
(29)
h
(
λ
)
=
β
m
λ
+
σ
1
+
q
+
μ
+
η
1
β
m
q
(
λ
+
σ
1
+
q
+
μ
)
(
λ
+
σ
2
+
μ
)
+
η
2
β
m
q
σ
2
(
λ
+
σ
1
+
q
+
μ
)
(
λ
+
σ
2
+
μ
)
(
λ
+
μ
+
ψ
)
+
η
2
β
m
σ
1
(
λ
+
σ
1
+
q
+
μ
)
(
λ
+
μ
+
ψ
)
=
1
.
If the real parts of the roots of the equationh
(
λ
)
=
1 are nonnegative, that is, ℜ
(
λ
)
≥
0, then [17]
(30)
|
h
(
λ
)
|
≤
h
(
0
)
=
R
0
.
Hence, if R
0
<
1, ∀
λ such that ℜ
(
λ
)
≥
0, then |
h
(
λ
)
|
≤
R
0
<
1, showing that there are no solutions to h
(
λ
)
=
1 with positive real part. Hence, E
0 is locally asymptotically stable if R
0
<
1. This proof is completed.Then, using Lyapunov function we can get global stability ofE
0.Theorem 3.
The disease-free equilibriumE
0 of system (1) is globally asymptotically stable if R
0
<
min
⁡
{
R
c
,
(
σ
1
+
μ
)
m
/
(
σ
1
+
q
+
μ
)
}.Proof.
We note that no endemic equilibrium exists forR
0
<
R
c. Then, E
0 is a unique equilibrium of system (1). We now consider a Lyapunov function:
(31)
V
=
I
+
β
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
D
+
η
2
β
μ
+
ψ
A
.
The time derivative of V is given by
(32)
V
˙
=
I
˙
+
β
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
D
˙
+
η
2
β
μ
+
ψ
A
˙
=
{
β
(
I
+
η
1
D
+
η
2
A
)
S
+
(
1
-
ϵ
)
V
N
-
σ
1
I
-
q
I
1
+
h
I
-
μ
I
}
+
β
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
(
q
I
1
+
h
I
-
σ
2
D
-
μ
D
)
+
η
2
β
μ
+
ψ
(
σ
1
I
+
σ
2
D
-
(
μ
+
ψ
)
A
)
≤
{
β
(
I
+
η
1
D
+
η
2
A
)
-
σ
1
I
-
q
I
1
+
h
I
-
μ
I
}
+
β
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
(
q
I
-
σ
2
D
-
μ
D
)
+
η
2
β
μ
+
ψ
(
σ
1
I
+
σ
2
D
-
(
μ
+
ψ
)
A
)
=
I
{
β
-
σ
1
-
q
1
+
h
I
-
μ
+
β
σ
2
+
μ
(
η
1
+
σ
2
η
2
μ
+
ψ
)
q
hii
+
η
2
β
μ
+
ψ
σ
1
}
≤
I
{
R
0
σ
1
+
q
+
μ
m
-
(
σ
1
+
μ
)
}
=
(
σ
1
+
q
+
μ
)
I
m
{
R
0
-
(
σ
1
+
μ
)
m
σ
1
+
q
+
μ
}
.
Note thatV
˙
≤
0 if R
0
<
(
σ
1
+
μ
)
m
/
(
σ
1
+
q
+
μ
). Furthermore, V
˙
=
0 if and only if I
=
0. Therefore, the largest compact invariant set in Ω: V
˙
=
0, when R
0
<
min
⁡
{
R
c
,
(
σ
1
+
μ
)
m
/
(
σ
1
+
q
+
μ
)
}, is the singleton E
0. Thus, E
0 is globally asymptotically stable if R
0
<
min
⁡
{
R
c
,
(
σ
1
+
μ
)
m
/
(
σ
1
+
q
+
μ
)
}. This completes the proof.


## 4. Persistence of the Model

In this section, we will prove that system (1) is uniformly persistent. First, we present the following definition that is similar to that in [18, 19].Definition 4.
Model (1) is said to be uniformly persistent if there exists a positive constant ɛ
>
0 (independent of initial data) such that every solution with positive initial conditions satisfies liminf
⁡
t
→
∞
S
(
t
)
≥
ɛ, liminf
⁡
t
→
∞
V
(
t
)
≥
ɛ, liminf
⁡
t
→
∞
I
(
t
)
≥
ɛ, liminf
⁡
t
→
∞
D
(
t
)
≥
ɛ, and liminf
⁡
t
→
∞
A
(
t
)
≥
ɛ.Theorem 5.
IfR
0
>
1, system (1) is uniformly persistent; that is, there exists a positive constant ɛ, such that, for all initial values
(33)
(
S
(
0
)
,
V
(
0
)
,
I
(
0
)
,
D
(
0
)
,
A
(
0
)
)
∈
R
+
2
×
I
n
t
(
R
+
3
)
,
the solutions of system (1) satisfy liminf
⁡
t
→
∞
S
(
t
)
≥
ɛ, liminf
⁡
t
→
∞
V
(
t
)
≥
ɛ, liminf
⁡
t
→
∞
I
(
t
)
≥
ɛ, liminf
⁡
t
→
∞
D
(
t
)
≥
ɛ, and liminf
⁡
t
→
∞
A
(
t
)
≥
ɛ.Proof.
Define(34)
X
=
{
(
S
,
V
,
I
,
D
,
A
)
∣
S
≥
0
,
V
≥
0
,
I
≥
0
,
D
≥
0
,
A
≥
0
}
,
X
0
=
{
(
S
,
V
,
I
,
D
,
A
)
∣
S
≥
0
,
V
≥
0
,
I
>
0
,
D
>
0
,
A
>
0
}
,
∂
X
0
=
X
∖
X
0
.
It then suffices to show that system (1) is uniformly persistent with respect to (
X
0
,
∂
X
0
). First, by the form of (1), it is easy to see that both X and X
0 are positively invariant. Clearly, ∂
X
0 is relatively closed in X. Furthermore, system (1) is point dissipative. Set
(35)
M
∂
=
{
(
S
(
0
)
,
V
(
0
)
,
I
(
0
)
,
D
(
0
)
,
A
(
0
)
)
:
(
S
(
t
)
,
V
(
t
)
,
I
(
t
)
,
D
(
t
)
,
A
(
t
)
)
satisfies
(
1
)
,
(
S
(
t
)
,
V
(
t
)
,
I
(
t
)
,
D
(
t
)
,
A
(
t
)
)
∈
∂
X
0
,
∀
t
≥
0
}
.
We now prove that
(36)
M
∂
=
{
(
S
,
V
,
0,0
,
0
)
:
S
≥
0
,
V
≥
0
}
.
Assume that (
S
(
0
)
,
V
(
0
)
,
I
(
0
)
,
D
(
0
)
,
A
(
0
)
)
∈
M
∂. It suffices to show that
(37)
I
(
t
)
=
0
,
D
(
t
)
=
0
,
A
(
t
)
=
0
,
∀
t
≥
0
.
If this is not true, then there exists a t
0
≥
0 such that
(38)
I
(
t
0
)
>
0
,
D
(
t
0
)
=
0
,
A
(
t
0
)
=
0
,
or
(39)
I
(
t
0
)
>
0
,
D
(
t
0
)
>
0
,
A
(
t
0
)
=
0
,
or
(40)
I
(
t
0
)
>
0
,
D
(
t
0
)
=
0
,
A
(
t
0
)
>
0
,
or
(41)
I
(
t
0
)
=
0
,
D
(
t
0
)
>
0
,
A
(
t
0
)
=
0
,
or
(42)
I
(
t
0
)
=
0
,
D
(
t
0
)
=
0
,
A
(
t
0
)
>
0
,
or
(43)
I
(
t
0
)
=
0
,
D
(
t
0
)
>
0
,
A
(
t
0
)
>
0
.
For I
(
t
0
)
>
0, D
(
t
0
)
=
0, and A
(
t
0
)
=
0, we get
(44)
D
′
(
t
0
)
=
q
I
(
t
0
)
1
+
h
I
(
t
0
)
>
0
,
A
′
(
t
0
)
=
σ
1
I
(
t
0
)
>
0
.
It follows that there is an ε
0
>
0 such that D
(
t
)
>
0, A
(
t
)
>
0, for t
0
<
t
<
t
0
+
ε
0. This proves that
(45)
(
S
(
t
)
,
V
(
t
)
,
I
(
t
)
,
D
(
t
)
,
A
(
t
)
)
∉
∂
X
0
for
t
0
<
t
<
t
0
+
ε
0
,
which contradicts the assumption that (
S
(
0
)
,
V
(
0
)
,
I
(
0
)
,
D
(
0
)
,
A
(
0
)
)
∈
M
∂. Similarly, we can obtain contradictions for other cases. This proves that
(46)
M
∂
=
{
(
S
,
V
,
0,0
,
0
)
:
S
≥
0
,
V
≥
0
}
.
Note thatE
0 is globally asymptotically stable in I
n
t
M
∂, and E
0 is an isolated invariant set in X. That is to say, W
s
(
E
0
)
∩
X
0
=
∅. Every orbit in M
∂ converges to E
0, and E
0 is acyclic in M
∂. We claim that W
s
(
E
0
)
∩
X
0
=
∅ for R
0
>
1. If this is false, then we have W
s
(
E
0
)
∩
X
0
≠
∅. The system has a positive solution
(47)
(
S
(
t
)
,
V
(
t
)
,
I
(
t
)
,
D
(
t
)
,
A
(
t
)
)
,
where (
S
(
0
)
,
V
(
0
)
,
I
(
0
)
,
D
(
0
)
,
A
(
0
)
)
∈
X
0. Then
(48)
(
S
(
t
)
,
V
(
t
)
,
I
(
t
)
,
D
(
t
)
,
A
(
t
)
)
⟶
E
0
hhhhhhhhh
as
t
⟶
∞
for
R
0
>
1
.
For R
0
>
1, we can choose an η
>
0 small enough such that R
0
(
1
-
η
)
>
1. Then, when t is sufficiently large, we have
(49)
m
-
η
m
≤
S
(
t
)
+
(
1
-
ϵ
)
V
(
t
)
N
(
t
)
≤
m
+
η
m
,
(50)
I
′
≥
β
m
(
1
-
η
)
(
I
+
η
1
D
+
η
2
A
)
-
σ
1
I
-
q
I
1
+
h
I
-
μ
I
,
D
′
≥
q
I
1
+
h
I
-
σ
2
D
-
μ
D
,
A
′
≥
σ
1
I
+
σ
2
D
-
(
μ
+
ψ
)
A
.
Define
(51)
M
=
(
β
m
(
1
-
η
)
-
σ
1
-
q
-
μ
η
1
β
m
(
1
-
η
)
η
2
β
m
(
1
-
η
)
q
-
σ
2
-
μ
0
σ
1
σ
2
-
μ
-
ψ
)
.
Recall that the stability modulus of an n
×
n matrix M, denoted by s
(
M
), is defined as
(52)
s
(
M
)
=
max
⁡
{
Re
λ
:
λ
is
an
eigenvalue
of
M
}
.
Note thatM is irreducible and has nonnegative off-diagonal elements. It then follows that s
(
M
) is a simple eigenvalue of M with a (componentwise) positive eigenvector. Thus,
(53)
|
λ
I
-
M
|
=
λ
3
+
a
1
λ
2
+
a
2
λ
+
a
3
,
where
(54)
a
1
=
(
σ
1
+
q
+
μ
)
+
(
σ
2
+
μ
)
+
(
μ
+
ψ
)
-
β
m
,
a
2
=
(
σ
1
+
q
+
μ
-
β
m
)
(
σ
2
+
μ
)
+
(
σ
2
+
μ
)
(
μ
+
ψ
)
+
(
σ
1
+
q
+
μ
-
β
m
)
(
μ
+
ψ
)
-
η
1
β
m
q
-
η
2
β
m
σ
1
,
a
3
=
(
σ
1
+
q
+
μ
)
(
σ
2
+
μ
)
(
μ
+
ψ
)
(
1
-
R
0
(
1
-
η
)
)
.
For R
0
(
1
-
η
)
>
1, we obtain a
3
<
0. Thus, s
(
M
) is a simple positive eigenvalue of M with a (componentwise) positive eigenvector. By comparison theorem, we get
(55)
I
(
t
)
⟶
∞
,
D
(
t
)
⟶
∞
,
A
(
t
)
⟶
∞
,
hhhhhhhhhhhhhhhhhhhhhhhhhh
as
t
⟶
∞
,
which contradicts the assumption that
(56)
(
S
(
t
)
,
V
(
t
)
,
I
(
t
)
,
D
(
t
)
,
A
(
t
)
)
⟶
E
0
as
t
⟶
∞
.
This proves that W
s
(
E
0
)
∩
X
0
=
∅ for R
0
>
1. By [20], system (1) is uniformly persistent. Thus, the proof of the theorem is completed.


## 5. Numerical Simulations

### 5.1. Numerical Results

We initially investigate variation inR
0 with different vaccine efficacy, vaccination rate, and diagnosis rate to compare the impact of these intervention measures on HIV transmission. The parameter values in Table 1 are chosen based on HIV/AIDS transmission among IDUs in Yunnan province, China. For simplicity, we choose h
=
1. Differentiating partially R
0 with respect to ξ and ϵ, respectively, we obtain
(57)
∂
R
0
∂
ξ
=
-
β
σ
1
+
q
+
μ
ϵ
(
μ
+
ω
)
(
μ
+
ω
+
ξ
)
2
×
(
1
+
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
)
<
0
,
∂
R
0
∂
ϵ
=
-
β
σ
1
+
q
+
μ
ξ
μ
+
ω
+
ξ
×
(
1
+
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
)
<
0
,
which implies that an increase of vaccination rate and vaccine efficacy leads to the basic reproduction number decline, as shown in Figure 2(a), in which the contour plots of R
0 versus vaccine efficacy ϵ and vaccination rate ξ were plotted. It also shows that the basic reproduction number is more sensitive to vaccine efficacy than vaccination rate. Figure 2(b) shows the contour plot of R
0 with diagnosis rate and vaccination rate, which implies a decrease in R
0 with increasing diagnosis rate q and vaccination rate ξ. Furthermore, when 50
% of HIV individuals are diagnosed, vaccination level of at least 60
% would be needed to achieve R
0
<
1. This suggests that the strategies of diagnosis and vaccination should be stringent enough to reduce R
0.Table 1
Parameter description and values.
Parameters
Description
Estimated values
Source
Π
Recruitment rate
4348
[21]
β
Transmission coefficient
0.304
[26]
w
Per capita waning rate of vaccine
1
/
20
[5]
ξ
Per capita vaccination rate
0.4
Variable
ε
Vaccine efficacy
0.4
Variable
μ
Natural death rate
0.0246
[27]
ψ
Disease-induced death rate
0.7114
[21]
σ
1
Progression rate to AIDS stage for the infection stage
0.0413
[9]
σ
2
Progression rate to AIDS stage for the diagnosed stage
0.116
[9]
q
Diagnosis rate
0.304
[9]
η
1
Modification factor in transmission coefficient of diagnosed HIV-positive individuals
0.491
[9]
η
2
Modification factor in transmission coefficient of AIDS patients
0.1
VariableContour plots ofR
0 versus (a) vaccine efficacy ϵ and vaccination rate ξ and (b) diagnosis rate q and vaccination rate ξ.
(a)
(b)Next, we consider the effect of different transmission rate, vaccination rate, vaccine efficacy, and recruitment rate on transmission of HIV/AIDS. We take the year 2004 as starting time; since then the policy of diagnosis is consistent. In [21], we get that the number of diagnosed HIV-positive individuals and AIDS patients in Yunnan province was 27168 and 1223 in year 2004, respectively. Besides, 22.6
% of these HIV individuals were transmitted by share injectors [22]. Hence, D
(
0
)
=
27168
×
22.6
%
=
6140, A
(
0
)
=
1223
×
22.6
%
=
276. Note that the diagnosis rate is estimated to be 0.304 [9]; then we have I
(
0
)
=
D
(
0
)
/
0.304
=
20197. We have no reliable data on the number of susceptible individuals, that is, number of IDUs in Yunnan province. However, we know that 3.2 million blood samples were tested in Yunnan in [23]. We then assume that in these blood samples the fraction of share injectors is the same as fraction of transmission via share injectors (i.e., 22.6
%). Then the number of susceptible individuals who share injectors is S
(
0
)
+
V
(
0
)
=
3.2
×
10
6
×
22.6
%
=
723200. If the vaccination rate is assumed to be 0.4, then V
(
0
)
=
S
(
0
)
×
0.4. We obtain the initial values
(58)
S
(
0
)
=
516571
,
V
(
0
)
=
206628
,
I
(
0
)
=
20197
,
D
(
0
)
=
6140
,
A
(
0
)
=
276
.
Figure 3 shows the variation in the number of HIV infected individuals with different transmission rates, vaccination rates, vaccine efficacy, and recruitment rates. It follows from Figure 3(a) that decreasing transmission rate could lead to the number of HIV-positive individuals decline. The effect of increasing vaccination rate on HIV transmission is shown in Figure 3(b) and it is seen that the number of HIV-positive individuals becomes much smaller if vaccination rate increases more. Figure 3(c) illustrates that, with increasing vaccine efficacy, the number of HIV-positive individuals decreases. Figure 3(d) shows that if the inflow of susceptible individuals into the community is restricted due to education, the disease spread will slow down.Variation in the number of HIV-positive individuals with different (a) transmission coefficient, (b) vaccination rate, (c) vaccine efficacy, and (d) recruitment rate.
(a)
(b)
(c)
(d)


### 5.2. Sensitivity Analysis

In this section, we use sensitivity analysis method [24] to investigate the impact of various intervention measures on HIV transmission in Yunnan province, China. We hope that these results obtained here could improve the knowledge of the effects of different interventions.Figures4(a) and 4(b) show the comparison of sensitivity coefficients of new cases and prevalence against parameters β, ϵ, ξ, q, and Π, separately. Note that the sensitivity coefficient of new cases and prevalence can be interpreted as the percentage change in the number of new cases and prevalence for 1% decline in the parameters β and Π or 1% increase in ϵ, ξ, and q, respectively [25]. In particular, let function f be new cases or prevalence; the sensitivity coefficients (
SC
) of new cases and prevalence are given by
(59)
SC
=
f
(
perturbed
variables
)

-

f
(
original
variable
)
f
(
original
variable
)

×

100
%
.Sensitivity coefficients of new cases (a) and prevalence (b) onβ, ϵ, ξ, q, and Π over time t. All other parameters are shown in Table 1.
(a)
(b)It follows from Figure4 that a decrease in transmission coefficient β causes new cases and prevalence decline substantially. Besides, an increase in vaccine efficacy ϵ and vaccination rate ξ can lead to a decrease in new cases and prevalence, whereas the change of both diagnosis rate q and recruitment rate Π slightly affects the new cases or prevalence. Thus, new cases and prevalence are sensitive to transmission coefficient, vaccine efficacy, and vaccination rate. Then, reducing the transmission coefficient and increasing the vaccine efficacy and vaccination rate can greatly reduce new cases and prevalence.


## 5.1. Numerical Results

We initially investigate variation inR
0 with different vaccine efficacy, vaccination rate, and diagnosis rate to compare the impact of these intervention measures on HIV transmission. The parameter values in Table 1 are chosen based on HIV/AIDS transmission among IDUs in Yunnan province, China. For simplicity, we choose h
=
1. Differentiating partially R
0 with respect to ξ and ϵ, respectively, we obtain
(57)
∂
R
0
∂
ξ
=
-
β
σ
1
+
q
+
μ
ϵ
(
μ
+
ω
)
(
μ
+
ω
+
ξ
)
2
×
(
1
+
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
)
<
0
,
∂
R
0
∂
ϵ
=
-
β
σ
1
+
q
+
μ
ξ
μ
+
ω
+
ξ
×
(
1
+
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
)
<
0
,
which implies that an increase of vaccination rate and vaccine efficacy leads to the basic reproduction number decline, as shown in Figure 2(a), in which the contour plots of R
0 versus vaccine efficacy ϵ and vaccination rate ξ were plotted. It also shows that the basic reproduction number is more sensitive to vaccine efficacy than vaccination rate. Figure 2(b) shows the contour plot of R
0 with diagnosis rate and vaccination rate, which implies a decrease in R
0 with increasing diagnosis rate q and vaccination rate ξ. Furthermore, when 50
% of HIV individuals are diagnosed, vaccination level of at least 60
% would be needed to achieve R
0
<
1. This suggests that the strategies of diagnosis and vaccination should be stringent enough to reduce R
0.Table 1
Parameter description and values.
Parameters
Description
Estimated values
Source
Π
Recruitment rate
4348
[21]
β
Transmission coefficient
0.304
[26]
w
Per capita waning rate of vaccine
1
/
20
[5]
ξ
Per capita vaccination rate
0.4
Variable
ε
Vaccine efficacy
0.4
Variable
μ
Natural death rate
0.0246
[27]
ψ
Disease-induced death rate
0.7114
[21]
σ
1
Progression rate to AIDS stage for the infection stage
0.0413
[9]
σ
2
Progression rate to AIDS stage for the diagnosed stage
0.116
[9]
q
Diagnosis rate
0.304
[9]
η
1
Modification factor in transmission coefficient of diagnosed HIV-positive individuals
0.491
[9]
η
2
Modification factor in transmission coefficient of AIDS patients
0.1
VariableContour plots ofR
0 versus (a) vaccine efficacy ϵ and vaccination rate ξ and (b) diagnosis rate q and vaccination rate ξ.
(a)
(b)Next, we consider the effect of different transmission rate, vaccination rate, vaccine efficacy, and recruitment rate on transmission of HIV/AIDS. We take the year 2004 as starting time; since then the policy of diagnosis is consistent. In [21], we get that the number of diagnosed HIV-positive individuals and AIDS patients in Yunnan province was 27168 and 1223 in year 2004, respectively. Besides, 22.6
% of these HIV individuals were transmitted by share injectors [22]. Hence, D
(
0
)
=
27168
×
22.6
%
=
6140, A
(
0
)
=
1223
×
22.6
%
=
276. Note that the diagnosis rate is estimated to be 0.304 [9]; then we have I
(
0
)
=
D
(
0
)
/
0.304
=
20197. We have no reliable data on the number of susceptible individuals, that is, number of IDUs in Yunnan province. However, we know that 3.2 million blood samples were tested in Yunnan in [23]. We then assume that in these blood samples the fraction of share injectors is the same as fraction of transmission via share injectors (i.e., 22.6
%). Then the number of susceptible individuals who share injectors is S
(
0
)
+
V
(
0
)
=
3.2
×
10
6
×
22.6
%
=
723200. If the vaccination rate is assumed to be 0.4, then V
(
0
)
=
S
(
0
)
×
0.4. We obtain the initial values
(58)
S
(
0
)
=
516571
,
V
(
0
)
=
206628
,
I
(
0
)
=
20197
,
D
(
0
)
=
6140
,
A
(
0
)
=
276
.
Figure 3 shows the variation in the number of HIV infected individuals with different transmission rates, vaccination rates, vaccine efficacy, and recruitment rates. It follows from Figure 3(a) that decreasing transmission rate could lead to the number of HIV-positive individuals decline. The effect of increasing vaccination rate on HIV transmission is shown in Figure 3(b) and it is seen that the number of HIV-positive individuals becomes much smaller if vaccination rate increases more. Figure 3(c) illustrates that, with increasing vaccine efficacy, the number of HIV-positive individuals decreases. Figure 3(d) shows that if the inflow of susceptible individuals into the community is restricted due to education, the disease spread will slow down.Variation in the number of HIV-positive individuals with different (a) transmission coefficient, (b) vaccination rate, (c) vaccine efficacy, and (d) recruitment rate.
(a)
(b)
(c)
(d)


## 5.2. Sensitivity Analysis

In this section, we use sensitivity analysis method [24] to investigate the impact of various intervention measures on HIV transmission in Yunnan province, China. We hope that these results obtained here could improve the knowledge of the effects of different interventions.Figures4(a) and 4(b) show the comparison of sensitivity coefficients of new cases and prevalence against parameters β, ϵ, ξ, q, and Π, separately. Note that the sensitivity coefficient of new cases and prevalence can be interpreted as the percentage change in the number of new cases and prevalence for 1% decline in the parameters β and Π or 1% increase in ϵ, ξ, and q, respectively [25]. In particular, let function f be new cases or prevalence; the sensitivity coefficients (
SC
) of new cases and prevalence are given by
(59)
SC
=
f
(
perturbed
variables
)

-

f
(
original
variable
)
f
(
original
variable
)

×

100
%
.Sensitivity coefficients of new cases (a) and prevalence (b) onβ, ϵ, ξ, q, and Π over time t. All other parameters are shown in Table 1.
(a)
(b)It follows from Figure4 that a decrease in transmission coefficient β causes new cases and prevalence decline substantially. Besides, an increase in vaccine efficacy ϵ and vaccination rate ξ can lead to a decrease in new cases and prevalence, whereas the change of both diagnosis rate q and recruitment rate Π slightly affects the new cases or prevalence. Thus, new cases and prevalence are sensitive to transmission coefficient, vaccine efficacy, and vaccination rate. Then, reducing the transmission coefficient and increasing the vaccine efficacy and vaccination rate can greatly reduce new cases and prevalence.


## 6. Conclusion

In this paper, we established an epidemic model to investigate effects of saturated diagnosis and vaccination on HIV transmission. It proved that backward bifurcation occurs by employing center manifold theory, which causes the disease-free equilibrium to be locally asymptotically stable instead of globally asymptotically stable forR
0
<
1. Thus, making the basic reproduction number less than unity is not enough to eliminate the HIV infection. We note that R
0
<
1 is equivalent to
(60)
ξ
>
(
(
μ
+
ω
)
[
β
(
1
+
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
)
hhhhhhhhhh
-
(
σ
1
+
q
+
μ
)
(
1
+
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
)
]
)
×
(
(
1
+
η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)
)
(
σ
1
+
q
+
μ
)
hihhhh
-
β
(
1
-
ϵ
)
(

1

+

η
1
q
σ
2
+
μ
+
η
2
μ
+
ψ
(
σ
2
q
σ
2
+
μ
+
σ
1
)

)

)
-
1
:
=
ξ
c
,
which means that only the vaccination rate is greater than ξ
c; HIV infection might be eliminated, depending on initial data. There exists the critical threshold R
c, which cannot be explicitly expressed due to nonlinearity, such that when R
0
<
min
⁡
{
R
c
,
(
σ
1
+
μ
)
m
/
(
σ
1
+
q
+
μ
)
}
<
1, the disease-free equilibrium is globally asymptotically stable. However, if R
0
>
1, the disease uniformly persists.It is interesting to note that if the diagnosis is described linearly backward bifurcation does not happen. This implies that nonlinear diagnosis due to limited medical resources leads to backward bifurcation, and consequently complete elimination of HIV infection becomes difficult. That is, HIV infection might be extinct only by improving integrated interventions, which ensures thatR
0 is less than R
c and (
σ
1
+
μ
)
m
/
(
σ
1
+
q
+
μ
).Since several candidate HIV vaccines are in development, it is useful to study the effectiveness. Moreover, the detection of HIV-positive individuals is limited due to medical resources. We then applied the proposed model with nonlinear diagnosis and vaccination to examine HIV infection among IDUs in Yunnan province, China. Sensitivity analysis shows that new cases and prevalence are sensitive to transmission rate, vaccine efficacy, and vaccination rate, whereas diagnosis rate and recruitment rate slightly affect both of them. Therefore, enlarging vaccination rate, improving vaccine efficacy, and lowering transmission rate by reducing sterile injecting equipment are beneficial to reduce transmission of HIV infection. In order to efficiently reduce HIV transmission, combined intervention strategies are suggested to be implemented simultaneously.Effective antiretroviral therapy (ART) is an important strategy to slow down the progression to AIDS due to great reduction in viral loads and is not included in our model. Note that when HIV infected individuals are diagnosed and CD4 T cell counts decrease to 350 copies/μL, they will accept treatment. We will include treatment strategy to construct HIV/AIDS models to investigate the transmission of HIV/AIDS in the future work and provide policy makers with effective suggestions.



---

*Source: 414383-2014-03-30.xml*