# Modelling HIV/AIDS Epidemic among Men Who Have Sex with Men in China

Authors: Xiaodan Sun; Yanni Xiao; Zhihang Peng; Ning Wang
Journal: BioMed Research International
 (2013)

---

## Abstract

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.


---

## Body

## 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 Section2, 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 (I1), individuals with 200<l<350 (I2), 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 I1 and I2 will enter into T1 and T2 if they receive antiviral therapy, respectively, and finally progress to the AIDS stage (TA). 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 TA. The high-risk population size is represented by N; that is, N=S+I+I1+I2+A+T1+T2+TA,
(1)S′=U-λ(t)-(d+μ)S,I′=λ(t)-(d+αI+δ)I,I1′=ρ1δI-(d+αI+τ1+ξ1)I1,I2′=ρ2δI+ξ1I1-(d+αI+τ2+ξ2)I2,A′=(1-ρ)δI+ξ2I2-(d+αA)A,T1′=τ1I1-(d+αT+η1)T1,T2′=τ2I2+η1T1-(d+αT+η2)T2,TA′=η2T2-(d+αA)TA.Figure 1
The flow diagram of model with antiviral therapy.We assume that people enter into the susceptible class at a rateU, 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 εI1, εI2, εT1, εT2, and εA, respectively. Thus, susceptible people become HIV infected at a rate S(I+εI1I1+εI2I2+εT1T1+εT2T2+εA(A+TA))/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 I1 to I2 and from I2 to A, respectively. η1 and η2 denote the progression rates from T1 to T2 and from T2 to TA, respectively. αi(i=I,T,A) and τi(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+εI1I1+εI2I2+εT1T1+εT2T2+εA(A+TA))/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,I1,I2,A,T1,T2,TA)∈R+8:S+I+I1+I2+A+T1+T2+TA≤U/d} is positively invariant.


## 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 bed=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/40year-1. Zhou et al. [10] estimated the progression rate from I1 to I2 is ξ1=1/6 and the progression rate from I2 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 I2 and two-thirds are in stage I1 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 αI=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 αA=0.393*(1-65.3%)=0.136. The disease-related death rate for HIV-positive individuals with antiviral therapy is αT=0.172*(1-65.3%)=0.06.Recently, antiviral therapy is not provided to HIV-positive individuals withl>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: I1(0)=99, I2(0)=49, and A(0)=53. Considering that there are no HIV-positive individuals receiving antiviral therapy in year 2005, thus we have T1(0)=T2(0)=TA(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 εI1 as follows: the number of new infections in stage I1 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 I1 to stage I2 and from stage I2 to stage I3. By similar method we get that the numbers of new infections in I and I2 are 0.0407 and 0.017, respectively. Then, we have 1:εI1:εI2=0.0407:0.0294:0.017, which implies that εI1=0.7224, εI2=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 εT2, εA, and initial value of undiagnosed infected individuals I(0). Here we assume εT1=εT2 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 R0 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.


