library(magrittr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

1 Introduction:

HIV may have reached humans from one species of chimpanzees in central Africa as early as the late 1800s. This disease originated in Africa and later spread to other parts of the world (CDC, 2022). This virus targets the human immune system and is primarily spread through sexual contact with an infected individual. As the immune system weakens, the body becomes increasingly vulnerable to any diseases. Although there is no cure for HIV, antiretroviral therapy (ART) can help manage HIV by controlling the infection and preventing disease progression. Without ART, HIV will lead to AIDS, which will result in death (CDC, 2022). There are three stages of HIV infection: 1) acute phase, typically occurring 2 to 6 weeks after infection with HIV, is the earliest stage of HIV infection; 2) chronic HIV infection, a prolonged-10 years or longer incubation period , is the second stage of the infection; 3) without ART, patients will inevitably enter AIDS, the final and most severe stage (Department of Health. Victoria, 2015). Without treatment, people with AIDS typically survive about 3 years (HIVinfo,2021). As long as people get infected with HIV and with no treatment, individuals stay infected and infectious throughout their life, and remain in contact with the susceptible population until they die in the AIDS stage. Therefore, we decided to use the SI model to simulate HIV patients who were infected without using any treatment.

There are several modes within sexual transmission’s classification. For example, Simon and Ho (2006) that heterosexual transmission is the main mode of transmission of HIV, accounting for approximately 85% of all HIV-1 infections. Penile-vaginal transmission from males to females appears to be 2–3 times more efficient than the reverse, with evidence suggesting that initial sexual intercourse for females may carry particularly high transmission risks (Morrison, 2001). Similarly, male-male couples face elevated risks of HIV transmission through MSM (men who have sex with men) encounters. The CDC (2022) reported that anal sex is the riskiest type of sex for getting or transmitting HIV. According to Powers’ research in 2008, penile-anal transmission can occur as frequently as 1 in 3 exposures, while penile-vaginal transmission was reported at 1 in 10 exposures. The high risk of MSM transmission drove us to focus on this particular population-male-male couples. From 2007, the MSM transmission became significant in China. The proportion of new infections caused by transmission among MSM increased dramatically from 3.4% in 2007 to 25.5% in 2017 (Citation need). The HIV epidemic among MSM raises serious concerns in Chongqing, one of the four municipalities in China: Chongqing has the highest HIV prevalence among MSM in the nation with a rate of 13.8% and ranging from 10% to 23% between 2001 and 2018 (Citation need). Therefore, we narrowed down our population to male who lived in Chongqing and had sex with men in the past.

##explain the population we chose in the second paragraph of the introduction

2 Simulations

Given that we are studying HIV without the intervention of ART treatment, all models consist of compartments for Susceptible, Infectious (referring to individuals infected with HIV, including those who progress to AIDS), and Death (referring to death due to HIV or AIDS). In the upcoming simulation, we will simulate the development of HIV transmission among MSM in Chongqing from 2011 to 2018, in order to compare it with the real HIV infection situation caused by MSM transmission during this period in Chongqing. Therefore, our N represents the estimated homosexual population in Chongqing in 2011, and I represents the estimated population infected by MSM transmission in that year to be the initial infectious population.

The three models presented later include one deterministic and two stochastic models.

Notation Definition:

  1. lambda_t: Probability of individuals transitioning from the Susceptible (S) to the Infectious (I) stage, or the probability of each susceptible individual having sexual contact and then getting infected with an infectious individual.

  2. ks: The number of sexual contacts per unit time. The data on the number of sexual contacts per unit time (month) for the MSM population engaging in unprotected anal sexual conduct, ranging from 0.66 to 1.16 with a mean of 0.91 per month, is estimated based on academic sources. (Lu et. al 2020)

  3. pc_t: Probability of infection given contact between a susceptible and infectious individual after having an unprotected anal sexual conduct, ranging from 0.01 to 0.2, with a mean of 0.105 per month. (Wilton)

  4. Mortality rate: The mortality rate throughout the simulation process ranges from 0.057 to 0.082 per year for the first two years and from 0.082 to 0.143 per year for the subsequent years. We calculated a weighted mean mortality rate of 0.11 per year in both the deterministic and the first stochastic model (Li et. al 2016). Because we will be conducting monthly HIV progression simulations thereafter, we will divide the annual mortality rate by 12 to use later.

