Modelling HIV/AIDS epidemic among men who have sex with men in China.

A compartmental model with antiviral therapy was proposed to identify the important factors that influence HIV infection among gay men in China and suggest some effective control strategies. We proved that the disease will be eradicated if the reproduction number is less than one. Based on the number of annual reported HIV/AIDS among MSM we used the Markov-Chain Monte-Carlo (MCMC) simulation to estimate the unknown parameters. We estimated a mean reproduction number of 3.88 (95% CI: 3.69-4.07). The estimation results showed that there were a higher transmission rate and a lower diagnose rate among MSM than those for another high-risk population. We compared the current treatment policy and immediate therapy once people are diagnosed with HIV, and numerical studies indicated that immediate antiviral therapy would lead to few HIV new infections conditional upon relatively low infectiousness; otherwise the current treatment policy would result in low HIV new infection. Further, increasing treatment coverage rate may lead to decline in HIV new infections and be beneficial to disease control, depending on the infectiousness of the infected individuals with antiviral therapy. The finding suggested that treatment efficacy (directly affecting infectiousness), behavior changes, and interventions greatly affect HIV new infection; strengthening intensity will contribute to the disease control.

1. Introduction

After the detection of the first acquired immunodeficiency syndrome (AIDS) patient, human immunodeficiency virus (HIV) spreads at an alarming rate worldwide. In recent years, though the number of people newly infected with HIV is decreasing, the prevalence of HIV among men who have sex with men (MSM) has increased significantly in China since 2005. In China, the first AIDS patient among MSM was found in 1989 and the proportion of reported cases resulting from homosexual contact in year 2005 to 2011 is 0.4%, 2.5%, 3.4%, 5.9%, 8.6%, 10.8%, and 13.0%, respectively [1]. The prevalence of HIV among MSM is increasing year by year, has reached an average of 5% in large and medium cities, and is greater than 10% in the main cities of the Southwest [1]. Among all people living with HIV, the proportion of people infected by homosexual contact increases from 14.7% in 2009 to 17.4% in 2011 [2]. In mainland China MSM bear a high burden of HIV since the cultural norm may cause them to wish to be hidden, which makes them very likely to be exposed to HIV, and moreover causes much more difficulties to implement intervention measures. At the same time, anal intercourse between men, if unprotected, carries a high risk of HIV transmission [3]. There are studies that showed that HIV infection rate via homosexual contact is much higher than heterosexual contact [4]. So, how to identify the important factors which greatly influence HIV infection among MSM and control the quick increase in HIV epidemic among MSM became urgent in recent years. The primary objective of our study is to understand the HIV epidemic among MSM and to try to suggest some control strategies.

Discovered in 1981, HIV is one of the few things that draw attention from both mathematicians, medical scientists and behavioral scientists. Many models have been proposed in order to predict and control the spread of HIV effectively. A basic SIR model was formulated by Gran et al. [5] in 2008. Considering the feature of chronic disease, the authors generalized the basic model to a Markov model to represent different infection rates in different stages [5]. In such model, HIV/AIDS infected individuals are divided into 13 different stages according to CD4+ level, and at each stage people have a specific infection rate. The idea of dividing the infected group to different stages has been widely used in the study of HIV dynamics [6-10]. In 2008, Zhou et al. [10] proposed a staged discrete model to investigate HIV epidemic in mainland China. We will also utilize the idea in this paper to formulate our mathematical model.

In mainland China, the government initiated a large scale of antiviral therapy since 2003, and free treatment has been expanded to all HIV-positive individuals whose CD4+ T count is less than 350 cells per μL (microliter) of blood from 2007. In recent years, antiviral therapy has become a measure to control the epidemic since there are studies demonstrating that antiviral therapy can extend the life expectancy of HIV/AIDS infected individuals and reduce their infection rates [4, 13, 18, 19]. Granich et al. [12] put forward a new therapy, starting antiviral therapy as soon as people are diagnosed with HIV. And he proved that this therapy can reduce the infection rate and death rate greatly. Note that lower infection rate is beneficial to the disease control, but extended life span increases the risk to infect others. For this high-risk group, whether immediate antiviral therapy once people are diagnosed with HIV is beneficial to disease control or not is controversial since early treatment will inevitably result in early occurrence of drug resistance due to bad adherence and no additional drugs are available [20]. So, to what extent antiviral therapy will be implemented or when to initiate antiviral therapy is unclear. All these fall within the scope of this study.

This paper is divided into 6 sections. In Section 2, a mathematical model with antiviral therapy is formulated following a modeling approach for stratification of the population according to the clinical progression of the disease and epidemiological status of the individuals [6]. In Section 3, the unknown parameters and initial values involved in the model are estimated. In Section 4, the main results of this paper are given which include the threshold dynamics, prediction for the disease, the effects of antiviral therapy, and the effects of intervention measures. In Section 5, some sensitivity analyses are conducted. At last, discussions are given in Section 6.

2. Mathematical Model