## 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 bed=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/40year-1. Zhou et al. [10] estimated the progression rate from I1 to I2 is ξ1=1/6 and the progression rate from I2 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 I2 and two-thirds are in stage I1 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 αI=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 αA=0.393*(1-65.3%)=0.136. The disease-related death rate for HIV-positive individuals with antiviral therapy is αT=0.172*(1-65.3%)=0.06.Recently, antiviral therapy is not provided to HIV-positive individuals withl>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: I1(0)=99, I2(0)=49, and A(0)=53. Considering that there are no HIV-positive individuals receiving antiviral therapy in year 2005, thus we have T1(0)=T2(0)=TA(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 εI1 as follows: the number of new infections in stage I1 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 I1 to stage I2 and from stage I2 to stage I3. By similar method we get that the numbers of new infections in I and I2 are 0.0407 and 0.017, respectively. Then, we have 1:εI1:εI2=0.0407:0.0294:0.017, which implies that εI1=0.7224, εI2=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 εT2, εA, and initial value of undiagnosed infected individuals I(0). Here we assume εT1=εT2 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 R0 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 numberR0, 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 R0 for system (1) (see Appendix B for details). For convenience, we let β0=βc(1-π), ω1=d+αI+δ, ω2=d+αI+τ1+ξ1, ω3=d+αI+τ2+ξ2, ω4=d+αA, ω5=d+αT+η1 and ω6=d+αT+η2. The expression of R0 is given as follows:
(2)R0=β0ω1[1+ρ1δεI1ω2+ρ1δξ1εI2ω2ω3+ρ2δεI2ω3+ρ1δξ1ξ2εAω2ω3ω42222+ρ2δξ2εAω3ω4+(1-ρ)δεAω4+ρ1δτ1εT1ω2ω52222+ρ1δτ1η1εT2ω2ω5ω6+ρ1δτ2ξ1εT2ω2ω3ω6+τ2ρ2δεT2ω3ω62222+ρ1δτ1η1η2εAω2ω5ω6ω4+ρ1δξ1τ2η2εAω2ω3ω6ω4+ρ2δτ2η2εAω3ω6ω4ρ1δεI1ω2].As the model shows, the primary infected individual will stay at the compartmentsI, I1, I2, A, T1, T2, and TA (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εI1, β0εI2, β0εA, β0εT1, β0εT2, and β0εTA. It is easy to see that β0/ω1, β0εI1/ω2, β0εI2/ω3, β0εA/ω4, β0εT1/ω5, β0εT2/ω6, and β0εTA/ω4 represent the numbers of individuals the primary case infects at each stage, respectively. Using the total probability formula we get the expression of R0.Meanwhile, system (1) has a disease-free equilibrium (DFE) E0=(S0,0,0,0,0,0,0,0), where S0=U/(d+μ). From [29] we can get that when R0<1, the disease-free equilibrium of system (1) is locally stable but unstable when R0>1. We can also prove that E0 is globally attractive when R0<1 (see details in Appendix C.1). Thus, we get the following theorem.Theorem 1.
WhenR0<1, the disease-free equilibrium of system (1) is globally asymptotically stable, but unstable when R0>1.Besides the disease-free equilibrium the system has an endemic equilibriumE+ when R0>1. It is not difficult to get the expression of E+=(S*,I*,I1*,I2*,A*,T1*,T2*,TA*), where S*=(U-ω1I*)/(d+μ), I1*=e1I*, I2*=e2I*, A*=e3I*, T1*=e4I*, T2*=e5I*, TA*=e6I*,  I*=U(R0-1)/(ω1R0+(d+μ)(1+e1+e2+e3+e4+e5+e6-(ω1/(d+μ)))), e1=ρ1δ/ω2, e2=(ρ2δ/ω3)+(ρ2δξ1/ω2ω3), e3=((1-ρ)δ/ω4)+(ρ2δξ2/ω3ω4)+(ρ1δξ1ξ2/ω2ω3ω4), e4=ρ1δ(1+σ)τ1/ω2ω5, e5=(ρ2δ(1+σ)τ2/ω3ω6)+(ρ1δξ1(1+σ)τ2/ω2ω3ω6)+(ρ1δ(1+σ)τ1η1/ω2ω5ω6), and e6=(ρ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 R0>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 R0>1, in the sense that there is a positive number δ0 such that for all initial values (S0,I0,I10,I20,A0,T10,T20,TA0)∈R+×Int(R+7), the solution of system (1) satisfies
(3)liminft→∞(S(t),I(t),I1(t),I2(t),A(t),T1(t),T2(t),TA(t))≥(δ0,δ0,δ0,δ0,δ0,δ0,δ0,δ0).


### 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 Table2. 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 R0, 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].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.
(a)
(b)Figure 3
MCMC plots forR0.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].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 Table2. Antiviral therapy started when CD4+ counts are less than 350 cells per μL.
(a)
(b)
(c)
(d)


### 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+αI+ξ1; we get the expression of basic reproduction number in this situation, which is denoted by R0*; then
(4)R0*=β0ω1[1+ρ1δεI1ω2*+ρ1δξ1εI2ω2*ω3+ρ2δεI2ω3+ρ1δξ1ξ2εAω2*ω3ω42222+ρ2δξ2εAω3ω4+(1-ρ)δεAω4+ρ1δτ2ξ1εT2ω2*ω3ω62222+τ2ρ2δεT2ω3ω6+ρ1δξ1τ2η2εAω2*ω3ω6ω4+ρ2δτ2η2εAω3ω6ω4].Situation 2. Start antiviral therapy immediately once people are diagnosed with HIV. In this situation, the expression of basic reproduction number R0 is given in Section 4.1.Using the parameter values listed in Table2 we get that R0=3.93, R0*=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 ofR0 is greater than R0* 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 R0 and R0* depends on many factors such as treatment coverage, infectiousness, and so forth. In the following we investigate the relation between R0 and R0*:
(5)R0-R0*=β0ρ1δω1[(εI1+ξ1εI2ω3+ξ1ξ2εAω3ω4+τ2ξ1εT2ω3ω6+ξ1τ2η2εAω3ω6ω4)22222222222×(1ω2-1ω2*)+τ1εT1ω2ω5+τ1η1εT2ω2ω5ω622222222222+τ1η1η2εAω2ω5ω6ω4].
In order to determine the sign of R0-R0*, we turn to consider
(6)ω1ω2(R0-R0*)β0ρ1δτ1=-1ω2*(εI1+ξ1εI2ω3+ξ1ξ2εAω3ω4+τ2ξ1εT2ω3ω6+ξ1τ2η2εAω3ω6ω4)+εT1ω5+η1εT2ω5ω6+η1η2εAω5ω6ω4=-1ω2*(εI1+ξ1εI2ω3)+1ω5εT1+(η1ω5-τ2ξ1ω2*ω3)(εT2+η2ω4ω6εA).
It is obvious that sign(R0-R0*)=sign(ω1ω2(R0-R0*)/β0ρ1δτ1), where sign(a) denotes the sign function and indicates the sign of a. So if εT1 and η1 satisfy the following:
(7)εT1ω5+(η1ω5-τ2ξ1ω2*ω3)(εT2+η2εAω4ω6)<1ω2*(εI1+ξ1εI2ω3),
we have R0-R0*<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 R0 and R0* 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 εT1, εT2, and εA; we then analyze the value of R0* and R0 as a function of v. Thus, v=0 denotes that the infected individuals are not infectious anymore and v=1 means the modification factors εT1, εT2, and εA 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 R0* and R0 in Figure 5(a). It shows that there is a critical value v*=0.81 such that R0<R0* for v<v* and R0>R0* 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.(a) Plots ofR0* and R0 against factor v. (b) Contour plot of R0 varies the factor v and treatment coverage rate, where τ1=τ2=τ. Parameters used are shown in Table 2.
(a)
(b)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 R0 with treatment coverage rate. We simply suppose τ1=τ2=τ. Solving the partial derivative of R0 with respect to τ we have
(8)∂R0∂τ=β0ω1[-ρ1δω22εI1-(ρ2δω32+ρ1δξ1(ω2+ω3)ω22ω32)εI2222+AεT1+BεT2+CεA(ρ2δω32+ρ1δξ1(ω2+ω3)ω22ω32)],
where
(9)A=ρ1δ(d+αI+ξ1)ω22ω5,(10)B=ρ1δη1(d+αI+ξ1)ω22ω5ω6+ρ2δ(d+αI+ξ2)ω32ω6+ρ1δξ1ω2ω3ω6-ρ1δξ1τ(ω2+ω3)ω22ω32ω6,(11)C=-ρ2δξ2ω32ω4-ρ1δξ1ξ2(ω2+ω3)ω22ω32ω4+ρ1δη1η2(d+αI+ξ1)ω22ω4ω5ω6+ρ2δη2(d+αI+ξ2)ω32ω4ω6+ρ1δξ1η2ω2ω3ω4ω6-ρ1δξ1η2τ(ω2+ω3)ω22ω32ω4ω6.We can verify thatA>0,B>0 (see Appendix D for details about the proof of B>0). Let ∂R0/∂τ=0; we have
(12)AεT1+BεT2=ρ1δω22+[ρ2δ+ρ1δξ1(ω2+ω3)ω22]εI2ω22-CεA≜E-.Thus, there exists a critical valueE- for a combination of parameters η1, η2, αT, and so forth such that ∂R0*/∂τ2<0 for AεT1+BεT2<E¯ and ∂R0*/∂τ2>0 for AεT1+BεT2>E¯. 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 R0 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 rateU, 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. Figure6 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].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.
(a)
(b)
(c)
(d)Plots of estimated number of HIV/AIDS cases vary with constant recruitmentU 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.
(a)
(b)
(c)
(d)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.
(a)
(b)
(c)
(d)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 (ϵT1, ϵT2, ϵA), disease progression rates (η1, η2), and disease-related death rates (αT and αA). 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 factorv when antiviral therapy started immediately after diagnose. τ1=τ2=0.8. Other parameters used are shown in Table 2.