2.1 2a. Deterministic Model

This simulation models the development of HIV under the assumption that all parameters and initial conditions are fixed, and no consideration is given to the influence of randomness or uncertainty factors, which describes that happens on average in the population. So all the value of variables being chosed are mean values. In this section, I simulated the development of HIV infection over 96 months, plotted the data frame and graphs, and summarized the annual number of new HIV infections.

The equations are: \[\lambda_t=1-\exp\left\{- k_s\times pc_t\times \frac{I_{t-1}}{N}\right\}\] \[[S[t] = S[t - 1] - \lambda_t[t] \cdot S[t - 1] ]\] \[[I[t] = I[t - 1] + \lambda_t[t] \cdot S[t - 1] ]\] \[[\text{Dea}[t] = \text{Dea}[t - 1] + \text{round}( \text{MortalityRate} \times I[t - 1] ) ]\] \[\lambda_t*S[t-1]\] represents the number of people moving from the susceptible to infectious compartment

\[\texttt{round}(MortalityRate \times I_{t - 1})\] refers the number of people moving from the infectious to death compartment

HIV_deterministic_simulation <- function(N, S0, I0,MortalityRate, ks, pc_t, TotalTime) {
    
    S <- numeric(TotalTime + 1)
    I <- numeric(TotalTime + 1)
    Dea <- numeric(TotalTime + 1) 
    lambda_t <- numeric(TotalTime + 1) 

    I[1] <- I0
    S[1] <- S0
    Dea[1] <- 0
    lambda_t[1] <- 1 - exp((-ks*pc_t * I0) / N) 
    
    for (t in 2:(TotalTime + 1)) {
        lambda_t[t] <- 1 - exp(-ks *pc_t*I[t - 1] / N)
        S[t] <- S[t - 1] - lambda_t[t] * S[t - 1] 
        I[t] <- I[t - 1] + lambda_t[t] * S[t - 1] - round(MortalityRate * I[t - 1])
        Dea[t] <- Dea[t - 1] + round(MortalityRate * I[t - 1])
    }
    HIV_sim_df <- data.frame(Time = 1:TotalTime,
                              S = S[1:TotalTime],
                              I = I[1:TotalTime],
                              Dea = Dea[1:TotalTime],
                              lambda_t = lambda_t[1:TotalTime])
    return(HIV_sim_df)}

set.seed(257)
HIV_sim_df <- HIV_deterministic_simulation(16767, 14503, 2263, (0.11)/12, 0.91, 0.105, 96)

HIV_plot <- ggplot(data = HIV_sim_df, mapping = aes(x = Time)) +
        geom_line(aes(y = S, color = "Susceptible")) +
        geom_line(aes(y = I, color = "Infected")) +
        geom_line(aes(y = Dea, color = "Death")) +
        geom_line(aes(y = lambda_t, color = "Lambda_t")) +
        labs(title = "HIV Model Simulation", x = "Time (Months)", y = "Population") +
        scale_color_manual(values = c("Susceptible" = "red", "Infected" = "green", "Death" = "purple", "Lambda_t" = "orange"))

print(HIV_plot)

# Calculate newly infectious every 12 months
HIV_deterministic_simulation_summary <- function(N, S0, I0, MortalityRate, ks, pc_t, TotalTime) {
    HIV_sim_df <- HIV_deterministic_simulation(N, S0, I0, MortalityRate, ks, pc_t, TotalTime)
    HIV_sim_df$NewInfected <- HIV_sim_df$lambda_t * HIV_sim_df$S
    
    NewInfected <- tapply(HIV_sim_df$NewInfected, rep(1:ceiling(TotalTime / 12), each = 12, length.out = TotalTime), sum)
    summary_df <- data.frame(Year = seq(1, length(NewInfected)),
                              NewInfected = NewInfected)
    
    return(summary_df)
}

