MODELING AND DYNAMICS OF HIV TRANSMISSION AMONG HIGH-RISK GROUPS IN GUANGZHOU CITY, CHINA∗

In this paper, a multicompartmental model is formulated to study how HIV is transmitted among different HIV high-risk groups, including MSM (men who have sex with men), FRs (foreigner residents), FSWs (female sex workers), and IDUs (injection drug users). The explicit expression for the basic reproduction number is obtained via the next generation matrix approach. We show that the disease free equilibrium is locally as well as globally asymptotically stable (the disease goes to extinction) when the basic reproduction number is less than unity, and the disease is always present when the basic reproduction number is larger than unity. As an illustration of our theoretical results, we conduct numerical simulations. We also conduct a case study where model parameters are estimated from the demographic and epidemiological data from Guangzhou. Using the parameter estimates, we predict the HIV/AIDS trend for each high-risk group. Furthermore, our study suggests that reducing the transmission routes of the disease and increasing condom use will be useful for control of HIV transmission.


Introduction
Human Immunodeficiency Virus (HIV) is a sexually transmitted disease (STD) that causes acquired immunodeficiency syndrome (AIDS). AIDS causes gradual failure of the immune system, leaving victims vulnerable to a variety of infections. Since the discovery of HIV-1 in the early 1980s, the disease has spread in successive waves to most regions around the globe. In the world, about 36.9 million people are infected with HIV, and an estimated 1.1 million people died due to AIDS in 2017 [44]. It is one of the top ten infectious diseases and a leading cause of death in mainland China, as reported by the China Center for Disease Control and Prevention (China, CDC). The cumulative total number of reported HIV infection was 89,067 as of December 2004 [20], a figure that increased to 758, 610 as of December 2017 [28].
On the one hand, the main factors by which HIV/AIDS has been spreading in China include needle-sharing among intravenous drug users and unsafe sexual activities in the past decade. By 2002, HIV was present amongst IDUs in all Chinese inland provinces. It is believed that IDUs may have been the core source for all later sub-epidemics in China [42]. However, in recent years, the largest number of HIV/AIDS cases have been among MSM, followed by FSWs and IDUs [12]. The new HIV/AIDs infected cases among MSM and IDUs are respective 34,358 (25.5% of new total HIV/AIDS cases) and 4,296 (3.2% of new total HIV/AIDS cases) in 2017, China [28]. Hence, the Chinese government has strengthened its intervention and control efforts to MSM group, developed national working policies and guidelines on HIV/AIDS prevention and control among MSM, and convened national technical workshops on comprehensive HIV/AIDS prevention interventions among MSM.
On the other hand, statistics in recent years indicate that there are three main reasons for the continued surge in the number of foreigners living with HIV in China: the removal of restrictions on the entry of AIDS patients by the Chinese government. Unlike other parts of China, Guangzhou benefits from the reform and opening policy and has a close linkage with the rest of the world. Guangzhou ranks the third highest in the nation for the number of FRs, increasing from 19,000 in 2008 to 49,798 in 2017 [43]. The number of HIV infection cases among FRs increased significantly [17] after the Chinese government canceled entry requirements for restrictions on AIDS and HIV infected foreigners in 2010 [4,45]. Abroad HIV/AIDS cases found in Guangzhou are mainly young and middle-aged people, and sexual behavior is the primary route of HIV infection in these FRs HIV/AIDS cases [1]. Above all, we conclude that FRs play a crucial role in HIV transmission in Guangzhou. Based on these data and comparisons, HIV high-risk groups (groups of people with high-risk to infect with HIV) in Guangzhou include FRs, MSM, FSWs, and IDUs.
In recent time, mathematical models have been applied to the study of the spread of HIV/AIDS in the high-risk group. Hethcote etc [18] (1992) first formulated a model for MSM, studied biosynthetic properties and estimated parameters and incidences for MSM in San Francisco. Tan etc [22] (1993) proposed a stochastic model to present HIV/AIDS epidemic and infection in MSM. Tan etc [33] (1995) considered MSM groups, studied the effects of the randomness of risk factors on the HIV epidemic. Soon afterward, Tan etc [34] (1996) formulated a specific model for the HIV epidemic and the effects of age and race on MSM groups. In addition, a state-space model of the HIV epidemic in MSM was established in [35]. Zindoga etc [51] (2009) established a model with MSM (homosexuals and bisexuals) HIV infection. The model in this paper considered the emigration influence on the spread of HIV in MSM. Xu etc [46] (2011) presented a model concerning sexual transmission of HIV/AIDS among FSWs and MSM in Jiangsu province, China. Sun etc [29] (2013) proposed a mathematical model for HIV/AIDS epidemics among MSM in China that consider antiviral therapy. Zhang etc [52] (2016) developed a mathematical model on the transmission dynamics of HIV among MSM, FSWs, and IDUs groups. Li etc [23] (2019) established a mathematical model to study the persistence of HIV-1 spreading in MSM population in China. However the existing mathematical models only considered HIV transmission among the above three groups. How HIV transmission among high-risk groups when FRs group is presented remains unclear and is the subject of this study.
Based on the simplified flow chart of potential bridges for HIV infection to MSM, Guangzhou, China (Figure 1 in [20]), we give a more realistic and reasonable transmission flow diagram in Figure 1. The purpose of this study is to extend the existing mathematical models by including FRs group in HIV high-risk groups.
Moreover, we consider the interconnection by sexual behaviors among these groups. Hence, the total population will be divided into four groups: (1) MSM, (2) FRs, (3) FSWs, and (4) IDUs. Then sensitivity analysis of the basic reproduction numbers will be carried out in terms of the model parameters, from which effective control measures will be discussed for the spread of HIV among HIV high-risk groups in Guangzhou, China. (red imaginary line), the incidence rate from I3 to S1 is (orange imaginary line), the incidence rates from I2 to S1, S3 are α 21 I 2 N 1 and α 23 I 2 N 3 (purple imaginary line). Red solid lines show the moving rate from S1, S3 to S4 and green solid lines show the moving rate from S4 to other susceptible groups.
The paper is structured as follows. The model is formulated in Section 2. Theoretic analyses of model (2.1) (the basic reproduction ratio and the global asymptotical stability of disease-free equilibrium) are presented in Section 3 and Section 4. In Section 5, we prove the permanence of the disease. Numerical simulations to extend found analytical results are provided in Section 6. In Section 7, a case study is applied to Guangzhou, China. Conclusion and discussion round up the paper in Section 8.

