Trace the History of HIV and Predict Its Future through Genetic Sequences.

Traditional methods of quantifying epidemic spread are based on surveillance data. The most widely used surveillance data are normally incidence data from case reports and hospital records, which are normally susceptible to human error, and sometimes, they even can be seriously error-prone and incomplete when collected during a destructive epidemic. In this manuscript, we introduce a new method to study the spread of infectious disease. We gave an example of how to use this method to predict the virus spreading using the HIV gene sequences data of China. First, we applied Bayesian inference to gene sequences of two main subtypes of the HIV virus to infer the effective reproduction number (GRe(t)) to trace the history of HIV transmission. Second, a dynamic model was established to forecast the spread of HIV medication resistance in the future and also obtain its effective reproduction number (MRe(t)). Through fitting the two effective reproduction numbers obtained from the two separate ways above, some crucial parameters for the dynamic model were obtained. Simply raising the treatment rate has no impact on lowering the infection rate, according to the dynamics model research, but would instead increase the rate of medication resistance. The negative relationship between the prevalence of HIV and the survivorship of infected individuals following treatment may be to blame for this. Reducing the MSM population's number of sexual partners is a more efficient strategy to reduce transmission per the sensitivity analysis.

1. Introduction

AIDS (Acquired Immunodeficiency Syndrome) is a disease caused by HIV (Human Immunodeficiency Virus) and has become a disease that seriously threatens human health in the world today. In order to alleviate the AIDS epidemic in China, the Chinese government began to implement the “Four Frees and One Care” HIV/AIDS free treatment policy in 2003 [1,2]. This policy has greatly improved the patient’s survival rate and the patient’s life quality. On the other hand, the high mutation rate of HIV virus has led to the generation of HIV-resistant strains [3]. The spread of drug-resistant strains causes a waste of national medical resources.

We want to know the effect of the free treatment policies, especially the drug resistance situation of HIV/AIDS by dynamic models. It is hard to avoid parameter acquisition when performing prediction with dynamics models. Some parameters are general parameters, which are common to both China and abroad. We can assign these values by searching relative literature. Some might be greatly different due to different countries or different treatment strategies, so they need be fitted by monitoring data. However, such areas happen to lack monitoring data except for some gene sequences of HIV patients being detected before treatment. In this manuscript, we will show how to combine the gene sequences of HIV and the dynamic model to study the dynamics of HIV spreading.

HIV is an RNA virus. RNA viruses have a short reproduction time and a high mutation rate, which can easily lead to a large number of genetic mutations during the transmission process. Therefore, during the spread of HIV, relevant signals will be left in the sampled virus sequence. In recent years, researchers have begun to conduct retrospective studies on the spread of the HIV virus on the genetic level. Most of these studies are based on genetic sequences to build molecular networks of transmission using Bayesian evolution analysis [4,5,6,7,8,9,10,11,12]. For example, in 2008, Lewis et al. used the Bayesian framework to conduct a phylogenetic analysis of 2126 HIV gene sequences in London, and they constructed a network for HIV transmission, revealing that the local HIV/AIDS epidemic can be traced back to the late 1990s [13]. In 2018, Zhao et al. constructed a molecular network of the spread of MSM in Beijing based on the HIV gene sequence, and they proposed the need to strengthen education and intervention for people with potential high risk. These findings have important implications for the parameterization of epidemiological models and the design of intervention strategies.

Traditional dynamic models can be used to predict the spread of virus in populations [14,15,16,17,18,19,20,21,22], and the Bayesian evolution analysis of virus separated from the infected people can trace the historical dynamics of virus transmission [13,23,24]. However, the work connecting the gene sequences and the dynamic model has not yet been seen. This manuscript will use Anhui province as an example, combine genetic sequence molecular evolution analysis and system dynamics models to conduct an in-depth research on the HIV/AIDS epidemic in the content of the HIV-free treatment policy in China, and thus predict the drug resistance spreading in the next few years.

2. Research Methodology

2.1. Gene Sequence Analysis

To trace the history of transmission dynamics of HIV, we sampled totally 360 HIV- gene sequences (CRF 01AE (150), CRF 07BC (190), CRF 08BC (9) and CRF 5501B (11)) from 360 HIV patients who were infected with HIV through sexual behavior in Anhui Province. This data set is sampled between October 2017 and September 2018. Considering that the proportions of the CRF 08BC and CRF 5501B are too small, we only analyzed the subtypes CRF 01AE and CRF 07BC, comprising a total of 340 HIV virus gene sequences. We applied Mega for sequence alignment. The total length of the gene that we studied was 1075 bp.

The selection of a nucleotide substitution model is an important step during the phylogenetic analysis. Common nucleotide substitution models include the HKY () model, GTR () model, and so on [25]. These different substitution models determined different evolutionary distances and different phylogenetic trees [26]. Here, we used software Phylosuite to select the most suitable substitution model according to the Akaike Information Criterion (AIC) [27], where was defined as: , in which k is the number of corresponding parameters and L is the likelihood value. Method-of-moments was used to estimate the method. The model with lower value is the more suitable model. After calculating the AIC score, the HKY + I model and the uncorrelated exponential relaxed molecular clock are determined finally.