## 4.1. Threshold Dynamics

The basic reproduction numberR0, 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 R0 for system (1) (see Appendix B for details). For convenience, we let β0=βc(1-π), ω1=d+αI+δ, ω2=d+αI+τ1+ξ1, ω3=d+αI+τ2+ξ2, ω4=d+αA, ω5=d+αT+η1 and ω6=d+αT+η2. The expression of R0 is given as follows:
(2)R0=β0ω1[1+ρ1δεI1ω2+ρ1δξ1εI2ω2ω3+ρ2δεI2ω3+ρ1δξ1ξ2εAω2ω3ω42222+ρ2δξ2εAω3ω4+(1-ρ)δεAω4+ρ1δτ1εT1ω2ω52222+ρ1δτ1η1εT2ω2ω5ω6+ρ1δτ2ξ1εT2ω2ω3ω6+τ2ρ2δεT2ω3ω62222+ρ1δτ1η1η2εAω2ω5ω6ω4+ρ1δξ1τ2η2εAω2ω3ω6ω4+ρ2δτ2η2εAω3ω6ω4ρ1δεI1ω2].As the model shows, the primary infected individual will stay at the compartmentsI, I1, I2, A, T1, T2, and TA (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εI1, β0εI2, β0εA, β0εT1, β0εT2, and β0εTA. It is easy to see that β0/ω1, β0εI1/ω2, β0εI2/ω3, β0εA/ω4, β0εT1/ω5, β0εT2/ω6, and β0εTA/ω4 represent the numbers of individuals the primary case infects at each stage, respectively. Using the total probability formula we get the expression of R0.Meanwhile, system (1) has a disease-free equilibrium (DFE) E0=(S0,0,0,0,0,0,0,0), where S0=U/(d+μ). From [29] we can get that when R0<1, the disease-free equilibrium of system (1) is locally stable but unstable when R0>1. We can also prove that E0 is globally attractive when R0<1 (see details in Appendix C.1). Thus, we get the following theorem.Theorem 1.
WhenR0<1, the disease-free equilibrium of system (1) is globally asymptotically stable, but unstable when R0>1.Besides the disease-free equilibrium the system has an endemic equilibriumE+ when R0>1. It is not difficult to get the expression of E+=(S*,I*,I1*,I2*,A*,T1*,T2*,TA*), where S*=(U-ω1I*)/(d+μ), I1*=e1I*, I2*=e2I*, A*=e3I*, T1*=e4I*, T2*=e5I*, TA*=e6I*,  I*=U(R0-1)/(ω1R0+(d+μ)(1+e1+e2+e3+e4+e5+e6-(ω1/(d+μ)))), e1=ρ1δ/ω2, e2=(ρ2δ/ω3)+(ρ2δξ1/ω2ω3), e3=((1-ρ)δ/ω4)+(ρ2δξ2/ω3ω4)+(ρ1δξ1ξ2/ω2ω3ω4), e4=ρ1δ(1+σ)τ1/ω2ω5, e5=(ρ2δ(1+σ)τ2/ω3ω6)+(ρ1δξ1(1+σ)τ2/ω2ω3ω6)+(ρ1δ(1+σ)τ1η1/ω2ω5ω6), and e6=(ρ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 R0>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 R0>1, in the sense that there is a positive number δ0 such that for all initial values (S0,I0,I10,I20,A0,T10,T20,TA0)∈R+×Int(R+7), the solution of system (1) satisfies
(3)liminft→∞(S(t),I(t),I1(t),I2(t),A(t),T1(t),T2(t),TA(t))≥(δ0,δ0,δ0,δ0,δ0,δ0,δ0,δ0).


## 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 Table2. 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 R0, 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].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.
(a)
(b)Figure 3
MCMC plots forR0.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].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 Table2. Antiviral therapy started when CD4+ counts are less than 350 cells per μL.
(a)
(b)
(c)
(d)


## 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+αI+ξ1; we get the expression of basic reproduction number in this situation, which is denoted by R0*; then
(4)R0*=β0ω1[1+ρ1δεI1ω2*+ρ1δξ1εI2ω2*ω3+ρ2δεI2ω3+ρ1δξ1ξ2εAω2*ω3ω42222+ρ2δξ2εAω3ω4+(1-ρ)δεAω4+ρ1δτ2ξ1εT2ω2*ω3ω62222+τ2ρ2δεT2ω3ω6+ρ1δξ1τ2η2εAω2*ω3ω6ω4+ρ2δτ2η2εAω3ω6ω4].Situation 2. Start antiviral therapy immediately once people are diagnosed with HIV. In this situation, the expression of basic reproduction number R0 is given in Section 4.1.Using the parameter values listed in Table2 we get that R0=3.93, R0*=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 ofR0 is greater than R0* 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 R0 and R0* depends on many factors such as treatment coverage, infectiousness, and so forth. In the following we investigate the relation between R0 and R0*:
(5)R0-R0*=β0ρ1δω1[(εI1+ξ1εI2ω3+ξ1ξ2εAω3ω4+τ2ξ1εT2ω3ω6+ξ1τ2η2εAω3ω6ω4)22222222222×(1ω2-1ω2*)+τ1εT1ω2ω5+τ1η1εT2ω2ω5ω622222222222+τ1η1η2εAω2ω5ω6ω4].
In order to determine the sign of R0-R0*, we turn to consider
(6)ω1ω2(R0-R0*)β0ρ1δτ1=-1ω2*(εI1+ξ1εI2ω3+ξ1ξ2εAω3ω4+τ2ξ1εT2ω3ω6+ξ1τ2η2εAω3ω6ω4)+εT1ω5+η1εT2ω5ω6+η1η2εAω5ω6ω4=-1ω2*(εI1+ξ1εI2ω3)+1ω5εT1+(η1ω5-τ2ξ1ω2*ω3)(εT2+η2ω4ω6εA).
It is obvious that sign(R0-R0*)=sign(ω1ω2(R0-R0*)/β0ρ1δτ1), where sign(a) denotes the sign function and indicates the sign of a. So if εT1 and η1 satisfy the following:
(7)εT1ω5+(η1ω5-τ2ξ1ω2*ω3)(εT2+η2εAω4ω6)<1ω2*(εI1+ξ1εI2ω3),
we have R0-R0*<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 R0 and R0* 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 εT1, εT2, and εA; we then analyze the value of R0* and R0 as a function of v. Thus, v=0 denotes that the infected individuals are not infectious anymore and v=1 means the modification factors εT1, εT2, and εA 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 R0* and R0 in Figure 5(a). It shows that there is a critical value v*=0.81 such that R0<R0* for v<v* and R0>R0* 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.(a) Plots ofR0* and R0 against factor v. (b) Contour plot of R0 varies the factor v and treatment coverage rate, where τ1=τ2=τ. Parameters used are shown in Table 2.
(a)
(b)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 R0 with treatment coverage rate. We simply suppose τ1=τ2=τ. Solving the partial derivative of R0 with respect to τ we have
(8)∂R0∂τ=β0ω1[-ρ1δω22εI1-(ρ2δω32+ρ1δξ1(ω2+ω3)ω22ω32)εI2222+AεT1+BεT2+CεA(ρ2δω32+ρ1δξ1(ω2+ω3)ω22ω32)],
where
(9)A=ρ1δ(d+αI+ξ1)ω22ω5,(10)B=ρ1δη1(d+αI+ξ1)ω22ω5ω6+ρ2δ(d+αI+ξ2)ω32ω6+ρ1δξ1ω2ω3ω6-ρ1δξ1τ(ω2+ω3)ω22ω32ω6,(11)C=-ρ2δξ2ω32ω4-ρ1δξ1ξ2(ω2+ω3)ω22ω32ω4+ρ1δη1η2(d+αI+ξ1)ω22ω4ω5ω6+ρ2δη2(d+αI+ξ2)ω32ω4ω6+ρ1δξ1η2ω2ω3ω4ω6-ρ1δξ1η2τ(ω2+ω3)ω22ω32ω4ω6.We can verify thatA>0,B>0 (see Appendix D for details about the proof of B>0). Let ∂R0/∂τ=0; we have
(12)AεT1+BεT2=ρ1δω22+[ρ2δ+ρ1δξ1(ω2+ω3)ω22]εI2ω22-CεA≜E-.Thus, there exists a critical valueE- for a combination of parameters η1, η2, αT, and so forth such that ∂R0*/∂τ2<0 for AεT1+BεT2<E¯ and ∂R0*/∂τ2>0 for AεT1+BεT2>E¯. 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 R0 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 rateU, 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. Figure6 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].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.
(a)
(b)
(c)
(d)Plots of estimated number of HIV/AIDS cases vary with constant recruitmentU 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.
(a)
(b)
(c)
(d)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.
(a)
(b)
(c)
(d)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 (ϵT1, ϵT2, ϵA), disease progression rates (η1, η2), and disease-related death rates (αT and αA). 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 factorv 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 R0 and the expected number of total HIV/AIDS individuals in 2015 on each parameter.We initially examine the sensitivity of basic reproduction numberR0 to parameter variations. We choose the sample size n=4000, parameters interested as the input variables, and the value of R0 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, εI1, εI2, εA, εT1, and εT2 have positive impact upon R0, whilst δ, τ1, τ2, αI, αT, and αA have negative impact. From the PRCC results we know that R0 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, αI 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
ϵ
I
1
[0.01, 1]
0.7224
ϵ
I
2
[0.01, 1]
0.4177
ϵ
A
[0.01, 1]
0.2495
α
I
[0.02, 1]
0.172
ϵ
T
1
[0.01, 1]
0.4177
α
T
[0.02, 1]
0.06
ϵ
T
2
[0.01, 1]
0.3507
α
A
[0.02, 1]
0.1364Table 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
ϵ
I
1
Modification factor for HIV infected individuals withl>350
0.7224
—
[5]
ϵ
I
2
Modification  factor for HIV infected individuals with200<l<350
0.4177
—
[5]
ϵ
T
1
Modification factor for HIV infected individuals withl>350 and receiving antiviral therapy
0.3507
—
—
ϵ
T
2
Modification factor for HIV infected individuals with200<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 withl>350
0.5879
—
[10]
ρ
2
Proportion of diagnosed HIV-positive individuals with200<l<350
0.2939
—
[10]
ξ
1
Progression rate fromI1 to I2
1/6
—
[10]
ξ
2
Progression rate fromI2 to A
1/3
—
[10]
η
1
Progression rate fromT1 to T2
1/12
—
[12]
η
2
Progression rate fromT2 to TA
1/6
—
[12]
τ
1
Antiviral therapy coverage rate for HIV-positive individuals withl>350
0
—
—
τ
2
Antiviral therapy coverage rate for HIV-positive individuals with200<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 ofS (high-risk susceptible)
763240
87285
MCMC
I
(
0
)
Initial value ofI (individuals living with undiagnosed HIV)
5988
1274
MCMC
I
1
(
0
)
Initial value ofI1 (HIV-positive individuals with l>350)
99
—
Database
I
2
(
0
)
Initial value ofI2 (HIV-positive individuals with 200<l<350)
49
—
Database
A
(
0
)
Initial value ofA (AIDS patients diagnosed at onset of AIDS symptoms)
53
—
Database
T
1
(
0
)
Initial value ofT1 (HIV-positive individuals with l>350 and receiving therapy)
0
—
Database
I
2
(
0
)
Initial value ofT2 (HIV-positive individuals with 200<l<350 and receiving therapy)
0
—
Database
T
A
(
0
)
Initial value ofTA (AIDS patients diagnosed before onset of AIDS symptoms)
0
—
Database
R
0
The basic reproduction number
3.8840
0.097
Calculated(a) Partial rank correlation coefficients (PRCC) results for the dependence ofR0 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.
(a)
(b)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 Table1 except for ϵT1, τ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 R0.


## 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 numberR0) 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 R0<1, whilst the system is uniformly persistence when R0>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 αI. 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.



---

*Source: 413260-2013-09-30.xml*