Model development
The model sub-divides the total population at time t into the following groups: susceptible groups S i (t), infected groups I i (t), AIDS groups (the people in this group refer to first instance of AIDS) A i (t) (i = 1, 2, 3, 4), where 1,2,3, and 4 represent MSM, FRs, FSWs, and IDUs groups, respectively. In order to construct the corresponding model, we make the following assumptions: (1) The susceptible among MSM, FSWs, and IDUs are classified by their behaviors. That is, if a susceptible individual drugs with others by sharing needles at time t, then he/she belongs to IDUs group S 4 (t); (2) Since bisexuality is difficult to distinguish among MSM, we assume that MSM group includes bisexuality in our model, which means that there is sexual transmission between MSM and FSWs; (3) AIDS cases in the MSM and FSWs are assumed to be sexually inactive, and so do not influence the dynamics of the virus thus are not distinguished; (4) We do not consider sexual transmission of the disease between IDUs and other three groups, and just consider the infection of FRs to other groups in our model. Since MSM, IDUs, and FSWs in this paper respective refer to the domestic men who have sex with men, injection drug users, and female sex workers, we ignore the movement between FRs and IDUs.
Flow diagram of HIV infected among high-risk groups is demonstrated in Figure  1. The model is expressed through the following system of ordinary differential equations: represents the total number of the ith group. Some parameters are given by the similar form to [52] as follows where d 11 , d 13 44 , and r i , Λ i (i = 1, 2, 3, 4) are described in Table 8 detailedly, µ and δ represent natural death rate and AIDS related death rate, respectively.