2.1.1. Bayesian Inference

The birth death skyline (BDSKY) model in BEAST2 is a random description of population change, which allows the extraction of virus transmission dynamics information from a phylogenetic tree [11]. This model is established on Bayesian inferences, which can estimate the effective reproduction number of HIV/AIDS transmission from gene sequences. The Bayesian approach allows us to continuously update our estimate of a set of parameters, , as data become available.

, the prior distribution, represents our prior estimates about the true value of . is the likelihood distribution. It is also often written as , which means the probability of observing the data given . For the method to work, it is necessary to calculate the likelihood distribution for all possible values of . is the model evidence, and it is the same for all possible hypotheses (values of ) being considered. is the posterior distribution and represents our updated estimate of the value of given the observed data.

The main objective of Bayesian inference is to calculate the posterior distribution of our parameters using our prior beliefs updated with our likelihood. From the posterior distribution, we can determine the most likely values of given the observed data. Since we are usually only interested in relative probabilities of different hypotheses, can be left out of the calculation, and we write the model form of Bayes’ theorem as where ∝ means “proportional to”.

2.1.2. The Effective Reproduction Number Infered from Gene Sequences

Since we will study two different subscript gene sequences, we use for the subtype CRF 07BC and for the subtype CRF 01AE. We will deduce the propagation dynamics of the two subtypes separately. For each subtype j, the BDSKY process can conclude three parameters: transmission rate , removal rate and sampling probability . This process allows individuals to “birth” () or “death” ( and ) at any point in time. A “birth event” corresponds to an individual’s infection, and a “death event” corresponds to an infected individual becoming non-infected (i.e., removed from the infected compartment because of individual death, successful treatment, or individual behavior change). Then, the effective reproduction number derived from the gene sequences can be defined as [11] in which the left-hand subscript G denotes that the effective reproduction number is derived from the gene sequences. For estimating of subtype j or , the Bayes’ theorem that we use is

In order to catch the characteristics of transmission in different historical stages, we divided the period from the origin of each subtype HIV sequences to the sampling time of the last sample into a total of six segments. Then, from each subtype gene sequence of j and , we can infer the effective reproduction number during each segment t, through Equation (1). Hence, is a piece wise function for each of the two subtypes, which is a constant during the same time segment t. It refers to the average number of new infections caused by an infected person at a certain time point t during the outbreak. We use the Markov Chain Monte Carlo (MCMC) method to achieve Bayesian phylogeny inference.

2.2. Transmission Dynamics Model

2.2.1. A Nested Model and Its Four Submodels

To predict the prevalence of HIV drug-resistant strains, we considered the transmission dynamics of HIV/AIDS through sexual route, including gay, straight and bisexual. The whole population was divided into four subgroups: heterosexual women who only have sex with men (marked as w), heterosexual men who only have sex with women (marked with ), homosexual men who only have sex with men (marked with ), and bisexual men who have sex with both women and men (marked as b). The population was divided into susceptible people S, infected people I (infected but not receiving treatment), treated people T, and drug-resistant people R according to the infection. Considering that the HIV infection process has obvious stage characteristics, according to the level of CD4 cells in the patient’s body, we divide the infected people into three stages: the first infection stage (CD4 count cells/mm); the second infection stage (cells/mm< CD4 counts cells/mm); and the third stage of infection (AIDS stage, CD4 count cells/mm). Similarly, and , can be explained. The details of the model can be found in Figure 1. When the superscript g in the transmission process diagram Figure 1 is taken as and b, it represents the four types of subgroups. See Appendix A.1 for the corresponding dynamic model.

Figure 1

Transmission process diagram. The superscript g can be and b, respectively, which represents the four types of subgroup populations. Parameter a is an input of susceptible individuals; are the disease progression rates in infection stage 1, 2 and 3, respectively; are the treatment rates of patients in stage 1, 2 and 3, respectively; is the percentage of patients who give up HAART; is rate of drug resistance after HAART; is the percentage of patients who are no longer resistant.

The standard treatment of HIV/AIDS in China has gone through four periods: the period before 2002 featured no treatment; the period from 2003 to 2011 was when the treatment was started when the CD4 count was less than cells/mm; the period from 2012 to 2015 was when the treatment was started when the CD4 count was less than cells/mm and the period after 2015 was when the treatment was started soon after discovery. We call these periods I (1992–2002), II (2003–2011), III (2012–2015) and IV (2016–now), respectively.

Due to different prevention and control policies being held in different historical stages, the dynamics of HIV transmission have different stage characteristics. So, our model (Figure 1) in fact is a nested model. For the four different historical treatment periods I, II, III and IV, we can simplify the model to be four submodels by letting some parameters be zero based on its history characters; see Table 1 for the four different submodels. According to the “next-generation operator” method in the literature [28], we can obtain the basic reproduction number for each submodel of corresponding historical period I, II, III and IV (see Appendix A.2). Here, the left-hand subscript M means that the basic reproduction number was obtained from I, II, III or IV, one of the four models respectively.

Table 1

Submodels of different historical periods and corresponding assumptions.