HIV_summary <- HIV_deterministic_simulation_summary(16767, 14503, 2263, (0.11)/12, 0.91, 0.105, 96)

# print(HIV_summary)

2.2 2b.Stochastic Simulations 1

In the first HIV stochastic model, we model the number of individuals moving between compartments at each time step according to a binomial distribution. Consequently, the simulation is stochastic, resulting in random fluctuations in the population moving between compartments, thereby better capturing the dynamics of disease progression in the real world.

After simulating the progress of HIV over 96 months, I first calculated the annual number of new HIV infections. Then, I conducted the same simulation 100 times with different seeds and plotted them on a single graph.

Notation Definition: 1. \[(I_{\text{SI}})^{(t)} \sim \text{Binomial}(S, \lambda t)\] This variable denotes the number of individuals who transition from the Infectious (S) compartment to the infected (I) compartment at time t. It is generated from a binomial distribution.

  1. \[(D_{\text{ID}})^{(t)} \sim \text{Binomial}(I, \ Mortalityrate)\] This variable denotes the number of deaths occurring in the infected (I) compartment at time t. It captures the individuals who succumb to HIV-related mortality. Similarly, it is generated from a binomial distribution.

  2. Binomial distribution: is a probability distribution that describes the number of successes in a fixed number of independent Bernoulli trials, each with the same probability of success. It is characterized by two parameters: the number of trials n and the probability of success p on each trial.

The differential equations we used: \[S[t] = S[t - 1] - I_{SI}[t] \\ I[t] = I[t - 1] + I_{SI}[t] - D_{ID}[t] \\ Dea[t] = Dea[t - 1] + D_{ID}[t]\]

HIV_stochastic_simulation <- function(N, S0, I0, MortalityRate, ks, pc_t, TotalTime) {
    # Create the compartments   
    S <- numeric(TotalTime + 1)
    I <- numeric(TotalTime + 1)
    Dea <- numeric(TotalTime + 1)
    lambda_t <- numeric(TotalTime + 1)
    I_SI <- numeric(TotalTime + 1)
    D_ID <- numeric(TotalTime + 1)

    # Initialize time step
    I[1] <- I0
    S[1] <- S0
    Dea[1] <- 0
    lambda_t[1] <- 1 - exp((-ks * pc_t* I0) / N) 
    
    for (t in 2:(TotalTime + 1)) {
        # Update transmission rate lambda_t
        lambda_t[t] <- 1 - exp(-ks * pc_t*I[t - 1] / N)
        
        # Update other compartments
        I_SI[t] <- rbinom(n = 1, size = S[t - 1], prob = lambda_t[t])
        D_ID[t] <- rbinom(n = 1, size = I[t - 1], prob = MortalityRate)

        S[t] <- S[t - 1] - I_SI[t]
        I[t] <- I[t - 1] + I_SI[t] - D_ID[t]
        Dea[t] <- Dea[t - 1] + D_ID[t]
    }
  
    HIV_stochastic_df <- data.frame(Time = 1:(TotalTime + 1),
                                 S = S,
                                 I = I,
                                 Dea = Dea,
                                 lambda_t = lambda_t,
                                 I_SI = I_SI,
                                 D_ID = D_ID,
                                 MortalityRate = MortalityRate)
  
    return(HIV_stochastic_df)
}

set.seed(256)

HIV_stochastic_df <- HIV_stochastic_simulation(16767, 14504, 2263, (0.11)/12, 0.91, 0.105, 96)
# print(HIV_stochastic_df)