Our model is formulated based on the key epidemiological properties of HIV/AIDS and some implemented public heath interventions such as condom use and antiviral therapy. The underlying structure of the model comprises classes of individuals among MSM who are high-risk susceptibles (S), HIV infected but not yet diagnosed (I), diagnosed HIV-positive individuals who have not yet progressed to AIDS stage, and those with AIDS clinical symptoms (A). Further, according to the CD4+ T cell counts in the blood we divide the diagnosed HIV-positive individuals without receiving treatment into two different stages: the HIV-positive individuals with l > 350 (I 1), individuals with 200 < l < 350 (I 2), where l denotes the CD4+ T level in the blood (i.e., the number of CD4+ T cells per microliter of blood). HIV-positive individuals in group I 1 and I 2 will enter into T 1 and T 2 if they receive antiviral therapy, respectively, and finally progress to the AIDS stage (T ). The flow diagram is described in Figure 1. Let A be AIDS patients who are diagnosed after onset of AIDS. Since all AIDS patients are given treatment at the present stage in China, here we do not distinguish between AIDS patients in A and those in T . The high-risk population size is represented by N; that is, N = S + I + I 1 + I 2 + A + T 1 + T 2 + T ,

Figure 1

The flow diagram of model with antiviral therapy.

We assume that people enter into the susceptible class at a rate U, exit high-risk population at a constant rate μ. Each one has a natural death rate (d). It is assumed that the transmission probability for the undiagnosed HIV-positive individuals per high-risk behavior (i.e., male-male sex) is a constant, denoted by β. c represents the contact rate per year, and π represents the protection rate by interventions such as condom use and antiviral therapy. Since HIV/AIDS infected individuals in different stages may have different transmission probabilities per high-risk behavior due to different CD4 cell count, antiviral therapy, behavior changes, and so forth, we then introduce a modification factor for each infectious class, denoted by ε , ε , ε , ε , and ε , respectively. Thus, susceptible people become HIV infected at a rate S(I + ε I 1 + ε I 2 + ε T 1 + ε T 2 + ε (A + T ))/N. ρ 1 and ρ 2 are the proportions of HIV-positive individuals with l > 350 and 200 < l < 350 when diagnosed, respectively. Let ρ = ρ 1 + ρ 2, and then 1 − ρ is the proportion of individuals who have already progressed to AIDS stage when diagnosed. Constants ξ 1 and ξ 2 represent the progression rates from I 1 to I 2 and from I 2 to A, respectively. η 1 and η 2 denote the progression rates from T 1 to T 2 and from T 2 to T , respectively. α   (i = I, T, A) and τ   (i = 1,2) represent the disease related death rate and antiviral therapy coverage rate for each stage. The model equations are given in (1), where λ(t) = βc(1 − π)S(I + ε I 1 + ε I 2 + ε T 1 + ε T 2 + ε (A + T ))/N. The model is illustrated in Figure 1, and definitions of parameters are given in Table 2. We can verify that any solution of system (1) with nonnegative initial value is nonnegative. And it is easy to get that G = {(S, I, I 1, I 2, A, T 1, T 2, T ) ∈ R + 8 : S + I + I 1 + I 2 + A + T 1 + T 2 + T ≤ U/d} is positively invariant.

Table 2

Parameters and initial values.

Parameters	Definition	Value	Std	Source
β	Transmission probability of HIV per high-risk behavior	—	—	—
c	Contact rate per year	—	—	—
π	Protection rate by intervention measures (condom use)	—	—	—
β 0	Transmission coefficient, β 0 = βc(1 − π)	0.8578	0.013	MCMC
ϵ I1	Modification factor for HIV infected individuals with l > 350	0.7224	—	[5]
ϵ I2	Modification factor for HIV infected individuals with 200 < l < 350	0.4177	—	[5]
ϵ T1	Modification factor for HIV infected individuals with l > 350 and receiving antiviral therapy	0.3507	—	—
ϵ T2	Modification factor for HIV infected individuals with 200 < l < 350 and receiving antiviral therapy	0.3507	0.028	MCMC
ϵ A	Modification factor for AIDS patients with antiviral therapy	0.2495	0.028	MCMC
U	Recruitment rate of susceptible	402450	20764	MCMC
d	Natural death rate	0.0149	—	[11]
μ	Exit rate of susceptible	0.025	—	—
δ	Diagnose rate	0.0799	0.020	MCMC
ρ	Proportion of diagnosed HIV-positive individuals ρ = ρ 1 + ρ 2	0.8820	0.006	MCMC
ρ 1	Proportion of diagnosed HIV-positive individuals with l > 350	0.5879	—	[10]
ρ 2	Proportion of diagnosed HIV-positive individuals with 200 < l < 350	0.2939	—	[10]
ξ 1	Progression rate from I 1 to I 2	1/6	—	[10]
ξ 2	Progression rate from I 2 to A	1/3	—	[10]
η 1	Progression rate from T 1 to T 2	1/12	—	[12]
η 2	Progression rate from T 2 to T A	1/6	—	[12]
τ 1	Antiviral therapy coverage rate for HIV-positive individuals with l > 350	0	—	—
τ 2	Antiviral therapy coverage rate for HIV-positive individuals with 200 < l < 350	0.2	—	[13, 14]
α I	Disease-related death rate for HIV infected individuals without receiving antiviral therapy	0.172	—	[15]
α T	Disease-related death rate for HIV infected individuals with antiviral therapy	0.06	—	[15, 16]
α A	Disease-related death for AIDS patients with antiviral therapy	0.136	—	[16, 17]
S(0)	Initial value of S (high-risk susceptible)	763240	87285	MCMC
I(0)	Initial value of I (individuals living with undiagnosed HIV)	5988	1274	MCMC
I 1(0)	Initial value of I 1 (HIV-positive individuals with l > 350)	99	—	Database
I 2(0)	Initial value of I 2 (HIV-positive individuals with 200 < l < 350)	49	—	Database
A(0)	Initial value of A (AIDS patients diagnosed at onset of AIDS symptoms)	53	—	Database
T 1(0)	Initial value of T 1 (HIV-positive individuals with l > 350 and receiving therapy)	0	—	Database
I 2(0)	Initial value of T 2 (HIV-positive individuals with 200 < l < 350 and receiving therapy)	0	—	Database
T A(0)	Initial value of T A (AIDS patients diagnosed before onset of AIDS symptoms)	0	—	Database
R 0	The basic reproduction number	3.8840	0.097	Calculated