Submodel	Assumptions	Historical Period
I	δ_i=0 and Tig=0,Rig=0, i=1,2,3, g=w,mw,mm,b	1992–2002
II	δ_i=0 and Tig=0, i=1,2, g=w,mw,mm,b	2003–2011
III	δ_1=0 and T1g=0, g=w,mw,mm,b	2012–2015
IV	The full model	2016–now

2.2.2. The Effective Reproduction Number Inferred from the Dynamic Submodels

For a general dynamic model of an infectious disease, the basic reproduction number is the expected number of infections directly generated by one case given that all individuals are equally susceptible. As the infection spreads, the susceptibility of the population decreases. The effective reproduction number, , is related to the basic reproduction number , by where is the average susceptibility of the population. is often used as an indicator of the effectiveness of interventions, such as social distancing measures, to contain the spread of a virus. If is greater than 1.0, the infection is growing at an exponential rate. If is at 1.0, the spread is sustained at a linear rate. If is less than 1.0, the infection is spreading at a slower pace and will eventually die out.

Specific to our study on HIV/AIDs spreading during four different historical periods I (1992–2002), II (2003–2011), III (2012–2015), and IV (2016–now), let us make where the basic reproduction number was obtained from Section 2.2.1 and and represent the susceptible population and the total target population of each historical period I, II, III and IV, respectively [29].

3. Results

3.1. Tracing Back the Dynamic History of HIV/AIDS through Bayesian Inference

The MCMC converged after 20 million iterations for both subtypes. We discarded the first 2 million iterations. The phylogenetic tree and model parameters were sampled every 1000 iterations. All parameter estimates showed the effective sample size (ESS) of more than 200. Some results are shown in Table 2.

Table 2

Epidemiological parameter estimations and their 95% HPD of subtype CRF 01AE and CRF 07BC.

Subtype		MRCA	The Effective Reproduction Number j=01AE,07BC
		Median	GjR_e(1)	GjR_e(2)	GjR_e(3)	GjR_e(4)	GjR_e(5)	GjR_e(6)
CFR01AE	95% HPD lower95% HPD upper	1982	1.64	1.47	1.42	1.37	1.23	2.17
		1973	1.08	0.97	0.94	0.91	0.83	1.61
		1990	2.30	2.08	2.00	1.94	1.72	2.87
CFR07BC	95% HPD lower95% HPD upper	1992	1.58	1.52	1.51	1.50	1.40	2.21
		1976	1.15	1.11	1.10	1.10	1.03	1.72
		1998	2.10	2.00	1.98	1.97	1.81	2.75

For each data set and substitution model, we analyze them with an uncorrelated exponential relaxed molecular clock and a strict clock model. In the study of virus evolution, we used nucleotide substitution models including the HKY () model and GTR () model [25]. The most appropriate combination was selected according to the Akaike Information Criterion (AIC) [27], where was defined as , in which k is the number of corresponding parameters and L is the likelihood value. Method-of-moments was used to estimate the method. Finally, the uncorrelated exponential relaxed molecular clock and HKY + I model is determined by calculating the value. The important epidemiological parameter estimation and their 95% HPD of subtype CRF 01AE and CRF 07BC are shown in Table 2.

The results in Table 2 show that the most recent common ancestor (MRCA) to the 150 CRF 01AE subtype gene sequences was around 1982, and the 95% highest posterior density (HPD) is between 1973 and 1990. Correspondingly, the MRCA to these 190 CRF 07BC subtype gene sequences was relatively late, around 1992 (95% HPD: 1976-1998). Most importantly, we obtained the effective reproduction number over time segment t from the two subtypes, respectively, in which represents the time segment. The effective reproduction number of the two subtypes showed some similar transmission characteristics in terms of time: both of them were in a relatively stable state in the early stage, kept around 1.6 with slight fluctuations, and showed a slight downward trend from 2008 to 2012. However, things changed a lot around 2012: for some reason, the effective reproduction number increased rapidly and then remained at about 2.2.

When we study the history of HIV transmission dynamics in the entire population through the two effective reproduction numbers , we will have to integrate them together. The two effective reproduction numbers originated at different times but ended at the same time. We took the common historical stage of them. So, the research period was from 1992 to 2018. During this common historical period, we denote the average of the effective reproduction number for the entire population as , which changes along with time t.

In order to figure out the formula for , we used the proportion of the number of sequences of each subtype in the total gene sequence as its weight in the effective reproduction number, since subtypes widely disseminated in the population should have more contribution for . Let be the sequence number of gene subtype j (, ). Then, the weighted average effective reproduction number of the entire population at time t can be described by Formula (3).

Use R software to visualize ; then, the result can be found in Figure 2. Figure 2 shows the median line of and its 95% HPD interval from 1992 to 2018. Before 2011, we can see that the epidemic showed a relatively stable transmission trend through sexual route. The value of remained stable and slightly decreased, with a median value of about 1.33 (95% HPD: 0.94–1.8). However, from 2011 to 2014, rose rapidly to 2.13 (95% HPD: 1.32–2.83) for some reason. After 2015, remained stable, with its median value maintained at approximately 2.20 (95% HPD: 1.67–2.85). Some explanations will be given later on the possibility of rising in this period.