HIV_stochastic_plot <- ggplot(data = HIV_stochastic_df, mapping = aes(x = Time)) +
  geom_line(aes(y = S, color = "Susceptible")) +
  geom_line(aes(y = I, color = "Infected")) +
  geom_line(aes(y = Dea, color = "Death")) +
  geom_line(aes(y = lambda_t, color = "Lambda_t")) +
  labs(title = "HIV Model Simulation", x = "Year", y = "Population") +
  scale_color_manual(values = c("Susceptible" = "red", "Infected" = "green", "Death" = "purple", "Lambda_t" = "orange"))

print(HIV_stochastic_plot)

set.seed(556)
HIV_stochastic_df_summary <- HIV_stochastic_df %>%
  mutate(Year = ifelse(Time <= 96, ((Time - 1) %/% 12) + 1, NA)) %>%
  filter(!is.na(Year)) %>%
  group_by(Year) %>%
  summarize(NewInfected = sum(I_SI))

# print(HIV_stochastic_df_summary)
SID_sim_list <- list()
for (i in 1:100) {
  set.seed(123 + i)  
  SID_sim_df <- HIV_stochastic_simulation(16767, 14504, 2263, (0.11)/12, 0.91, 0.105, 96)
  SID_sim_list[[i]] <- SID_sim_df
}

# Calculate summary statistics for the simulations
summary_df <- bind_rows(SID_sim_list) %>%
  group_by(Time) %>%
  summarize(
    S_lower = quantile(S, 0.25),
    S_upper = quantile(S, 0.75),
    I_lower = quantile(I, 0.25),
    I_upper = quantile(I, 0.75),
    D_lower = quantile(Dea, 0.25),
    D_upper = quantile(Dea, 0.75),
    S_mean = mean(S),
    I_mean = mean(I),
    D_mean = mean(Dea)
  )

# Plot the combined graph
combined_plot <- ggplot() +
  labs(title = "Combined HIV Model Simulation", x = "Time", y = "Population", color = "Group") +  
  scale_color_manual(values = c("Susceptible" = "#FF5733", "Infected" = "#33FF57", "Death" = "#3366FF")) 

for (i in 1:100) {
  sim_df_subset <- SID_sim_list[[i]]
  
  combined_plot <- combined_plot +
    geom_line(data = sim_df_subset, aes(x = Time, y = S, color = "Susceptible"), alpha = 0.05) + 
    geom_line(data = sim_df_subset, aes(x = Time, y = I, color = "Infected"), alpha = 0.05) + 
    geom_line(data = sim_df_subset, aes(x = Time, y = Dea, color = "Death"), alpha = 0.05) 
}

combined_plot <- combined_plot +
  geom_ribbon(data = summary_df, aes(x = Time, ymin = S_lower, ymax = S_upper, fill = "red"), alpha = 0.5, color = NA) +
  geom_ribbon(data = summary_df, aes(x = Time, ymin = I_lower, ymax = I_upper, fill = "green"), alpha = 0.5, color = NA) +
  geom_ribbon(data = summary_df, aes(x = Time, ymin = D_lower, ymax = D_upper, fill = "blue"), alpha = 0.5, color = NA) +
  geom_line(data = summary_df, aes(x = Time, y = S_mean, color = "Susceptible"), linetype = "dashed", size = 1) +  
  geom_line(data = summary_df, aes(x = Time, y = I_mean, color = "Infected"), linetype = "dashed", size = 1) + 
  geom_line(data = summary_df, aes(x = Time, y = D_mean, color = "Death"), linetype = "dashed", size = 1)  
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(combined_plot)

2.3 Stochastic simulation 2:

In this stochastic 2 model of HIV simulation, we introduced uncertainty by incorporating parameters as random variables. Unlike the stochastic 1 model, where parameters were fixed, here we assign probability distributions to parameters that determine movement between compartments. This approach aims to better reflect the uncertainty inherent in real-world scenarios, making the simulation more realistic. By treating parameters governing movement between compartments as random variables, we simulate a range of possible scenarios, capturing the variability observed in HIV transmission dynamics. We conducted simulations for 96 months and then calculated the annual number of new HIV infections.