The basic reproduction number R 0
The objective of this section is to calculate the basic reproduction number of system (2.1). To this purpose, we first denote the disease free equilibrium E 0 =(S 0 1 , 0, 0, S 0 2 , 0, 0, S 0 3 , 0, 0, S 0 4 , 0, 0), then substituting E 0 into system (2.1), yields It follows that For convenience, we denoteĒ 0 =(I 1 (t), ). According to the concepts of the next generation matrix [5,53], F and V are given as follows Then the matrices F (related to new infections) and V (related to transfers between classes) are given by Hence, the next generation matrix is its eigenvalues satisfy the following characteristic equation . Thus, we obtain that the basic reproduction number R 0 as follows where R 1 0 , R 2 0 respective represent the basic reproduction number of PRA and IDUs, R 3 0 represents the basic reproduction number of MSM and FSWs, and they can be expressed as follows For convenience analysis later, we denote Then we obtain that Next, we present the relationship between R 3 0 , R * 0 and 1, respectively. Define It follows from (3.2) that y(R 3 0 ) = 0 and y(1) = 1 − R * 0 . One obtains the following equivalence relations:

Global stability of E 0
The main focus of this section is to analyze the local and global behavior of the disease-free equilibrium E 0 of system (2.1). First, we state and prove the result about local asymptotic stability for the disease-free equilibrium E 0 . The 12 × 12 Jacobian matrix J(E 0 ) can be presented as follows where where The eigenvalues of matrix (4.1) are determined by those of matrices (4.2), (4.4), and (4.5), it suffices to prove that all eigenvalues of the above matrices must have negative real parts when To achieve this goal, we have to show the following Lemma Proof. Part (1) follows directly from K 4×4 is real column diagonal dominant matrix and all the diagonal entries are negative. To prove part (2), the characteristic equation of H 2×2 is expressed as Then we have a 1 < 0, a 2 < 0, which implies all of eigenvalues of H 2×2 have negative real parts. On the contrary, the matrix H 2×2 has at least one eigenvalue with positive real part as R 1 0 > 1, which proves part (2). For part (3), the proof of this result is quite similar to that given earlier for part (2) and so is omitted.
Using these analysis together with Lemma 4.1 we obtain that the following theorem.
To complete the proof of the global stability of E 0 , we first show the following Theorem.
According to the standard comparison principle, there exists t 1 > 0 such that for any ϵ > 0, Proof. We divide our proof in three steps. As to continuous and bounded function Based on Theorem 4.2, we obtain that any solution of system (2.1), By the fluctuation lemma (Lemma 4.2 in [16]), we know there exists a sequence {t n }, with t n → ∞ (n → ∞), such that Substituting the sequence {t n } into the infected compartments of system (2.1) and taking superior limit for both sides of those equations, we obtain that It follows from the inequality (4.10) that that is to say that lim t→∞ I 2 (t) = 0. Thus, we obtain lim t→∞ A 2 (t) = 0 from the inequality (4.10). The same proof still goes for lim t→∞ A 4 (t) = lim t→∞ I 4 (t) = 0 when R 2 0 < 1. The next thing to do this proof is show that if R * 0 < 1, then I j (t) → 0, A j (t) → 0, j = 1, 3. We assume I ∞ 1 > 0, then it follows from inequalities (4.9), (4.11) that The inequality (4.13) holds for R * 0 > 1, which contradicts with R * 0 < 1. Hence, , it may be concluded that the limit system of

Permanence of the disease
In this section, to show permanence of the disease of system (2.1), we give the following Theroem.
Proof. By the similar arguments as those in the paper [36]. We first introduce some notations which will be used throughout the proof process.
We will show that system (2.1) is uniformly strong persistent with respect to (X 1 , X 2 ).
Denote Ω = ∪ , where ω (Z(0)) is the omega limit set of solution . Thus the disease-free equilibrium E 0 is the unique equilibrium of system (2.1) in M ∂ , it is easy to see that E 0 is locally asymptotically stable, which implies that E 0 is globally asymptotically stable (since system (5.2) is a linear system). Therefore, we have Ω = {E 0 }, and E 0 is covering of Ω, which is isolated and is acyclic (since there exists no solution in M ∂ which links E 0 to itself). Finally, the proof will be done if we show E 0 is a week repeller for X 1 i.e.
From Theorem 4.4 together with Theorem 5.1, we can claim that the basic reproduction number R 0 is a threshold parameter, which determines the outcome of disease, namely if R 0 < 1, the disease-free equilibrium E 0 is globally asymptotically stable (the disease dies out). While R 0 > 1, the disease will permanence.