3. Data and Parameter Estimation

3.1. Data

Based on the current surveillance system we can get the number of annual reported HIV-positive cases and AIDS patients among MSM from year 1985 to 2009 in mainland China. The first HIV infected case among MSM in China was reported in 1989. However, the disease has spread very quickly among this population. The proportion of reported cases resulting from homosexual contact increases as follows: 2.5% (2006), 3.4% (2007), 5.9% (2008), and 8.6% (2009) [1]. Sentinel surveillance results have shown that the HIV-positive rates among MSM population in different regions have been consistently greater than 1% and are increasing year by year. Homosexual transmission is becoming one of the most important drivers of the AIDS epidemic. We choose year 2005 as a starting point since the consistent surveillance and testing policy has been implemented in mainland China since then [15].

3.2. Parameter Values Obtained from Data or the Literature

The natural death rate was estimated to be d = 0.0149 [11]. Generally speaking, if a person is a homoeroticism he tends to be a homoeroticism, through all his life. So we roughly choose μ = 1/40  year−1. Zhou et al. [10] estimated the progression rate from I 1 to I 2 is ξ 1 = 1/6 and the progression rate from I 2 to A is ξ 2 = 1/3. Gran et al. [5] estimated that the scale of progression rate reduced by antiviral therapy is 1/2. So we choose η 1 = 1/2ξ 1 = 1/12 and η 2 = 1/2ξ 2 = 1/6. Because ξ 1 = 1/2∗ξ 2, we can suppose that one-third of the HIV-positive individuals are in stage I 2 and two-thirds are in stage I 1 when diagnosed; that is, ρ 1 = 2ρ/3 and ρ 2 = ρ/3. Xiao et al. [15] estimated the additional death rate for the diagnosed HIV-positive individuals to be 0.172; then we choose α = 0.172. The disease related death rate for AIDS patients without receiving treatment is estimated to be 0.393 by Zhang et al. [17]. The ministry of health, China [16], estimated that the disease related death rate can be decreased by 65.3% after antiviral therapy, so the disease related death rate for AIDS patients with antiviral therapy is α = 0.393∗(1 − 65.3%) = 0.136. The disease-related death rate for HIV-positive individuals with antiviral therapy is α = 0.172∗(1 − 65.3%) = 0.06.

Recently, antiviral therapy is not provided to HIV-positive individuals with l > 350, we then set τ 1 = 0. According to the treatment policy of China, the HIV-positive individuals with l > 350; are given treatment since early 2007. Lou et al. [13] and Xu et al. [14] estimated the treatment coverage rate for HIV-positive individuals with 200 < l < 350 nowadays to be 0.2. Then we choose time-dependent function τ 2 = 0.2 from year 2007 to 2009 and τ 2 = 0 from year 2005 to 2006. From the database we can get some initial values for the system: I 1(0) = 99, I 2(0) = 49, and A(0) = 53. Considering that there are no HIV-positive individuals receiving antiviral therapy in year 2005, thus we have T 1(0) = T 2(0) = T (0) = 0.

In order to determine the relative infectiousness for each class we follow the principle employed by Gran et al. [5]. They divided the diagnosed HIV-positive individuals into 5 stages based on the CD4+ T cell counts in the blood (above 500, 350–499, 200–349, below 200 copies, and AIDS stage). The numbers of new infections caused by one individual, in each stage, per unit time in a totally susceptible population are 0.031, 0.025, 0.017, 0.013, and 0 for these 5 stages, respectively. Note that we combine the first two stages as one stage in our model. We calculate ε as follows: the number of new infections in stage I 1 is ((1/0.0455)/(1/0.0455 + 1/0.125)) × 0.031 + ((1/0.125)/(1/0.0455 + 1/0.125)) × 0.025 = 0.0294, where 0.0455 and 0.125 are the progression rates from stage I 1 to stage I 2 and from stage I 2 to stage I 3. By similar method we get that the numbers of new infections in I and I 2 are 0.0407 and 0.017, respectively. Then, we have 1 : ε : ε = 0.0407 : 0.0294 : 0.017, which implies that ε = 0.7224, ε = 0.4177. Other unknown parameters are estimated on the basis of the real data.

3.3. MCMC Procedure