Figure 2

The weighted average of the effective reproduction number and its 95% HPD interval.

3.2. Dynamic Model Analysis: Parameter Estimations and Predicting the Transmission of HIV/AIDS

3.2.1. Parameter Estimations

As we all know, an important step in the prediction of infectious disease spreading is to obtain (or estimate from data) parameters values. For most of the general parameters we treated them as constants from references. Details of the values and references can be found in Appendix A.4. However, there are still some epidemiological related parameters that we are not sure of and need to be estimated from data.

Unlike novel coronavirus patients, HIV-infected people do not have obvious infection symptoms. Many HIV infected people do not know that they have been infected for a long time. This prevents the government from obtaining information about people living with HIV as early as possible. Moreover, the Chinese government’s screening efforts for HIV-infected people are far less than those for infected novel coronavirus. As a result, recorded HIV infection data, if available, may not be very accurate and may not be suitable for parameter fitting. In Anhui province, which we studied, there is not much surveillance data on HIV/AIDS that can be used. One form of reliable data available are 340 genetic sequences of HIV-infected people taken before treatment.

In the previous subsections, we have extrapolated the dynamics of HIV transmission () from 1992 to 2018 using the 340 genetic sequences. In order to use these gene data to estimate unknown parameters of the model (Figure 1, which consists of four different submodels I (1992–2002), II (2003–2011), III (2012–2015) and IV (2016–now), depending on the historical period from 2003–now), we calculated the effective reproduction number for each submodel I, II, III and IV by Equation (2). Obviously, is functions of parameters and time t. Since the two effect reproduction numbers obtained through different methods should be the same at the same time t, i.e., , we fit them at the same historical period. In this way, we explored the humanistic and statistical features in different historical stages. We used MCMC to realize this. The algorithm ran for 2.00 million iterations, and we adapted the proposal distribution after 2 million iterations using Geweke’s method [30] to assess convergence. The fitted curves of the two effective reproduction numbers for four historical periods I, II, III and IV are shown in Figure A1 in Appendix A.4. The fitted parameters obtained from the above four periods are shown in Table 3. The value of these parameters will be referred to as baseline parameter values in the following research.

Figure A1

The fitted curves of the two effective reproduction numbers and for four historical periods I (a), II (b), III (c) and IV (d), respectively.

Table 3

Model parameters and the 95% HPD fitted by the two effective reproduction numbers.

Parameters	Description	Value	Fitted from Submodel
a^w	input of heterosexual women	343,780	IV
a^mw	input of heterosexual men	358,200	IV
a^mm	input of homosexual men	4974 (2936, 7012)	IV
a^b	input of bisexual men	4646 (2662, 6630)	IV
c^w	the average number of sexual partners per year for heterosexual women	1.9957 (0.846, 3.1454)	IV
c^mw	the average number of sexual partners per year for heterosexual men	2.0946 (0.8514, 3.1742)	IV
c^mm	the average number of sexual partners per year for homosexual men	5.4984 (2.0701, 7.8565)	IV
c^b	the average number of sexual partners per year for bisexual men	5.5765 (2.1816, 8.3334)	IV
ϵ	transmissibility ratio of drug-resistant to untreated wild strain carriers	0.43983	II
1/θ1T	course of the disease after treatment in stage 1	7.4030 (6.8823, 8.0000 )	IV
1/θ2T	course of the disease after treatment in stage 2	10.8948 (10.0000, 12.5000)	III
1/θ3T	course of the disease after treatment in stage 3	4.1969(3.3333, 4.7619 )	II
δ_1	treatment rate of stage 1	0.50188 (0.2155, 0.7883)	IV
δ_2	treatment rate of stage 2	0.75282 (0.6086, 0.8971)	IV
δ_3	treatment rate of stage 3	0.74938 (0.6057, 0.8930)	IV
v	the percentage of male partners of bisexual men	0.62567 (0.4549, 0.7965)	IV

3.2.2. HIV/AIDS Epidemic Trend under Baseline Parameters

We first studied the dynamic process of the total number of infected people and the proportion carrying drug-resistant strains over time from 2015 to 2025 in the entire population of Anhui province under baseline parameters (Figure 3). The results show that the total number of HIV/AIDS infections during this period appears to be an upward trend, and the total number of surviving infections will increase to 33,260 by 2025 (95% confidence interval (CI): 28,096–40,090, Figure 3a). The number of new infections is also increasing and will reach 4522 (95% CI: 3366–6130) by 2025. According to the Anhui Provincial Health and Family Planning Commission, from 2017 to 2019, the cumulative reports of HIV/AIDS surviving patients whose current address is Anhui province are 18619 [31], 17183 [32] and 19604 [33], respectively. They are indicated by red dots in Figure 3a. Considering that the detection rate of HIV infections in Anhui province has not reached 100%, the number of reported cases should be less than the actual number of HIV infections. Based on the reported HIV-positive people and the HIV testing rate [31] in the corresponding year, we have reason to say that our predicted results are basically consistent with the actual situation.

Figure 3