Numerical simulations
In this section, we implement numerical simulation to confirm and extend the analytic results, and illustrate the effect of the basic reproduction numbers on the disease dynamics. For this pattern, we choose r i = 1.25 and other parameters are the same as Table 8. The disease-free equilibrium E 0 is a unique stable state when R 0 < 1(theoretical result is obtained in Theorem 4.4), which is shown in Figure  2. The simulations are shown in Figure 3 extend our analytical result presented in Theorem 5.1. While we do not have results on the existence and stability of endemic equilibrium of the system (2.1), our numerical studies in Figure 3 indicate that endemic equilibrium may exist. It means that HIV/AIDS will never be eliminated in Guangzhou, China.  Table 8, R0 < 1.   Table 8, R 1 0 > 1, R 2 0 > 1, and R 3 0 > 1.

A case study
In this section, we use system (2.1) to predict the potential HIV transmission among HIV high risk groups in Guangzhou, China. Firstly, we collect and settle the data from Guangzhou. Secondly, we estimate the parameters of system (2.1) by the data fitting. Thirdly, we give some predictions for the HIV transmission among HIV high risk groups in Guangzhou. Last but not least, we study the effect of control measure with different interventions ( condom use, etc.). We obtain the following data of HIV/AIDS cases in Guangzhou from some papers and websites (see Table  1-6):   Other routes  2006  11  9  --2007  5  25  2  5  2008  6  25  20  3  2009  3  31  20  3  2010  7  24  20  8  Total  54  121  41  6  Proportion 16.2% 36.3% 25.2% -

Estimation of model parameters and initial values.
Here, we present a brief description of some parameters not already mentioned. The natural human death rate µ, is measured as the reciprocal of the average lifespan of human in Guangzhou. The value is selected based on data from the Guangzhou Statistics Bureau [15]. The parameters, 1/r 1 , 1/r 2 , 1/r 3 , and 1/r 4 , whose give the respective average time HIV virus carriers convert to AIDS patients in the four HIV infected human groups. They are estimated from paper [6] to be in the range of 5-12 years. Mean value of HIV-related death of 0.7114 years −1 has also been cited [57]. The key parameters, condom use rates d 1i , d 3i and d 44 are obtained from literature [19,21,25]. Those paper all focus on statistic and Mate analysis based on the real demographic and epidemiological data in Guangzhou city, hence our parameter values are reasonable and credible. Similarly, the infected ratio of different groups c 1i , c 2i , c 3i , and c 44 are estimated on the basis of the published medical research articles [3,40,54,55,58], respectively. Table 8 lists all parameters, their definitions, and their ranges. The simulation of system (7.1) shown in Figure  5 and the references used to determine the parameter values. Moving rates from S 1 to S 4 -0.0234 Fit p 34 Moving rates from S 3 to S 4 -0.011 Fit p 41 Moving rates from S 4 to S 1 -0.014 Fit p 43 Moving rates from S 4 to S 1 -0.0132 Fit d 11 , d 13 Condom using ratio for MSM From Table 4, there is sharp rise of MSM infected HIV/AIDS in 2008, when 117 new cases are seen, meanwhile , the total number of HIV/AIDS cases is 1584 in Table 7. Besides, the number of HIV/AIDS among foreigner is 85 in 2008, Guangzhou in Table 6. The initial time t = 0 is chosen in the year 2008. From the statistical data of four district, Guangzhou in 2008 (see Table 1-6, Figure 4), we estimate among those HIV/AIDS cases, 32.31% are infected via injecting drug use, 57.69% are through heterosexual transmission. Based on the statistical data in the literatures [11,50], we obtain the ratio of men to women is almost 3.5:1, and the ratio of HIV case is almost 75% of HIV/AIDS cases. Then we estimate five compartments populations and HIV/AIDS cases in 2008 (t = 0), respectively.
MSM compartment: From Table 4, we obtain initial value I 1 (0) + A 1 (0) = 117.  [13,43]. From Table 6  To get the estimates of p 14 , p 34 , p 41 , p 43 , and Λ i (i = 1, 2, 3, 4), we drop out the terms involving I i , A i and S i ∼ = N i (HIV/AIDS cases account for a very small proportion of the high-risk population). Then we obtain the following system from system (2.1) We use system (7.1) to fit the total numbers corresponding to the IDUs, MSM and FRs data, respectively (see Fig. 5), where µ = 0.0063 (see Table 8).
Then we obtain the recruitment Λ i and the moving rates p 14 , p 34 , p 41 , p 43 as follows  Table 8, and obtain R 2 0 = 0.847 < 1, R 1 0 = 1.135 > 1 and R 3 0 = 1.607 > 1. In general, numerical estimation results show in Figure 5 and Figure 6 indicate that our model provides a relatively good match with the annual new reported infected HIV/AIDS data of MSM, FRs, FSWs, IDUs HIV/AIDS cases in Guangzhou from 2008 to 2017, respectively ( the statistical data in Table 4, 7, 6). Moreover, we use our model to predict the tend of four high-risk groups HIV/AIDS infected cases in Guangzhou, it can be seen that the HIV/AIDS infected cases of MSM, FRs and FSWs will increase in the next few  Table 8.     years. Specifically, the HIV/AIDS infected populations of MSM increase rapidly. However, the HIV/AIDS infected populations of injecting drug users will decrease sharply. These results show that the sexual transmission is the main route among HIV high-risk groups, whose are in line with literature [12].