We utilize the Markov-Chain Monte-Carlo (MCMC) simulation to estimate the main parameters, and initial values of model (1). Metropolis-Hastings (M-H) algorithm [21-26] are employed for MCMC simulation. The algorithm runs for 210000 iterations with a burn-in of 10000 iterations. Using the rest 200000 samples, we obtain our estimates [27]. The estimation procedure is carried out in two steps. Firstly, in order to reduce the number of parameters that need to be estimated to the minimum, we use the method proposed by Tang et al. [28]. Note that the size of high-risk population is very large compared with HIV/AIDS infected population; S/N approximately equals 1. Assuming S = N, we can get a reduced model ((A.1)), given in Appendix A. By fitting the reduced model ((A.1)) to the data on annual reported HIV/AIDS cases from 2005 to 2009, we get the estimates for the transmission coefficient β 0, proportion of diagnosed individuals who are HIV-positive ρ, diagnose rate δ, modification factors ε , ε , and initial value of undiagnosed infected individuals I(0). Here we assume ε = ε since there are no significant difference between the infectiousness of people in these two groups. Using the simulation results of these parameters we can get the Markov Chain of basic reproduction number R 0 and its variation. Secondly, based on the estimation results we repeat the same procedure by fitting model (1) to our real data to estimate the rest of the unknown parameters S(0) and U. See Appendix A for details about parameter estimation.

4. Main Results

4.1. Threshold Dynamics

The basic reproduction number R 0, the average number of secondary cases generated by a single primary case in a fully susceptible population during its average infectious period [29, 30], is a threshold parameter for the infectious disease and can help determine whether an infectious disease will spread through a population. Following the next-generation matrix method [29] we can get R 0 for system (1) (see Appendix B for details). For convenience, we let β 0 = βc(1 − π), ω 1 = d + α + δ, ω 2 = d + α + τ 1 + ξ 1, ω 3 = d + α + τ 2 + ξ 2, ω 4 = d + α , ω 5 = d + α + η 1 and ω 6 = d + α + η 2. The expression of R 0 is given as follows:

As the model shows, the primary infected individual will stay at the compartments I, I 1, I 2, A, T 1, T 2, and T (if experiencing the stage) for mean times of 1/ω 1, 1/ω 2, 1/ω 3, 1/ω 4, 1/ω 5, 1/ω 6, and 1/ω 4, respectively. The probabilities of the primary patient to experience each stage are 1, ρ 1 δ/ω 1, (ρ 1 δξ 1/ω 1 ω 2) + (ρ 2 δ/ω 1), (ρ 1 δξ 1 ξ 2/ω 1 ω 2 ω 3) + (ρ 2 δξ 2/ω 1 ω 3) + ((1 − ρ)δ/ω 1), ρ 1 δτ 1/ω 1 ω 2, (ρ 1 δτ 1 η 1/ω 1 ω 2 ω 5) + (ρ 1 δτ 2 ξ 1/ω 1 ω 2 ω 3) + (τ 2 ρ 2 δ/ω 1 ω 3), and (ρ 1 δτ 1 η 1 η 2/ω 1 ω 2 ω 5 ω 6)+(ρ 1 δξ 1 τ 2 η 2/ω 1 ω 2 ω 3 ω 6) + (ρ 2 δτ 2 η 2/ω 1 ω 3 ω 6). And the transmission coefficients at each stage are β 0, β 0 ε , β 0 ε , β 0 ε , β 0 ε , β 0 ε , and β 0 ε . It is easy to see that β 0/ω 1, β 0 ε /ω 2, β 0 ε /ω 3, β 0 ε /ω 4, β 0 ε /ω 5, β 0 ε /ω 6, and β 0 ε /ω 4 represent the numbers of individuals the primary case infects at each stage, respectively. Using the total probability formula we get the expression of R 0.

Meanwhile, system (1) has a disease-free equilibrium (DFE) E 0 = (S 0, 0,0, 0,0, 0,0, 0), where S 0 = U/(d + μ). From [29] we can get that when R 0 < 1, the disease-free equilibrium of system (1) is locally stable but unstable when R 0 > 1. We can also prove that E 0 is globally attractive when R 0 < 1 (see details in Appendix C.1). Thus, we get the following theorem.

Theorem 1

When R 0 < 1, the disease-free equilibrium of system (1) is globally asymptotically stable, but unstable when R 0 > 1.

Besides the disease-free equilibrium the system has an endemic equilibrium E + when R 0 > 1. It is not difficult to get the expression of E + = (S*, I*, I 1*, I 2*, A*, T 1*, T 2*, T *), where S* = (U − ω 1 I*)/(d + μ), I 1* = e 1 I*, I 2* = e 2 I*, A* = e 3 I*, T 1* = e 4 I*, T 2* = e 5 I*, T * = e 6 I*, I* = U(R 0 − 1)/(ω 1 R 0 +(d + μ)(1 + e 1 + e 2 + e 3 + e 4 + e 5 + e 6 − (ω 1/(d + μ)))), e 1 = ρ 1 δ/ω 2, e 2 = (ρ 2 δ/ω 3) + (ρ 2 δξ 1/ω 2 ω 3), e 3 = ((1 − ρ)δ/ω 4) + (ρ 2 δξ 2/ω 3 ω 4) + (ρ 1 δξ 1 ξ 2/ω 2 ω 3 ω 4), e 4 = ρ 1 δ(1 + σ)τ 1/ω 2 ω 5, e 5 = (ρ 2 δ(1 + σ)τ 2/ω 3 ω 6) + (ρ 1 δξ 1(1 + σ)τ 2/ω 2 ω 3 ω 6) + (ρ 1 δ(1 + σ)τ 1 η 1/ω 2 ω 5 ω 6), and e 6 = (ρ 2 δ(1 + σ)τ 2 η 2/ω 3 ω 6 ω 4) +(ρ 1 δξ 1(1 + σ)τ 2 η 2/ω 2 ω 3 ω 6 ω 4) + (ρ 1 δ(1 + σ)τ 1 η 1 η 2/ω 2 ω 5 ω 6 ω 4).