Notation Definition: The only difference between this model and the previous one lies in the assignment of parameters:

\[\text{MortalityRate} = [0.057,0.082], \quad \text{for the first two years.} \\ \text{MortalityRate} = [0.082,0.143], \quad \text{for the years after.} \\ k_s \sim \text{Uniform}(0.66, 1.16) \quad (\text{Lu et al., 2020}) \\ k_s \sim \text{Uniform}(0.01, 0.2) \quad (\text{Wilson})\]

Here, we employ a uniform distribution to randomly sample values for these three parameters within the upper and lower bounds of the realistic data we are seeking. This randomization is performed within the loop to simulate various scenarios.

Uniform distribution, also known as rectangular distribution, is a continuous probability distribution where all outcomes are equally likely within a specified range. In this context, it ensures that parameter values are selected uniformly from the defined range, allowing for a comprehensive exploration of the parameter space and enhancing the robustness of the simulation results.

Differential Equations: Difference equations stay the same as in Stochastic 1

HIV_stochastic2_simulation <- function(N, S0, I0, TotalTime) {
  # Create the compartments   
  S <- numeric(TotalTime + 1)
  I <- numeric(TotalTime + 1)
  Dea <- numeric(TotalTime + 1)
  lambda_t <- numeric(TotalTime + 1)
  I_SI <- numeric(TotalTime + 1)
  D_ID <- numeric(TotalTime + 1)

  I[1] <- I0
  S[1] <- S0
  Dea[1] <- 0
  
  MortalityRate1 <- runif(1, (0.057)/12, (0.082)/12)
  MortalityRate2 <- runif(1, (0.082)/12, (0.142)/12)
  pc_t <- runif(1, 0.01, 0.2)  
  ks <- runif(1, 0.66, 1.66)
  
  for (t in 2:(TotalTime + 1)) {
    if (t <= 24) {
      lambda_t[t] <- 1 - exp(-ks * pc_t * I[t - 1] / N)
      I_SI[t] <- rbinom(n = 1, size = S[t - 1], prob = lambda_t[t])
      D_ID[t] <- rbinom(n = 1, size = I[t - 1], prob = MortalityRate1)
    } else {
      lambda_t[t] <- 1 - exp(-ks * pc_t * I[t - 1] / N)
      I_SI[t] <- rbinom(n = 1, size = S[t - 1], prob = lambda_t[t])
      D_ID[t] <- rbinom(n = 1, size = I[t - 1], prob = MortalityRate2)
    }
    
    S[t] <- S[t - 1] - I_SI[t]  
    I[t] <- I[t - 1] + I_SI[t] - D_ID[t]  
    Dea[t] <- Dea[t - 1] + D_ID[t]
  }
  
  # Create dataframe for simulation results
  HIV_stochastic2_df <- data.frame(Time = 1:(TotalTime + 1),
                                   S = S,
                                   I = I,
                                   Dea = Dea,
                                   lambda_t = lambda_t,
                                   I_SI = I_SI,
                                   D_ID = D_ID
  )
  
  return(HIV_stochastic2_df)
}
set.seed(156)  
sim_result <- HIV_stochastic2_simulation(16767, 14504, 2263, 96)
HIV_stochastic2_df <- sim_result

HIV_stochastic2_df_summary <- HIV_stochastic2_df %>%
  mutate(Year = ifelse(Time <= 96, ((Time - 1) %/% 12) + 1, NA)) %>%
  filter(!is.na(Year)) %>%
  group_by(Year) %>%
  summarize(NewInfected = sum(I_SI))

# print(HIV_stochastic2_df_summary)
SID_sim2_list <- list()
for (i in 1:100) {
  set.seed(123 + i)  
  SID_sim2_df <- HIV_stochastic2_simulation(16767, 14504, 2263, 96)
  SID_sim2_list[[i]] <- SID_sim2_df
}