Sensitivity analysis
In order to analyze the effects of the parameter values on the basic reproduction number of system (2.1), we perform sensitivity analysis by Latin square sampling and partial rank correlation coefficient (PRCC) methods [26]. In the absence of available data on the distribution functions, we choose a uniform distribution for all input parameters with the minimum and maximum values are shown in Table 8 and tested for significant PRCCs for all parameters of R 3 0 (PRCCs is given in Table  9). Figure 7 shows PRCCs values of the parameters against the basic reproduction number, which indicates that the parameters c 11 , d 11 and r 1 of the model (2.1) have significant impact on the basic reproduction number. These results indicate that decreasing the transmission rate of sexual or increasing the condom use proportion are two effective measures to reduce the basic reproduction number, that is, the development of highly condom use is a critical control measure. So, studying the effects of changing these two parameters on the basic reproduction number is of significant applied value.   Figure 8 (a). Similarly, in Figure 8 (b), it is shown that R 1 0 and R 2 0 are increasing functions of c 22 and c 44 , respectively. R 1 0 and R 2 0 are decreasing functions of d 22 and d 44 , respectively. Furthermore, from Figure  8 (c), (d) we can see that R 2 0 nearly less than 1, it means that the population of infected HIV among IDUs will decrease and eliminate ultimately, which implies that injecting drugs is not main route of HIV transmission any more. In Figure 9, it can be seen that R 3 0 almost more than 1, this result reveals that HIV infected cases will increase sharply among MSM (homosexual/bisexual) through sexual transmission and this trends will not be improved unless some control measures are taken.
We are also interested in the impact of the recruitment Λ i (i = 1, 2, 3, 4) on the total HIV/AIDS cases. For this purpose, we examine the variation of I T (t) + A T (t) with parameters Λ i in Figure 10, here I T (t) = I 1 (t) + I 2 (t) + I 3 (t) + I 4 (t), A T (t) = A 1 (t) + A 2 (t) + A 3 (t) + A 4 (t) and Λ i (i = 1, 2, 3, 4) being all varied in the interval [0, 4000]. From Figure 10 (a), it is revealed that I T (t) + A T (t) is an increasing function with respect to Λ 1 and the total HIV/AIDS cases will be decrease dramatically when Λ 1 decreases. Furthermore, it follows from Figure 10 (b), (c) that both Λ 2 and Λ 3 affect the number of HIV/AIDS cases in a positive manner, but to lesser extent than Λ 1 . Moreover, it can be seen that Λ 4 is barely influenced on the trends of the total HIV/AIDS cases in Figure 10 (d). These results indicated that the leading cause of the increase in HIV/AIDS cases is the increase in the number of MSM, followed by FSWs, FRs. From a practical point of view, we claim that the intervention of the high-risk groups is a very effective measure to control the spread of HIV disease.   Table 8.   Table 8.      Table 8.
From the numerical fitting and prediction results in Figure 6, it can be observed that the HIV/AIDS cases among MSM, FSWs, and PRs increased while the HIV/AIDS cases among IDUs decreased year by year. By analyzing the impact of the recruitment of four high-risk groups on the total cases (see Figure 10), we found that the total HIV/AIDS cases is most affected by MSM, followed by FSWs, FRs, and IDUs. Based on the above numerical results, we conclude that the order of significance of the four high-risk groups in HIV transmission is MSM, FSWs, FRs, IDUs. These facts indicated that MSM, FSWs, and FRs should be the key targets for HIV control to prevent the spread of the disease from high-risk groups to the general population.   Table 8.