The number of HIV-positive people in the entire population (a) and drug resistance rate (b), where the continuous curve is the median line. Red dots in (a) are reported cases of 2017, 2018 and 2019 respectively.

The target population of our study is the general population. We studied the proportion of each population in the total new infections under baseline parameters from 2015 to 2025. The results showed that the proportion of heterosexual HIV-positive people decreased from 30.28% in 2015 to 20.43% in 2025 (95% CI: 15.92–26.09%, of which female infected people accounted for 15.61%, while heterosexual males accounted for 4.82%). In contrast, the proportion of HIV-positive MSM (homosexual and bisexual men) has steadily increased. By 2025, 78.84% (95% CI: 68.97–88.70%) of new HIV infections will be infected through homosexual sex behavior. Among them, homosexual HIV-positive people accounted for 53.87%, and bisexual HIV-positive people accounted for 24.97%.

Regarding the subpopulations, although the HIV infection rate among heterosexual men has declined and tends to be stable, the rates in other populations have shown steady upward trends. By 2025, the HIV-positive rate among heterosexual women will reach 0.25% (95% CI: 0.21–0.28%), while the HIV-positive rate among heterosexual men will be lower, only about half of the positive rate of women (0.12%, 95% CI: 0.10–0.13%). In contrast, men who have sex with men belong to high-risk groups, and their HIV positive rate is dozens of times that of heterosexual men: the HIV infection rate of bisexual men is 2.44% (95% CI: 1.88–3.12%), while that of homosexual men is much higher, reaching 3.62% (95% CI: 2.81–4.64%) by 2025.

Considering that large-scale treatment for HIV-infected people in China has lasted for quite a long time, we are concerned that drug-resistant strains (secondary drug resistance) will appear in the patients due to poor medication compliance, which will lead to the wide spread of drug-resistant strains (primary drug resistance) in the population. Unsurprisingly, the proportion of HIV-positive people carrying drug-resistant strains in Anhui province has shown a rapid upward trend (Figure 3b) and will basically reach the highest point of 9.88% by 2025, which corresponds to 3282 (95% CI: 2862–3812) in terms of quantity. The proportion of primary drug-resistant patients in the total newly infected population will also increase over time. By 2025, the number of primary drug-resistant patients will account for 9.63% of total patients carrying drug-resistant strains (95% CI: 6.56–14.2%). In addition, the proportion of primary drug-resistant patients among total new infections first showed an upward trend, then tended to be stable each year, and finally will increase to 7.01% (95% CI: 6.61–7.44%) by 2025.

3.2.3. Influence of Intervention Measures on HIV/AIDS Epidemic Trend

In order to implement the “Healthy China 2030” program and deepen the medical and health system reformation, China’s “13th Five-Year Plan for Combating and Prevention of AIDS” in 2017 set the “90–90–90” goal. That is, the proportion of infected persons detected by testing should be above 90%; the proportion of diagnosed infected persons receiving antiviral treatment should be above 90%; the treatment success rate of infected persons receiving antiviral treatment should be above 90% [34] by 2020. In the baseline parameters, the treatment rates of the three different infection stages were 50.19% (95% CI: 21.55%, 78.83%), 75.28% (95% CI: 60.86%, 89.71%) and 74.94% (95% CI: 60.57%, 89.30%). These treatment rates are still far from the second “90%”. Next, we consider what will happen if the treatment rate increases to the target state in 2020.

The total HIV infection rate and drug resistance of Anhui province from 2020 to 2025 under the “90%” treatment rate are shown in Figure 4, which is a comparison chart between the baseline and the target situation. The results show that the total infection rate still keeps increasing (the green curve in Figure 4a) even under the “90%” treatment rate, which is only slightly lower than the baseline parameters (the red curve in Figure 4a). Even more surprising, the total drug resistance rate after intensified treatment will increase more rapidly than the baseline situation, reaching 10.92% (95% CI: 10.69–11.12%) by 2025, which is seen as the red curve in Figure 4b. In addition, the proportion of primary drug resistance in new infections has also shown an increasing trend compared with baseline parameters.

Figure 4

Impact of target treatment rates on the HIV/AIDS epidemic: the infected rate (a) and the rate of drug resistance (b). The green box plots are based on the baseline parameters, while the red ones are based on target treatment rate. The continuous green and red curves are their median lines respectively.

In short, when we increase the treatment rate to 90%, the total infection rate can be reduced only by 0.017% by 2025, but in turn, the total drug resistance rate and the proportion of primary drug resistance in new infections will increase 1.06% and 1.12%, respectively. What prevents ART treatment from playing an effective role in the epidemic? Is the current drug treatment effect not good enough? If the new first-line drugs can improve the treatment effect (including reducing the infectivity and prolonging the lifespan), how much will the epidemic be improved?