# Calculate summary statistics for the simulations
summary2_df <- bind_rows(SID_sim2_list) %>%
  group_by(Time) %>%
  summarize(
    S_lower = quantile(S, 0.25),
    S_upper = quantile(S, 0.75),
    I_lower = quantile(I, 0.25),
    I_upper = quantile(I, 0.75),
    D_lower = quantile(Dea, 0.25),
    D_upper = quantile(Dea, 0.75),
    S_mean = mean(S),
    I_mean = mean(I),
    D_mean = mean(Dea)
  )

# Plot the combined graph
combined2_plot <- ggplot() +
  labs(title = "Combined HIV Model Simulation", x = "Time", y = "Population", color = "Group") +  
  scale_color_manual(values = c("Susceptible" = "red", "Infected" = "green", "Death" = "blue")) 

for (i in 1:100) {
  sim2_df_subset <- SID_sim2_list[[i]]
  
}

combined2_plot <- combined2_plot +
  geom_ribbon(data = summary2_df, aes(x = Time, ymin = S_lower, ymax = S_upper, fill = "Susceptible"), alpha = 0.6, color = NA) +
  geom_ribbon(data = summary2_df, aes(x = Time, ymin = I_lower, ymax = I_upper, fill = "Infected"), alpha = 0.6, color = NA) +
  geom_ribbon(data = summary2_df, aes(x = Time, ymin = D_lower, ymax = D_upper, fill = "Death"), alpha = 0.6, color = NA) +
  geom_line(data = summary2_df, aes(x = Time, y = S_mean, color = "Susceptible"), linetype = "dashed", size = 0.4) +  
  geom_line(data = summary2_df, aes(x = Time, y = I_mean, color = "Infected"), linetype = "dashed", size = 0.4) + 
  geom_line(data = summary2_df, aes(x = Time, y = D_mean, color = "Death"), linetype = "dashed", size = 0.4) +
  labs(color = "Group")  # 添加图例标签

print(combined2_plot)

3 Data Comparison

The real data we are comparing comes from the annual number of new HIV cases caused by MSM in Chongqing from 2011 to 2018 (Yu et. al 2017).

data_df <- data.frame(
  Year = 1:8,
  NewInfected = c(660, 849, 1068, 1298, 1397, 1552, 1687, 1782) 
)
plot_real_data <- ggplot(data = data_df, aes(x = Year, y = NewInfected)) +
  geom_line(color = "blue") +
  labs(title = "Real Data - Newly Infectious Over Time", x = "Year", y = "New Infected") +
  theme_minimal()

plot_deterministic <- ggplot(data = HIV_summary, aes(x = Year, y = NewInfected)) +
  geom_line(color = "red") +
  labs(title = "Deterministic Model - Newly Infectious Over Time", x = "Year", y = "New Infected") +
  theme_minimal()

plot_stochastic1 <- ggplot(data = HIV_stochastic_df_summary, aes(x = Year, y = NewInfected)) +
  geom_line(color = "green") +
  labs(title = "Stochastic 1 Model - Newly Infectious Over Time", x = "Year", y = "New Infected") +
  theme_minimal()

plot_stochastic2 <- ggplot(data = HIV_stochastic2_df_summary, aes(x = Year, y = NewInfected)) +
  geom_line(color = "yellow") +
  labs(title = "Stochastic 2 Model - Newly Infectious Over Time", x = "Year", y = "New Infected") +
  theme_minimal()

combined3_plot <- grid.arrange(plot_real_data, plot_deterministic, plot_stochastic1, plot_stochastic2, nrow = 2, ncol = 2)

print(combined3_plot)
## TableGrob (2 x 2) "arrange": 4 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]

4 Discussion and Limitations