Intervention and control
Simulation results show that decreasing infection rate and increasing condom use can slow down the spread of HIV transmission (see in Figure 7) and the increase in MSM is primary factor of the increase in HIV/AIDS cases among high-risk group (see in Figure 10 ). Hence, we now analyze the infection rate and condom use influence HIV transmission among MSM, whose always have more risky sexual behavior and more likely to transmit HIV virus than other groups. In Figure 11 (a) − (c), it can be seen that increasing condom use when MSM have sex with other groups can reduce HIV transmission. In Figure 11 (d) − (f ), it is shown that decreasing infection rate between MSM and other groups can achieve the same goals. Hence, to control the disease among MSM group availably, the intervention measure for whole high-risk population behaviors must be taken.

Conclusion and Discussion
It is known that HIV high-risk groups play an important role in the spread of HIV epidemic. Hence, compartment models for HIV transmission among highrisk groups (IDUs, FSWs, MSM) have been discussed by many researchers, but these models mainly considered a subset of high risk groups. Yang formulated only two groups (FSWs and senior male clients) [48]. Sun analyzed an HIV/AIDS epidemic only among men who have sex with men [30]. Zhang and Zhou established one model including IDUs, FSWs, MSM, Senior male clients [52]. Although four high-risk groups are considered, the interconnection by sexual behaviors among these groups is not considered. Compared to existing compartment models on HIV infection [20,[30][31][32]48], the proposed model here studies potential bridges for HIV transmission among high risk groups. Specifically, the foreigner residents are considered.
Our main concerns are studying the potential bridge for HIV transmission among HIV high-risk groups (MSM, SFWs, IDUs, FRs) and their HIV/AIDS epidemic tendency. We have studied the qualitative behavior of the model. It has been shown that the dynamics of this model are simple i.e., the disease-free equilibrium is globally asymptotically stable when the basic reproduction number R 0 < 1, and the infection equilibrium will permanence when R 0 > 1. Furthermore, numerical analyses of our model based on data in Guangzhou show the following facts: (1) Foreigner residents play an important role in HIV transmission, and pose a big challenge for HIV control. (the HIV/AIDS cases among PRs in Guangzhou increased rapidly after 2008, see Table 6 and Figure 6 (c)), (2) Homosexual transmission among MSM has replaced injecting drugs as the main route of HIV transmission, which is in lines with reports [2,7,14,38], (3) The leading causes of HIV/AIDS are MSM, FSWs, and PRs and these high-risk groups should be the key targets for HIV control to prevent the spread of the disease to the general population. Sensitivity analyses of the basic reproduction number on various model parameters indicate that it is important, indeed crucial, to increase the condom use rate for the HIV high-risk groups to control the transmission of HIV in Guangzhou, China. Moreover, reduction in any one or more routes of disease transmission will be useful, although reduction in transmission of the high-risk groups will be more effective than other intervention measures.
From a practical point of view, the model established in our paper can be used to comprehend the transmission behaviors of the disease among HIV high-risk groups in Guangzhou and to predict the disease trends in a few years, which can help Guangzhou's and China's departments of health to implement preventive interventions in the HIV high-risk groups.