Moreover, according to the persistence theorem developed by Smith and Zhao in [31] we can prove the system is uniformly persistent if R 0 > 1. To this end we verify that invariant sets in the boundary of the feasible region are not attractors (see details in Appendix C.2.) Hence we have the following theorem.

Theorem 2

System (1) is uniformly persistent when R 0 > 1, in the sense that there is a positive number δ 0 such that for all initial values (S 0, I 0, I 1, I 2, A 0, T 1, T 2, T ) ∈ R + × Int(R + 7), the solution of system (1) satisfies

4.2. Prediction for the Disease

Based on the number of individuals living with HIV (not AIDS) or AIDS among MSM from year 2005 to year 2009, we estimate mean values of parameters and their standard deviations which are listed in Table 2. We also derive the goodness of fit with data together with the uncertainties (shown in Figure 2). The areas from light to dark mean the 50%, 90%, 95%, and 99% limits of the posterior uncertainty due to model parameters. In particular, we get the mean value of the basic reproduction number to be 3.8840 with a estimation error 0.097, and the 95% confidence interval is [3.69,4.07]. Figure 3 describes the Markov Chain for R 0, which has a good convergency. From the estimation results we find that the value of transmission coefficient β 0 is larger than the estimation for heterosexual transmission by Xu et al. [14] and that for the general high-risk population without considering transmission routs by Xiao et al. [15]. This implies that, on the one hand, HIV transmission probability among MSM are higher than that for other high-risk groups [4]; on the other hand, the condom use rate or other interventions among MSM are a bit lower [32]. What is more, the diagnose rate among MSM estimated is only 0.08, which is lower than that for other high-risk groups. That is almost due to the physiological specificity of MSM and the cultural norm of China. A study for Beijing was given in 2006, which pointed out that testing and awareness of HIV infection among MSM are very lacking [33].

Figure 2

Plots of data fitted results. (a) The number of annual reported AIDS patients. (b) The number of annual reported HIV-positive individuals. Squares denote the real data. Areas from light to dark mean the 50%, 90%, 95%, and 99% predictive probability limits due to parameter uncertainties.

Figure 3

MCMC plots for R 0.

We suppose that the treatment coverage rate will not change and the current treatment policy will remain the same (i.e., antiviral therapy starts when CD4+ T cell counts are less than 350 per μL of blood). Our estimation shows that the total number of HIV/AIDS individuals among MSM in 2015 will reach 1.46 × 106, with 1.41 × 106 HIV-positive individuals and 4.41 × 104 AIDS patients if the current surveillance, testing, and interventions are unchanged (shown in Figure 4(a)). In particular, our estimation gives that there are 6596 HIV/AIDS infected individuals among MSM receiving antiviral therapy in 2011, which takes up to 55% of the diagnosed individuals among MSM with CD4+ T count less than 350 cells per μL, whereas, estimation from Ministry of Health, people's republic of China stated that the antiviral therapy coverage rate has reached 73.5% in 2011 [2]. This difference agrees with the conclusion obtained by Tong that the treatment coverage rate for HIV/AIDS cases infected sexually is lower than that for HIV/AIDS individuals infected by other routes [34].

Figure 4

Prediction of HIV/AIDS among MSM in China from 2005 to 2020 and the uncertainties of the model. Areas from light to dark mean the 50%, 90%, 95%, and 99% predictive probability limits due to parameter uncertainties. (a) Total HIV/AIDS cases. (b) Total HIV-positive cases. (c) Total AIDS cases. (d) Total HIV/AIDS cases with antiviral therapy. Parameters and initial values used are shown in Table 2. Antiviral therapy started when CD4+ counts are less than 350 cells per μL.

4.3. Effects of Antiviral Therapy

Antiviral therapy can effectively decrease the viral load in the blood of HIV/AIDS infected individuals, thereby reduceing their infectivity. However, the decrease in viral load consequently alleviates the symptom of HIV/AIDS infected individuals; then they may become active or increase their high-risk behaviors. There are some studies showed that HIV infected individuals may increase their high-risk behaviors since they believe that antiviral therapy can make them more healthy [4, 35–37]. Moreover, antiviral therapy can prolong the lifespan of HIV/AIDS infected individuals, which consequently makes them have more opportunities to infect others. Some researchers supported the conclusion that immediate antiviral therapy once people are diagnosed with HIV is more effective for reducing the infection rate and death rate, and finally controlling the infection [12], whereas other researchers worried that early initiating treatment inevitably induces early occurrence of drug resistance due to poor drug adherence and side effect and consequently results in the decline in treatment efficacy hence it is harmful to control HIV infection [20]. So when antiviral therapy should be started is still controversial and needs a further study. In order to investigate this issue we consider the following two situations and compare the basic reproduction numbers of these two situations.