To answer these questions, we assume the following two situations: (a) the new drug can only prolong the lifespan after treatment by 20% (it can survive for 24 years after extension); (b) the new drug can not only prolong the life after treatment by 20% but also reduce the infection rate by 20%. The results are shown in Figure 5. We were surprised to find that assumption (a) will make the epidemic worse instead (Figure 5a, green line). This may be because patients will have more risky sexual behaviors due to longer lifespan. However, assumption (b) will slow down the trend (Figure 5a, red line). In terms of drug resistance, the drug resistance rate will decrease slightly under assumption (a) (Figure 5b, green line) and increase under assumption (b) (Figure 5b, red line). The main reason may be that treatment only suppresses the non-resistant strains, which may instead give drug-resistant strains more space for transmission. This result suggests that prolonging the lifespan of the patients is a double-edged sword. When developing new drugs, we should not only consider how to prolong the life span of patients but also consider how to reduce the infectivity. This conclusion will continue to be discussed in the sensitivity analysis.

Figure 5

Impacts of improving the treatment effect on the infected rate (a) and the rate of drug resistance (b). The blue box plots are based on the baseline parameters, the green ones are based on increase lifespan, and red ones are based on both increase lifespan and reduce infectivity. The continuous blue, green and red curves are their median lines respectively.

3.2.4. Sensitivity Analysis

In previous studies, we unexpectedly found that merely increasing the treatment rate to the target rate (90% treatment rate) could not lead to a significant improvement in HIV epidemic (Figure 5a). In order to find the reason, we did some sensitivity analysis. We find that only when shortening each infectious stage’s life span to less than 25% of its baseline values can the epidemic be controlled. Details can be found in Appendix A.5. In the process of new drug development, if we cannot effectively reduce the infection rate while continuing to pursue prolonging the life-span of patients after treatment, the epidemic will keep rising as a consequence. We speculate that this may be one of the reasons that the effective reproduction has been increasing since 2012 (Figure 2). In other words, it is difficult to control the epidemic only by developing new drugs in a short period. So, what are the key factors to control the HIV/AIDS epidemic?

To answer this question, we conducted a sensitivity analysis on the number of sexual partners in different populations. An interesting fact is that changing the number of heterosexual people’s sexual partners has little impact on the epidemic, and it almost can be negligible while comparing with the two subpopulations of homosexual sexual behavior who contributed a significant effect to the epidemic. As long as the number of bisexual men’s annual sex partners is controlled below 6, and meanwhile, the number of homosexual men’s annual sex partners is controlled below 4, the epidemic can be controlled. Details can be found in Appendix A.6.

4. Discussion and Conclusions

Generally speaking, predictive models for infectious diseases generally cannot avoid the problem of parameter estimation. After all, not all parameter values are easy to obtain. Even if some foreign literature mentions the value range of certain parameters, there may still exist differences due to different races. To some degree, parameter values determine the credibility of predictions. Facing the complex history of HIV treatment in China, different historical stages may have different values even for the same parameter. We innovatively carried out data mining on virus gene sequence detected from the infected people, so as to fit some parameter values which are difficult to obtain by social surveys or laboratory tests, making our prediction results more reliable. However, other non-Bayesian research methods are also good choices for the mathematical modeling of infectious diseases with a large number of reported data [35].

The idea of combining gene sequences with dynamic models came from the absence of monitoring data when doing model fitting. Unlike novel coronavirus patients, HIV-infected people do not have obvious infection symptoms. Many HIV-infected people do not know that they have been infected for a long time. This prevents the government from obtaining information about people living with HIV as early as possible. Moreover, the Chinese government’s screening efforts for HIV-infected people are far less than those for individuals infected with novel coronavirus. As a result, recorded HIV infection data, if available, may not be very accurate and may not be suitable for parameter fitting. In Anhui province, which we studied, there is not much surveillance data on HIV/AIDS that can be used. The only reliable data available for use are 340 genetic sequences of HIV-infected people taken before treatment. HIV is an RNA virus. This gives it the ability to mutate faster, and these mutations show up in the entire genetic sequence of the virus from generation to generation. In other words, the evolutionary history of HIV is hidden in its genetic sequences. Since both gene sequences and dynamic models can describe the transmission dynamics of viruses, it should be possible to combine them for the purpose of learning from each other. This is our motivation to finish this manuscript.

What is more, most dynamic models generally only consider heterosexuality and homosexuality when classifying the target population, and they will not further subdivide homosexual people. Our study divided the MSM population into homosexual men and bisexual men to analyze the importance of bisexual men in the spread of HIV/AIDS. These bisexual men may be biological and psychological bisexuals; also, it may be that homosexuals are forced to marry women under the pressure of public opinion to cover up their particular sexual orientation. We generally think that such people are a bridge that links HIV/AIDS from high-risk groups to the general population. Thus, cutting off this bridge will slow down the HIV epidemic. However, the result of our model shows an opposite result: if bisexual men with the unchanged number of sexual partners disconnect from men, the epidemic will slow down a lot; otherwise, it will lead to an increase in the epidemic. For more details, please read Appendix A.7.