For our three simulation results, they fail to accurately reflect the actual disease progression. In real-world data, new HIV infections have been consistently increasing annually, ranging from approximately 660 cases of HIV transmitted among MSM to 1,782 new cases in 2018. In contrast, the three simulated graphs show an initial peak in new infections in the first year, followed by a declining trend each subsequent year, resulting in only a few new cases, numbering in the single digits, by 2018. Several potential reasons may explain this disparity:

Distorted Real Data: HIV is a disease that often goes undetected during its latency period, making it challenging for health authorities to fully capture the actual number of infections in a given year. Additionally, we failed to find the actual annual HIV transmissions specifically attributable to MSM; instead, we calculated this based on the proportion of total new HIV cases transmitted via same-sex encounters each year.

Issues with Simulation: There are several issues with our simulation model:

  1. Limited Variation in ks and pc_t Data: The parameters ks (average number of sexual encounters per person per time step) and pc_t (probability of infection per anal intercourse with an HIV-positive individual) are based on limited data, derived from a single study on HIV in Chongqing and Stanford data respectively. Moreover, the likelihood of HIV transmission varies depending on factors such as active/passive anal intercourse and the stage of infection. However, our simulation fails to fully consider these variances.
  2. Selection of Mortality Rate: HIV is a chronic disease with varying mortality rates at different stages of infection. Therefore, using a single mortality rate for all individuals in the infectious compartment at each time step is not accurate and may prematurely move individuals from the infected compartment to the deceased compartment. However, the simulation does not accurately model the changing mortality rates at different stages.
  3. Homogeneous Population Assumption: Although our population size (N) represents all male homosexuals in Chongqing to the best of our ability, it fails to account for social and demographic differences that may result in unequal participation in unprotected anal intercourse and, consequently, varying levels of infection risk. These issues collectively affect the dynamics of the infected population (I) and the determination of lambda_t, resulting in a simulation model that inadequately captures the fluctuations in HIV infectivity observed in reality.

5 References:

Academic Researches & Journals 1. Wilton, J. (2017). Risk of exposure to HIV/AIDS. Retrieved from https://stanfordhealthcare.org/medical-conditions/sexual-and-reproductive-health/hiv-aids/causes/risk-of-exposure.html#:~:text=Therefore%2C%20unprotected%20sex%20with%20an,exposures)%20for%20receptive%20anal%20sex.

  1. Li, M., Tang, W., Bu, K., Mahapatra, T., Zhang, X., Feng, Y., … Wang, N. (2016). Mortality among people living with HIV and AIDS in China: Implications for enhancing linkage. Retrieved from https://www.nature.com/articles/srep28005

  2. Hongyue, Y., Wei, Z., Rongrong, L., Guohui, W., Kun, S., Yao, Z., Ying, X., Shu, Y., & Yafei, L. (2017, October 30). Epidemiological characteristics and long-term trends of HIV infection and AIDS in Chongqing from 1995 to 2015. Journal of Third Military Medical University. http://aammt.tmmu.edu.cn/Upload/rhtml/201705008.htm

  3. Lu, R., Zhang, X., Zhou, C., Zhang, W., Ouyang, L., Xing, H., … Wu, G. (2020). Trends of HIV, Syphilis and Hepatitis C infections among Men Who have Sex with Men (MSM) in Chongqing, China. Sexually Transmitted Diseases, 1. https://doi.org/10.1097/olq.0000000000001173

R Datasets 1. Bache S, Wickham H (2022). magrittr: A Forward-Pipe Operator for R. https://magrittr.tidyverse.org, https://github.com/tidyverse/magrittr.

  1. Wickham, H., François, R., Henry, L., & Müller, K. (2023). dplyr: A Grammar of Data Manipulation. R Foundation for Statistical Computing. https://cran.r-project.org/package=dplyr

  2. Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag. https://cran.r-project.org/package=ggplot2

  3. Auguie, B. (2017). gridExtra: Miscellaneous Functions for “Grid” Graphics. R Foundation for Statistical Computing. https://cran.r-project.org/package=gridExtra