Situation 1. Start antiviral therapy when CD4+ T counts are less than 350 cells per μL of blood (keep the current treatment policy). Then we have τ 1 = η 1 = 0. Let ω 2* = d + α + ξ 1; we get the expression of basic reproduction number in this situation, which is denoted by R 0*; then

Situation 2. Start antiviral therapy immediately once people are diagnosed with HIV. In this situation, the expression of basic reproduction number R 0 is given in Section 4.1.

Using the parameter values listed in Table 2 we get that R 0 = 3.93, R 0* = 3.88. This implies that immediate antiviral therapy once diagnosed may cause more new infections and is not beneficial to the disease control if parameter values are chosen as those listed in Table 2. In such scenarios, the effect of expanded lifespan may overweight the effect of reduced infectiousness caused by early treatment.

It is interesting to note that the value of R 0 is greater than R 0* based on the parameter values shown in Table 2. However, immediate antiviral therapy may not always lead to a greater value of the reproduction number. In fact, the relationship between R 0 and R 0* depends on many factors such as treatment coverage, infectiousness, and so forth. In the following we investigate the relation between R 0 and R 0*: In order to determine the sign of R 0 − R 0*, we turn to consider It is obvious that sign⁡(R 0 − R 0*) = sign⁡(ω 1 ω 2(R 0 − R 0*)/β 0 ρ 1 δτ 1), where sign⁡(a) denotes the sign function and indicates the sign of a. So if ε and η 1 satisfy the following: we have R 0 − R 0* < 0, which means that immediate antiviral therapy once people are diagnosed with HIV will result in a lower value of basic reproduction number and hence lead to fewer new infections.

From (7) we obtain that the difference between R 0 and R 0* is determined by many parameter values. For simplicity, we introduce a factor v which varies from 0 to 1 to let the modification factors vary simultaneously. To this end, let v multiply ε , ε , and ε ; we then analyze the value of R 0* and R 0 as a function of v. Thus, v = 0 denotes that the infected individuals are not infectious anymore and v = 1 means the modification factors ε , ε , and ε take values as those listed in Table 2. In order to compare the two situations (antiviral therapy started when CD4+ count less than 350 cells per μL, antiviral therapy started immediately once people are diagnosed with HIV) we draw R 0* and R 0 in Figure 5(a). It shows that there is a critical value v* = 0.81 such that R 0 < R 0* for v < v* and R 0 > R 0* for v > v*. That is to say, if we can reduce the modification factors to 81%, immediate antiviral therapy will be more beneficial in terms of inducing a less reproduction number. Summing up the above, if the modification factors for people with antiviral therapy are relatively small, the earlier to start antiviral therapy the better; in contrast, if the modification factors are relatively large, postponing antiviral therapy is better in terms of causing relatively few new infections.

Figure 5

(a) Plots of R 0* and R 0 against factor v. (b) Contour plot of R 0 varies the factor v and treatment coverage rate, where τ 1 = τ 2 = τ. Parameters used are shown in Table 2.

Should We Increase the Treatment Coverage Rate? In the following we investigate whether increasing the treatment coverage rate is beneficial or harmful to control HIV epidemic. To this end, we examine the variation in R 0 with treatment coverage rate. We simply suppose τ 1 = τ 2 = τ. Solving the partial derivative of R 0 with respect to τ we have where

We can verify that A > 0,  B > 0 (see Appendix D for details about the proof of B > 0). Let ∂R 0/∂τ = 0; we have

Thus, there exists a critical value for a combination of parameters η 1, η 2, α , and so forth such that ∂R 0*/∂τ 2 < 0 for and ∂R 0*/∂τ 2 > 0 for . Note that A and B are both positive, which implies that if the modification factors for HIV-positive individuals are relatively small, then a large treatment coverage rate will lead to decline in new infections and hence is beneficial to disease control; otherwise increasing treatment coverage rate would result in an increase in new infections. That is to say, given unchanged behavior of HIV-positive individuals, the infectiousness due to antiviral therapy plays a vital role in determining HIV new infections. In particular, when infectiousness due to antiviral therapy is below a critical value, increasing treatment coverage rate is effective to reduce the new infections. Whereas infectiousness is greater than a critical value, increasing treatment coverage rate induces more new infections, which implies the more therapy for the HIV infected individuals, the more new infections produced (see Figure 5(b) for details). It follows from Figure 5(b) that R 0 decreases as treatment uptake rate τ increases if v is less than a critical value (around 0.6), which implies that increasing treatment coverage is beneficial to disease control in the case of better treatment efficacy; otherwise, increasing treatment is harmful to it.

4.4. Effects of Intervention Measures

In recent years many intervention measures have been implemented to control the quick increase of HIV epidemic among MSM. The intervention measures include (1) strengthening education to the high-risk population, which decreases the constant recruitment rate U, the contact rate per year c and increases the condom use rate; (2) increasing surveillance and testing, which results in an increase in the diagnose rate δ. Note that both decreasing in contact rate c and increasing in condom use rate π can lead to the decline of transmission coefficient β 0.