In “China’s 13th Five-Year Plan for Containment and Prevention of AIDS”, China advocates that HIV/AIDS prevention and treatment should achieve three 90%: the proportion of infected people and patients who have been diagnosed and know about their infection status should reach 90% or more; the proportion of infected people and patients who meet the treatment conditions receiving antiviral treatment should be more than 90%; the success treatment rate of infected people and patients receiving antiviral treatment should be more than 90%. These three 90% involve screening, treatment and drug development. According to our findings, if no more effective drugs are available, simply increasing the treatment rate will cause a slight decrease in the number of new infections and the number of people with primary drug resistance. Moreover, the epidemic is still on the rise; that is to say, the increase of treatment rate cannot effectively control the epidemic, and the prevention and control effect is not ideal. A similar phenomenon was also found in the study [36] by Lou Jie et al. A plausible explanation is that after treatment, the life expectancy of HIV-infected patients is extended, but the current drug’s effect is not so good at reducing the infection rate, thus creating a greater possibility of transmission instead. In the meantime, free treatments of HIV have not yet covered the entire course of patients, so once drug resistance develops, the treatment led to failure.

Therefore, after the implementation of the “90–90–90” strategic measures, continuous improvements of the treatment rate did not bring the HIV epidemic under control but inversely caused more drug resistance. The Chinese government needs to intervene in various aspects of different groups of people in combination with other perspectives in order to curb the current domestic AIDS epidemic. The sensitivity analysis on the two critical parameters, infection rate after treatment and progress rate of disease, indicates neither one is the most critical factor to impact the epidemic. It would be easier to control the epidemic by reducing the number of sexual partners of the two subpopulations in the MSM group through publicity and other means. In addition, it is worth mentioning that considering the absolute leading role of homosexual people in the epidemic at the present stage, we should tilt most of the resources to this group so that we can benefit from spending money wisely.

In short, gene sequences can tell us the history of HIV/AIDS spreading and the dynamic model can tell us its future. Combining the two together is an innovative approach, especially for epidemics where reliable surveillance data are lacking. In addition, we believe that this new method is not only suitable for HIV but also for other RNA viruses, such as novel coronavirus. Of course, our method is not perfect, since the influencing factors considered by the Bayesian model are relatively simple. However, at the very least, genetic sequence data can be used as a supplement to macro surveillance data for epidemic prediction.

Table A1

The initial population in 2015.

	Population	Heterosexual Women	Heterosexual Men	Homosexual Men	Bisexual Men
Stage
S		S^w= 24,347,402	S^mw= 24,682,596	S^mm= 449,877	S^b= 330,210
I_1		I1w=1200	I1mw=1216	I1mm=1857	I1b=1365
T_1		T1w=0	T1mw=0	T1mm=0	T1b=0
R_1		R1w=5	R1mw=5	R1mm=5	R1b=5
I_2		I2w=335	I2mw=340	I2mm=527	I2b=389
T_2		T2w=810	T2mw=820	T2mm=1240	T2b=925
R_2		R2w=60	R2mw=61	R2mm=85	R2b=67
I_3		I3w=335	I3mw=340	I3mm=527	I3b=389
T_3		T3w=810	T3mw=820	T3mm=1240	T3b=925
R_3		R3w=60	R3mw=61	R3mm=85	R3b=67

Table A2

Infection rate: .

Parameters	Description	Value (95% CI)	Source
β1I	infection rate of wild strain carriers in stage 1	0.0428 (0.0160, 0.6590)	[15,37]
β2I	infection rate of wild strain carriers in stage 2	0.0500 (0.0200, 0.0700)	[37]
β3I	infection rate of wild strain carriers in stage 3	0.1000 (0.0500, 0.1400)	[37]
β1T	infection rate of wild strain carriers after HAART in stage 1	0.0282 (0.0106, 0.4349)	[14]
β2T	infection rate of wild strain carriers after HAART in stage 2	0.0330 (0.0132, 0.0462)	[14]
β3T	infection rate of wild strain carriers after HAART in stage 3	0.0120 (0.0060, 0.0168)	[14,38]
β1R	infection rate of drug-resistant carriers in stage 1	ϵβ1I
β2R	infection rate of drug-resistant carriers in stage 2	ϵβ2I
β3R	infection rate of drug-resistant carriers in stage 3	ϵβ3I
ϵ	transmissibility ratio of drug-resistant to untreated wild strain carriers	fit

Table A3

Disease progression rate: .

Parameters	Description	Value: Years	Source
1/θ1I	course of wild strain carriers in stage 1	4	[14]
1/θ2I	course of wild strain carriers in stage 2	6	[14]
1/θ3I	course of wild strain carriers in stage 3	2	[14]
1/θ1R	course of drug-resistant carriers in stage 1	4	[14,39,40]
1/θ2R	course of drug-resistant carriers in stage 2	6	[14,39,40]
1/θ3R	course of drug-resistant carriers in stage 3	2	[14,39,40]

Table A4

Other parameters.

Parameters	Description	Value	Source
η	the percentage of patients who give up HAART	0.0602	[41,42]
γ	rate of drug resistance after HAART	0.0545	[42,43]
ρ	the percentage of patients who are no longer resistant	0.2687	[42,44]
μ	natural mortality rate	0.0131	[45]
k	ratio of female to other and male to other transmission	0.5263	[36]

Table A5

Other fitted parameters.