To address the impact of each intervention measure on HIV infection among MSM, we investigate variation in number of HIV/AIDS infected individuals with varying transmission coefficients, treatment uptake rate. Figure 6 shows that strengthening education to high-risk population (i.e., smaller β 0) and increasing surveillance and testing (i.e., larger δ) slow down the spread of HIV significantly. In particular, if β 0 decreases by 30% and δ increases by 25%, the total number of HIV/AIDS infected individuals among MSM will decrease to 2.2959 × 105, which is decreased by 84.2%. Similarly, we could increase the diagnose rate δ only to examine the effect of improving HIV diagnose and testing strategy on HIV transmission among MSMs. In particular, if the diagnose rate δ increases to 0.3, 0.5, and 0.8, the total number of HIV/AIDS infected individuals will decrease by 32.74%, 52.14%, and 68.59% in year 2015, respectively, and the number of individuals with treatment will increase by 92.07%, 90.06%, and 63.00% in year 2015, as shown in Figure 8. It is interesting to note that relatively large testing rate leads to a slow increase of the number of individuals with treatment (shown in Figure 8(d)). Meanwhile, low recruitment rate U and large exit rate μ will also slow down the spread of HIV however their effects are very limited (as shown in Figure 7). This is because the number of HIV infected individuals is very small compared with the total high-risk population among MSM; thus the variation of S/N is very limited. This conclusion is in agreement with that for general high-risk population obtained by Xiao et al. [15].

Figure 6

Plots of estimated number of HIV/AIDS cases vary with transmission coefficient β 0 and diagnose rate δ. (a) Total number of HIV/AIDS cases. (b) Total HIV-positive cases. (c) Total AIDS cases. (d) Total HIV/AIDS with treatment. Other parameters used are shown in Table 2.

Figure 8

Plots of estimated number of HIV/AIDS cases vary with diagnose rate δ. (a) Total number of HIV/AIDS cases. (b) Total HIV-positive cases. (c) Total AIDS cases. (d) Total HIV/AIDS with treatment. Other parameters used are shown in Table 2.

Figure 7

Plots of estimated number of HIV/AIDS cases vary with constant recruitment U and exit rate μ. (a) Total number of HIV/AIDS cases. (b) Total HIV-positive cases. (c) Total AIDS cases. (d) Total HIV/AIDS with treatment. Other parameters used are shown in Table 2.

In order to investigate the effect of antiviral therapy, we study the variation in the incidence of HIV with different antiviral therapy efficacy represented by infectiousness (ϵ , ϵ , ϵ ), disease progression rates (η 1, η 2), and disease-related death rates (α and α ). In order to reduce the number of variables we introduce two factors v  (0 < v ≤ 1) and u  (0 < u ≤ 1) to describe variation in infectiousness and disease progression rates, respectively. During the analysis, v is set to be 0.75, 0.5 and 0.25 to represent that infectiousness is decreased to 75%, 50%, and 25%, respectively. We suppose that antiviral therapy decreased the progression rates to 5%. Figure 9 shows that HIV incidence will peak in around year 2016, and the better drug efficacy the lower HIV incidence. It indicates that although the drugs are persistently effective, HIV incidence will increase in the recent years. That is because effect of antiviral therapy on HIV epidemic between hosts is not instantaneous but delays some time.

Figure 9

Plots of incidence against factor v when antiviral therapy started immediately after diagnose. τ 1 = τ 2 = 0.8. Other parameters used are shown in Table 2.

5. Sensitivity Analysis

Note that in our model some parameters are known with uncertainties or have large variances, which may greatly affect outcomes. It is then necessary to do the uncertainty and sensitivity analysis such that the sensitive parameters can be detected. To examine the sensitivity of our results to parameter variations, we use latin hypercube sampling (LHS) and partial rank correlation coefficients (PRCCs) [38] to examine the dependence of the reproduction number R 0 and the expected number of total HIV/AIDS individuals in 2015 on each parameter.

We initially examine the sensitivity of basic reproduction number R 0 to parameter variations. We choose the sample size n = 4000, parameters interested as the input variables, and the value of R 0 as the output variable. Parameter values and ranges are listed in Table 1. Figure 10(a) shows the PRCC value of each parameter. Parameters with star above the bar are the significant ones, and the significance level we choose here is α = 0.05. It shows that β 0, ε , ε , ε , ε , and ε have positive impact upon R 0, whilst δ, τ 1, τ 2, α , α , and α have negative impact. From the PRCC results we know that R 0 is not sensitive to parameters ρ 1, ρ 2, η 1, and η 2. So, assumptions of ρ 1 = 2ρ/3, ρ 2 = ρ/3, and η 1, η 2 being half of ξ 1 and ξ 2 are justified and hence have little influence on our main results. Further, we can get that the most influential parameters are β 0, α since the PRCC values of them are larger than 0.8. In particular, simple calculation indicates that the basic reproduction number can decrease to less than 1 if β 0 reduces to 0.2.

Table 1

Parameter values and ranges.