Period	Parameter	Value/Function
I	a^w	509,832
	a^mw(t)	542,781 − a^mm(t) − a^b(t)
	a^mm(t)	1432.6 × t + 10,527
	a^b(t)	1607 × t + 10,431
	c^mm(t)	−0.036928 × t + 2.8484
	c^b(t)	−0.040123 × t + 3.1024
	c^w	1.5228
	c^mw	1.5659
	v	0.70368
II	a^w	613,031
	a^mw(t)	685,421 − a^mm(t) − a^b(t)
	a^mm(t)	1622.4 × t + 12,382
	a^b(t)	1398.7 × t + 10,591
	c^mm(t)	−0.042873 × t + 3.4122
	c^b(t)	0.092656 × t + 3.4909
	c^w	1.4985
	c^mw	1.5066
	v(t)	−0.0098732 × t + 0.59836
	δ_3(t)	0.03733 × t + 0.35124
III	a^w	449,161
	a^mw(t)	493,770 − a^mm(t) − a^b(t)
	a^mm(t)	1455.8 × t + 5545.3
	a^b(t)	2079.9 × t + 5093.9
	c^mm(t)	1/(0.61892 + exp(−2.6515 × t + 3.1585)) + 3.5015
	c^b(t)	1/(0.54199 + exp(−2.7309 × t + 3.2935)) + 3.4814
	c^w	1.5362
	c^mw	1.5631
	v(t)	0.07004 × t + 0.51912
	δ_2(t)	1/(1.6826 + exp(−1.809 × t + 0.53915))+ 0.2192
	δ_3(t)	1/(1.6585 + exp(−1.8287 × t + 0.55639)) + 0.19956

Table A6

New infection population under different measures.

	Years	2021	2022	2023	2024	2025
Measures
Baseline parameter (95% CI)		3194	3499	3831	4207	4619
		(2829, 3560)	(3031, 3970)	(3257, 4434)	(3511, 4955)	(3780, 5541)
(1) The treatment rate is increased by 10% (95% CI)		3150	3432	3752	4102	4497
		(2787, 3510)	(2981, 3884)	(3193, 4325)	(3424, 4820)	(3670, 5378)
(2) The treatment rate is increased by 20% (95% CI)		3007	3258	3554	3875	4234
		(2678, 3360)	(2855, 3714)	(3049, 4122)	(3269, 4585)	(3502, 5110)
(3) The withdrawal rate is reduced by 10% (95% CI)		3128	3416	3738	4094	4495
		(2787, 3500)	(2982, 3893)	(3201, 4346)	(3450, 4861)	(3708, 5437)
(4) The withdrawal rate is reduced by 20% (95% CI)		3069	3336	3656	4002	4392
		(2740, 3448)	(2983, 3937)	(3143, 4279)	(3381, 4781)	(3628, 5344)
(5) The number of sexual partners is reduced by 10% (95% CI)		2781	2991	3223	3474	3755
		(2474, 3113)	(2613, 3406)	(2768, 3740)	(2940, 4106)	(3123, 4502)
(6) The number of sexual partners is reduced by 20% (95% CI)		2326	2456	2602	2757	2925
		(2062, 2601)	(2142, 2788)	(2232, 2999)	(2332, 3228)	(2435, 3477)
(1)(3)(5) combination (95% CI)		2780	2927	3153	3401	3677
		(2502, 3145)	(2608, 3424)	(2766, 2744)	(2931, 4094)	(3104, 4494)
(2)(4)(6) combination (95% CI)		2309	2416	2540	2680	2831
		(2056, 2583)	(2119, 2743)	(2200, 2929)	(2291, 3138)	(2384, 3373)

Table A7

Primary drug resistance population under different measures.

	Years	2021	2022	2023	2024	2025
Measures
Baseline parameter (95% CI)		209	235	262	292	323
		(188, 230)	(209, 261)	(231, 294)	(253, 331)	(276, 371)
(1) The treatment rate is increased by 10% (95% CI)		208	234	261	291	321
		(187, 229)	(209, 260)	(231, 295)	(255, 333)	(277, 374)
(2) The treatment rate is increased by 20% (95% CI)		201	228	257	287	318
		(181, 222)	(204, 255)	(227, 189)	(251, 327)	(274, 368)
(3) The withdrawal rate is reduced by 10% (95% CI)		206	231	258	287	318
		(185, 227)	(207, 258)	(229, 291)	(251, 327)	(274, 367)
(4) The withdrawal rate is reduced by 20% (95% CI)		202	228	254	283	313
		(183, 224)	(205, 255)	(226, 288)	(249, 325)	(272, 365)
(5) The number of sexual partners is reduced by 10% (95% CI)		184	205	227	250	273
		(168, 199)	(183, 228)	(201, 256)	(218, 284)	(236, 315)
(6) The number of sexual partners is reduced by 20% (95% CI)		155	172	188	205	221
		(139, 172)	(153, 192)	(166, 212)	(179, 232)	(190, 254)
(1)(3)(5) combination (95% CI)		176	197	219	241	264
		(157, 194)	(175, 219)	(192, 246)	(209, 275)	(226, 305)
(2)(4)(6) combination (95% CI)		154	170	184	201	219
		(137, 174)	(151, 190)	(164, 210)	(178, 234)	(189, 250)