Parameters	Ranges	Initial values	Parameters	Ranges	Initial values
U	[10000, 100000]	402450	S(0)	[70000, 150000]	763240
β 0	[0.1, 1]	0.8578	δ	[0.01, 1]	0.0816
ρ 1	[0.01, 1]	0.5879	ρ 2	[0.01, 1]	0.2939
η 1	[0.02, 1]	1/12	η 2	[0.02, 1]	1/6
τ 1	[0.02, 1]	0.2	τ 2	[0.02, 1]	0.2
ϵ I1	[0.01, 1]	0.7224	ϵ I2	[0.01, 1]	0.4177
ϵ A	[0.01, 1]	0.2495	α I	[0.02, 1]	0.172
ϵ T1	[0.01, 1]	0.4177	α T	[0.02, 1]	0.06
ϵ T2	[0.01, 1]	0.3507	α A	[0.02, 1]	0.1364

Figure 10

(a) Partial rank correlation coefficients (PRCC) results for the dependence of R 0 on each parameter. (b) Partial rank correlation coefficients (PRCC) results for the dependence of total HIV/AIDS cases in year 2015 on each parameter. ∗ denotes the value of PRCC is not zero significantly, where the significance level is 0.05.

We also investigate the sensitivity of the expected number of people living with HIV/AIDS to parameter variations. Given the current treatment policy, parameters described in Table 1 except for ϵ , τ 1, and η 1 are chosen as the input variables. The total HIV/AIDS cases in year 2015 is output variable. The PRCC values are shown in Figure 10(b), from which we can find that both U and S(0) have little impacts on the number of people living with HIV/AIDS in year 2015. The influences of other parameters on total HIV/AIDS cases are similar to their influences on R 0.

6. Discussion

According to the CD4+ T cell counts in the blood, we divided the HIV-positive individuals to several stages and formulated a mathematical model with antiviral therapy. The unknown parameters involved in this model were estimated using the MCMC simulation basing on the real data (i.e., the number of annual reported HIV/AIDS among MSM). We defined the threshold value (the basic reproduction number R 0) which determines whether the epidemic goes to extinction or not. By using the stability theories and methods of ordinary differential equation, we proved that the disease-free equilibrium is globally asymptotically stable when R 0 < 1, whilst the system is uniformly persistence when R 0 > 1.

We studied the effect of antiviral therapy in two situations: antiviral therapy started immediately once people are diagnosed with HIV and antiviral therapy started when CD4+ T counts are less than 350 cells per μL (current policy). We found that, given unchanged behaviors after treatment, there exists a critical value for the infectiousness below which immediate treatment is better than the current policy in terms of the reproduction number, whereas current policy exhibits better results than immediate treatment if the infectiousness is greater than the critical level. It indicates that if the treatment efficacy is really good (i.e., relatively low infectiousness), our conclusion suggests immediate treatment; otherwise the current policy is recommended. Similarly, when the infectiousness is relatively low (relatively good treatment efficacy), increasing treatment coverage will decrease the reproduction number and lead to decline in new HIV infection, whilst increasing treatment coverage will result in an increase in new HIV infection for the relatively great infectiousness, which is in agreement with that for heterosexual transmission [14]. Summing up the above, if treatment efficacy is relatively good, our conclusions suggest immediate treatment with high uptake rate; otherwise the current policy is reasonable.

Using the data on the number of new reported HIV/AIDS infected individuals by year among MSM, we obtained estimates of the reproduction number, intervention parameter values, and the high-risk population size. Our estimated reproduction number is 3.88 (95% CI 3.69–4.07) which is in the ranges of estimates for Western Germany (3.43–4.08) and UK (3.38–3.96) [39]. From the estimated parameters we know that the transmission coefficient β 0 is much larger than the estimation for heterosexual transmission [14] and general high-risk population [28]. This result is associated with the conclusion by Lou et al. [13] that MSM are 19 times more likely to be infected with HIV than general population. In fact, there are a lot of money boys in China, also called male sex workers especially in bathhouse, bars, and clubs [34, 40]. Recent surveys showed that the condom use rate among MSM is very low, only less than 30 percent [32]. Anotherstudy showed that anal intercourse between men, if unprotected, carries a high-risk of HIV transmission [3]. All these factors lead to the high probability for HIV transmission among MSM, which is in agreement with our estimation. Our estimation also shows that the diagnose rate among MSM is much lower than that for other high-risk population [15]. Meanwhile, we estimated that the antiviral therapy coverage rate among MSM in 2011 is less than the estimation by Ministry of Health, People's Republic of China [2]. This agrees with the conclusions obtained by Tong [34] that the antiviral therapy coverage rate for individuals infected sexually is lower than other infected groups.

Simulation results show that strengthening education to high-risk population and increasing surveillance and testing can slow down the spread of disease. Further, sensitivity analysis implies that the most influential parameters are infection rate β 0 and disease related death rate for HIV-positive individuals α . Note that high efficacy drug can reduce the transmission probability of HIV per high-risk behavior [18], and the education may reduce the contact rate and increase the condom use rate. This means that a high effective drug and timely education may effectively control HIV epidemic.

In this paper, we concluded that if the infectiousness for HIV/AIDS infected cases is relatively small, treatment started immediately once diagnosed is more beneficial to disease control. It should be mentioned that we have not considered costs of antiviral therapy. However, early antiviral therapy will increase the financial burden of the government of China and may increase high risk of occurrence of drug resistance. We will consider these factors in the future study.