Longitudinal population-level HIV epidemiologic and genomic surveillance highlights growing gender disparity of HIV transmission in Uganda

Longitudinal surveillance

Between September 2003 and May 2018, nine consecutive survey rounds of the Rakai Community Cohort Study (RCCS) were conducted in 36 inland communities in south-central Uganda (Fig. 1, Supplementary Tables S1-S2, and Supplementary Fig. S1). The results presented in this paper derive from data collected through these surveys, including the population census, the RCCS survey participants, the incidence cohort, and the phylogenetic transmission cohort.

RCCS survey methods have been reported previously ^14,18. In brief, for each survey round, the RCCS did a household census, and subsequently invited all individuals that were of age 15-49 years and residents for at least 1 month to participate in the open, longitudinal RCCS survey; and so data collection was not randomized, and data collection was blind relative to previous interactions with individuals or any personal characteristics apart from age and residency status, and any research questions. Eligible individuals first attended group consent procedures, and individual consent was obtained privately by a trained RCCS interviewer. Following consent, participants reported in a private location, typically a tent at the survey hub, on demographics, behavior, health, and health service use. All participants were offered free voluntary counseling and HIV testing as part of the survey. Rapid tests at the time of the survey and confirmatory enzyme immunoassays were performed to determine HIV status. All participants were provided with pre-test and post-test counseling, and referrals of individuals who were HIV-positive for ART. Additionally, all consenting participants, irrespective of HIV status, were offered a venous blood sample for storage/future testing, including viral phylogenetic studies. Supplementary Table S1 summarises the characteristics of the RCCS participants and HIV-positive participants by age and gender. For the purpose of our analyses, we combined data from three pairs of geographically close areas in peri-urban settings into three communities, and 28 of 36 communities were continuously surveyed over all rounds (Supplementary Table S2). All epidemiologic data collected through RCCS are stored in a database running Microsoft SQL server 2019 and Microsoft Access version 2016.

Ethics declarations

The study was independently reviewed and approved by the Ugandan Virus Research Institute, Scientific Research and Ethics Committee, protocol GC/127/13/01/16; the Ugandan National Council of Science and Technology; and the Western Institutional Review Board, protocol 200313317. All study participants provided written informed consent at baseline and follow-up visits using institutional review board approved forms. This project was reviewed in accordance with CDC human research protection procedures and was determined to be research, but CDC investigators did not interact with human subjects or have access to identifiable data or specimens for research purposes. Participants in the RCCS received 10,000UGX (approximately 2.50USD) in compensation for the baseline and follow-up surveys.

Population size estimates

To characterize changes in population demography, individual-level data on the census-eligible individuals that were obtained during each census were aggregated by gender, 1-year age band (between 15 and 49 years) and survey round (Extended Data Fig. 1a-b, bars). The age reported by household heads in the census surveys tended to reflect grouping patterns towards multiples of 5, suggesting that household heads reported ages only approximately. For this reason, we smoothed population sizes across ages independently for every gender and survey round, using locally weighted running line smoother (LOESS) regression methods that fit multiple polynomial regressions in local neighborhoods as implemented in the R package stats version 3.6.2 with span argument set to 0.5 (Extended Data Fig. 1a-b, line). Model fit was assessed visually without a formal test, suggesting that the data met the assumptions of the statistical model.

Participation rates

To characterize participation rates, we calculated the proportion of RCCS participants in the census-eligible population by gender, 1-year age band and survey round (Extended Data Fig. 1c-d, bars). Following consent, participants reported either their birth date or current age themselves, and accompanying documentary evidence was requested. There were no obvious age grouping patterns of multiple of 5 among participants. Overall, participation rates were lower in men than women (63% vs. 75%). Participation rates also increased with age for both men and women, and were very similar across survey rounds. Considering the grouping patterns by age in the population count data, we again smoothed the participation rates across ages independently for every gender and survey round using LOESS regression as specified above for population size estimation (Extended Data Fig. 1c-d, line). Model fit was assessed visually without a formal test, suggesting that the data met the assumptions of the statistical model.

HIV status and prevalence

All RCCS participants were offered free HIV testing. Prior to October 2011, HIV testing was performed through enzyme immunoassays (EIAs) with confirmation via Western Blot and DNA PCR. After October 2011, testing was performed through a combination of three rapid tests with confirmation of positives, weakly positives and discordant results by at least two EIAs and Western Blot or DNA PCR ⁶². Overall, 99.7% participants took up the test offer across survey rounds, and Supplementary Table S1 documents the number of participants with HIV. From these survey data, we estimated HIV prevalence (i.e., probability for a participant to have HIV) with a non-parametric Bayesian model over the age of participants independently for both genders and survey round. Specifically, we used a binomial likelihood on the number of participants with HIV parameterized by the number of participants and HIV prevalence in each 1-year age band. The HIV prevalence parameter was modeled on the logit scale by the sum of a baseline term and a zero-mean Gaussian Process on the age space. The prior on the baseline was set to a zero-mean normal distribution with a standard deviation of 10. The covariance matrix of the Gaussian Process was defined with a squared exponential kernel, using a zero-mean half-normal prior with a standard deviation of 2 on the scale parameter of the squared exponential kernel and a zero-mean half-normal prior with a standard deviation of 11.3 ((49 − 15)/3) on the lengthscale of the squared exponential kernel. The model was fitted with Rstan release 2.21.0 using Stan’s adaptive Hamiltonian Monte Carlo (HMC) sampler ⁶³ with 10,000 iterations, including warm-up 500 iterations. Convergence and mixing were good, with highest R-hat value of 1.0029, and lowest effective sample size of 830). The model represented the data well, with 98.57% of data points inside 95% posterior predictive intervals, indicating that the data met the assumptions of the statistical model. For the mathematical modeling of transmission flows, we assumed that age- and gender-specific HIV prevalence were the same in non-participants in the RCCS communities as in the participants in these communities.

ART use

The RCCS measures ART use through participant reports since survey round 11. Self-reported ART use reflected viral suppression with high specificity and a sensitivity around 70% in the study population (Supplementary Table S9). We took the following pre-processing steps. For survey round 10, we assumed self-reported ART use to have been on the same levels as in round 11. Next, the ART use field was adjusted to “yes” for the participants with HIV who did not report ART use but who had a viral load measurement below 1,000 copies per milliliter (mL) plasma blood. Further, we considered it likely that with increasingly comprehensive care and changing treatment guidelines ^14,64, ART use in individuals with HIV who did not participate increased substantively over time, and this prompted us to consider as proxy of ART use in non-participants the observed ART use in first-time participants with HIV. Overall, first-time participants represented between 15.26% to 39.87% of all participants. Extended Data Fig. 8a-b exemplifies the self-reported ART use data in male participants and male first-time participants, along with the combined estimate of individuals with HIV in the study population who report ART use, summing over participants and non-participants. These estimates were obtained using the same Bayesian non-parametric model as for HIV prevalence. Convergence and mixing were good, with highest R-hat value of 1.0025 and lowest effective sample size of 978 for the participants and 1.0027, 521 respectively for first-time participants.

The model represented the data well, with 99.67% of data points for the participants inside the corresponding 95% posterior predictive intervals, and 99.24% for the first-time participants, indicating that the data met the assumptions of the statistical model. The resulting, estimated ART use rates in infected men and women are shown in Extended Data Fig. 8c.

Viral suppression

Since survey round 15, HIV-1 viral load was measured on stored serum/plasma specimens from infected participants using the Abbott real-time m2000 assay (Abbott Laboratories, IL, USA), which is able to detect a minimum of 40 copies/mL. Viral suppression was defined as a viral load measurement below 1,000 copies/mL plasma blood following recommendations of the World Health Organisation (WHO) ³⁵. To estimate virus suppression levels in the infected non-participants, we considered again as proxy data on infected first-time participants. Overall, viral load measurements were obtained from 19.3% of participants with HIV in survey round 15 and nearly all (>97.71%) participants with HIV since survey round 16^65,66,67. From these data we estimated the proportion of individuals in the study population with HIV who had suppressed virus, summing over participants and non-participants, using the same Bayesian non-parametric model as for HIV prevalence and ART use. Convergence and mixing were good with lowest R-hat value of 1.0016 and lowest effective sample size of 461 for the participants and 1.0052, 844 respectively for the first-time participants. The model represented the data well, with 98.19% of data points inside 95% posterior predictive intervals and 97.99% for the first-time participants, indicating that the data met the assumptions of the statistical model. For the purpose of mathematical modeling of transmission flows, we next considered the earlier survey rounds 10 to 14, for which viral load measurements were not available. On average, 93% of individuals reporting ART use also had suppressed virus (Supplementary Table S9), leading us to estimate the number of individuals with suppressed virus before 2011 from corresponding ART use data. Specifically, we estimated the proportion of the study population with HIV that was virally suppressed by adjusting the estimated ART use data with the sensitivity of being virally suppressed given self-reported ART use and the specificity of being virally suppressed given self-reported no ART use estimated from round 15 where available, and otherwise from round 16 (Supplementary Table S9). Specificity and sensitivity values by 1-year age bands were linearly interpolated between the midpoints of the age brackets in Supplementary Table S9. The resulting, estimated virus suppression levels in men and women with HIV are shown in Extended Data Fig. 8d, illustrating that the gap in virus suppression levels increased over time.

Sexual behavior

RCCS participants reported to interviewers in each round on aspects of sexual behavior, including the number of sexual partners in the past 12 months within the same community, the number of partners outside the community, and in round 15 also demographic characteristics of up to four partners (Supplementary Table S8). To interpret HIV transmission flows in the context of typical sexual contact networks, we focused on the detailed behavior data collected in round 15 and estimated sexual contact intensities between men and women by 1-year age band, defined as the expected number of sexual contacts of one individual of gender g and age a with the population of the opposite gender h and age b in the same community.

Estimates were obtained with the Bayesian rate consistency model, version 1.0.0, using default prior specifications ⁶⁸. We noted along with previous work ^69,70,71,72 that women tended to report considerably fewer contacts than men (Supplementary Table S8), prompting us to include in the linear predictor of contact rates additional age-specific random effects to capture under-reporting behavior in women. Further, community-specific baseline parameters were added to allow for variation in the average level of contact rates in each community, but the age-specific structure of contact rates was assumed to be identical across communities. The resulting model was fitted to all data pertaining to within-community sexual contacts in the last year, including reports of within-community contacts for which information on the partners remained unreported. Contacts reported with partners from outside the same community were excluded, because male-female contacts have to add up to femalemale contacts only in the same population denominator, and hence under-reporting could only be adjusted for when within-community contacts are considered. The model was fitted with CmdstanR version 0.5.1⁷³ using Stan’s adaptive HMC sampler ⁶³ with 4 chains, where each chain runs 2800 iterations, including 300 warm-up iterations. Convergence and mixing were good, with highest R-hat value of 1.003, and lowest effective sample size of 1,745. The model represented the data well, with > 99% of data points inside 95% posterior predictive intervals, indicating that the data met the assumptions of the statistical model. Supplementary Table S8 reports the estimated sexual contact intensities from men and women in survey round 15, and shows that the estimated, under-reporting adjusted sexual contact intensities in women were considerably higher than those directly reported. The table also shows that the estimated number of sexual contacts from men to women equal those from women to men, and the estimated age distribution of sexual contacts is shown in Fig. 2 and Extended Data Fig. 7.

Data and outcomes from the incidence cohort

The RCCS encompasses both a full census of the study communities and a population-based survey in each surveillance round, which enables identification and follow up of unique individuals over time, and thus provides a comprehensive sampling frame to measure HIV incidence. The RCCS incidence cohort comprises of all RCCS study participants who were HIV-negative at their first visit (baseline) and had at least one subsequent follow-up visit (Supplementary Fig. S1). Individuals in the incidence cohort were considered to be at risk of acquiring HIV after their first visit, and stopped accruing risk at the date of HIV acquisition or the date of last visit. Exposure times were estimated from data collected at survey visit times similarly as in ¹⁴. Individuals in the incidence cohort who remained negative until the last survey round contributed their time between the first and last survey visit to their exposure period. Individuals in the incidence cohort who were found to have acquired HIV must have done so between the visit date of the last round in which they were negative and the visit date of the current round, and the infection date was imputed at random between the two dates. This included incident cases who had no missed visit between the last negative and current visit (type 1) or one missed visit (type 2) as in ¹⁴, but also cases who had more than one missed visit (type 3). Unknown dates were imputed at random 50 times, and individual exposure periods and incident cases were then attributed to each survey round, summed over the cohort, and then averaged over imputations. Supplementary Table S3 and Extended Data Fig. 2 illustrate the age- and gender-specific exposure times and incidence events in each survey round. In sensitivity analyses, we considered only those individuals in the incidence cohort who resided in one of the 28 inland communities that were continuously surveyed across all rounds 10 to 18, and found similar incidence dynamics with slightly faster declines in incidence rates in younger men, although this difference was not statistically significant. No statistical methods were used to pre-determine sample sizes but our sample sizes are similar to those reported in previous publications ¹⁴.

Modelling and analysis

The primary statistical objective was to estimate longitudinal age-specific HIV incidence rates by 1-year age bands across (discrete) survey rounds, separately for each gender. We used a log-link mixed-effects Poisson regression model, with individual-level exposure times specified as offset on the log scale, common baseline fixed effect, and further random effects. The random effects comprised a one-dimensional smooth function on the age space, a one-dimensional smooth function on the survey round space, and an interaction term between age and survey round. The functions were specified as one-dimensional Gaussian processes, similar as in the model for estimating HIV prevalence. Alternative specifications, including two-dimensional functions over the participant’s age and survey round, and without interaction terms between age and survey rounds were also tried. We did not consider incidence trends in continuous calendar time because study communities were surveyed in turn, and so the incidence data within each round are structured by communities, which would require further modeling assumptions to account for. Due to the large number of individual observations, models were fitted using maximum-likelihood estimation (MLE) with the R package mgcv version 1.8-38 in the R language ⁷⁴, to each of the 50 data sets with imputed exposure times for each gender. Numerical convergence was examined with the gam.check function. Within and between sample uncertainties in parameter estimates, from the variability of the estimation procedure and the data imputation procedure, were incorporated in the age-, gender- and survey round-specific incidence rate estimates by drawing 1,000 replicate incidence rate estimates from the MLE model parameters and associated standard deviation obtained on each of the 50 imputation data sets, and then calculating median estimates and 95% prediction intervals over the 1, 000 × 50 Monte Carlo estimates (Fig. 1c). Model fits were evaluated by comparing predicted HIV incidence infections estimates to the empirical data. To assess model fit, incident cases were predicted using the Poisson model parameterised by replicate MLE incidence estimates. Overall, model fit was very good, with 98.80% [98.10-99.49] data points inside the 95% prediction intervals across all 50 imputed data sets and the fitted model was consistent with all the available data (Extended Data Fig. 6), indicating that the data met the assumptions of the statistical model. The Akaike information criterion was used to identify the best model for each gender, and the best model was as described above (Supplementary Table S4).

Data from the transmission cohort

Within the RCCS, we also performed population-based HIV deep-sequencing spanning a period of more than 6 years, from January 2010 to April 2018. The primary purpose of viral deep sequencing was to reconstruct transmission networks and identify the population-level sources of infections, thus complementing the data collected through the incidence cohort.

The RCCS viral phylogenetic transmission cohort comprises of all participants with HIV for whom at least one HIV deep sequence sample satisfying minimum quality criteria for deep-sequence phylogenetic analysis is available (Supplementary Fig. S1). For survey rounds 14 to 16 (PANGEA-HIV 1), viral sequencing was performed on plasma samples from participants with HIV who had no viral load measurement and self-reported being ART-naïve at the time of the survey, or who had a viral load measurement above 1,000 copies/mL plasma. We used this criterion because viral deep sequencing was not possible within our protocol on samples with virus less than 1,000 copies/mL plasma, and because self-reported ART use was in this population found to be a proxy of virus suppression with reasonable specificity and sensitivity ^14,22. Plasma samples were shipped to University College London Hospital, London, United Kingdom, for automated RNA sample extraction on QIAsymphony SP workstations with the QIAsymphony DSP Virus/ Pathogen Kit (Cat. No. 937036, 937055; Qiagen, Hilden, Germany), followed by one-step reverse transcription polymerase chain reaction (RT-PCR) ⁷⁵. Amplification was assessed through gel electrophoresis on a fraction of samples, and samples were shipped to the Wellcome Trust Sanger Institute, Hinxton, United Kingdom for HIV deep-sequencing on Illumina MiSeq and HiSeq platforms in the DNA pipelines core facility. Primers are publicly available ⁷⁵. For survey rounds 17 to 18 (PANGEA-HIV 2), viral load measurements were available for all infected participants and viral sequencing was performed on plasma samples of individuals who had not yet been sequenced and who had a viral load measurement above 1,000 copies/mL plasma. Plasma samples were shipped to the Oxford Genomics Centre, Oxford, United Kingdom, for automated RNA sample extraction on QIAsymphony SP workstations with the QIAsymphony DSP Virus/ Pathogen Kit (Cat. No. 937036, 937055; Qiagen, Hilden, Germany), followed by library preparation with the SMARTer Stranded Total RNA-Seq kit v2 - Pico Input Mammalian (Clontech, TaKaRa Bio), size selection on the captured pool to eliminate fragments shorter than 400 nucleotides (nt) with streptavidin-conjugated beads ⁷⁶ to enrich the library with fragments desirable for deep-sequence phylogenetic analysis, PCR amplification of the captured fragments, and purification with Agencourt AMPure XP (Beckman Coulter), as described in the veSEQ-HIV protocol ⁷⁷. Sequencing was performed on the Illumina NovaSeq 6000 platform at the Oxford Genomics Centre, generating 350 to 600 base pair (bp) paired-end reads. Sequencing probes are publicy available ⁷⁸. A subset of samples from survey rounds 14 to 16 with low quality read output under the PANGEA-HIV 1 procedure was re-sequenced with the veSEQ-HIV protocol. To enhance the genetic background used in our analyses, additional samples from the spatially neighboring MRC/UVRI/LSHTM surveillance cohorts and other RCCS communities were also included. For sequencing, the following software were used, QuantStudio Real-Time PCR System v1.3, Agilent TapeStation Software Analysis 4.1.1, Clarity Version 4.2.23.287, FreezerPro 7.4.0-r14598, and LabArchives Electronic Lab Notebook 2023. We restricted our analysis to samples from 2,172 individuals that satisfied minimum criteria on read length and depth for phylogeny reconstruction and sub-sequent inferences. Specifically, deep-sequencing reads were assembled with the shiver sequence assembly software, version 1.5.7⁷⁹. Next, phyloscanner version 1.8.1²⁶ was used to merge paired-end reads, and only merged reads of at least 250 bp in length were retained in order to generate 250bp deep-sequence alignments as established in earlier work ²².

Deep-sequencing was performed from 2010 (survey round 14) onwards, but because sequences provide information on past and present transmission events, we also obtained information on transmission in earlier rounds and calculated sequence coverage in participants that were ever deep-sequenced at minimum quality criteria for phylogenetic analysis. Specifically, we required that individuals had a depth of ≥30 reads over at least 3 non-overlapping 250bp genomic windows. Individuals who did not have sequencing output meeting these criteria were excluded from further analysis, and these were largely individuals sequenced only in PANGEA-HIV 1, and were primarily associated with low viral load samples ^77,80. In total, we deep-sequenced virus from 1,978 participants with HIV of who 559 were also in the incidence cohort. Supplementary Table S5 characterizes HIV deep-sequencing outcomes in more detail. No statistical methods were used to pre-determine sample sizes but our sample sizes are similar to those reported in previous publications ^20,28,51.

Reconstruction of transmission networks and source-recipient pairs

The HIV deep-sequencing pipeline provided sequence fragments that capture viral diversity within individuals, which enables phylogenetic inference into the direction of transmission from sequence data alone ^22,79,81. First, potential transmission networks were identified, and in the second step transmission networks were confirmed and the transmission directions in the networks were characterized as possible. In this study, the first step was modified from previous protocols ²² to ease computational burden, while the second step was as before performed with phyloscanner, using version 1.8.1.

In the first step ⁸², to identify potential transmission networks, HIV consensus sequences were generated as the most common nucleotide in the aligned deep-sequence fragments that were derived for each sample. We then calculated similarity scores between all possible combinations of consensus sequences in consecutive 500 bp genomic windows rather than the entire genome to account for the possibility of recombination events and divergent virus in parts of the genome. Similarity score thresholds to identify putative, genetically close pairs were derived from data of long-term sexual partners enrolled in the RCCS cohort similarly as in ^22,82, and then applied to the population-based sample of all possible combinations of successfully sequenced individuals. Overall, 2525 putative, genetically close individuals were identified, and these formed 305 potential transmission networks.

In the second step, we confirmed the potential transmission networks in phylogenetic deep-sequence analyses. We updated the background sequence alignment used in phyloscanner to a new sequence data set that included 113 representatives of all HIV subtypes and circulating recombinant forms and 200 near full-genome sequences from Kenya, Uganda, and Tanzania, obtained from the Los Alamos National Laboratory HIV Sequence Database (http://www.hiv.lanl.gov/). The deep-sequence alignment options were updated to using MAFFT version 7.475 with iterative refinement ⁸³, and additional iterative re-alignment using consistency scores in case a large proportion of gap-like columns in the first alignment was detected. Deep-sequence phylogeny reconstruction was updated to using IQ-TREE version 2.0.3 with GTR+F+R6 substitution model, resolving the previously documented deep-sequence phylogenetics branch length artefact ^20,84. Confirmatory analyses of the potential transmission networks were updated to using phyloscanner version 1.8.1 with input argument zeroLengthAdjustment set to TRUE. From phyloscanner output, we calculated pairwise linkage scores that summarise how frequently viral phylogenetic subgraphs of two individuals were adjacent and phylogenetically close in the deep-sequence phylogenies corresponding to all 250bp genomic windows that contained viral variants from both individuals ^22,26. Similarly we calculated pair-wise direction scores that summarise how frequently viral phylogenetic subgraphs of one individual were ancestral to the subgraphs of the other individual in the deep-sequence phylogenies corresponding to all 250bp genomic windows that contained viral variants from both individuals and in which subgraphs had either ancestral or descendant relationships ^22,26. Phylogenetically likely source-recipient pairs with linkage scores ≥ 0.5 and direction scores ≥ 0.5 were extracted, and only the most likely source-recipient pair with highest linkage score was retained if multiple likely sources were identified for a particular recipient. The resulting source-recipient pairs were checked further against sero-history data from both individuals where available. If sero-history data indicated the opposite direction of transmission, the estimated likely direction of transmission was set to that indicated by sero-history data.

Infection time estimates

The shape and depth of an individual’s subgraph in deep-sequence phylogenies also provide information on the time since infection, and since the sequence sampling date is known thus also on the infection time ⁸⁵ and the age of both individuals at the time of the infection event. We used the phyloTSI random forest estimation routine with default options, which was trained on HIV serocon-verter data from the RCCS and other cohorts, and uses as input the output of the phyloscanner software ²⁷. Individual-level time since infection estimates were associated with wide uncertainty (Extended Data Fig. 4a), and for this reason we refined estimates for the phylogenetically likely recipient in source-recipient pairs using the inferred transmission direction, age data, and where available longitudinal sero-history data. Specifically, we refined plausible infection ranges as indicated in the schema in Supplementary Fig. S2. Here, the dotted red rectangle illustrates the 2.5% and 97.5% quantiles of the phyloTSI infection time estimates for the phylogenetically likely recipient (x-axis) and transmitting partner (y-axis). We incorporated evidence on the direction of transmission by requiring that the date of infection of the phylogenetically likely recipient is after that of the transmitting partner (filled red triangle). Sero-history and demographic data were incorporated as follows. For both the recipient and the transmitting partner, the upper bound of the infection date was set as the 30^th day prior to the first positive test of the participant ⁸⁶.

The lower bound of the infection date was set to the largest of the following dates, the date of last negative test if available, the 15^th birthday, or the date corresponding to 15 years prior the upper bound ⁸⁷. The refined uncertainty range of the infection time estimates of the phylogenetically likely transmitting partner and recipient are illustrated as the purple triangle in the schema above, and obtained as follows. Firstly, we defined individual-level plausible ranges, by intersecting the range of dates consistent with the phyloTSI predictions and sero-history data. If the intersection was empty, we discarded the phyloTSI estimates. Then we intersected the rectangle given by the cartesian product of the plausible intervals for source and recipient with the half-plane consistent with the direction of transmission. Finally, infection dates were sampled at random from the refined uncertainty range, so that the median infection date estimates correspond to the center of gravity of the triangle (cross). In sensitivity analyses, we further integrated estimates of transmission risk by stage of infection ⁸⁸, though this had limited impact on the estimates (see Sensitivity analyses section below). In cases where the likely transmitting partner in one heterosexual pair was the recipient partner in another heterosexual pair, the above infection date refinement algorithm was applied recursively so that the refined infection date estimates were consistent across pairs. Finally, the transmission events captured by each source-recipient pair were attributed to the survey round into which the posterior median infection time estimate of the recipient fell, and in cases where the median estimate fell after the start time of a round and the end time of the preceding round, the event was attributed to the preceding round.

In total, we identified 539 source-recipient pairs that involved participants from the 36 survey communities and further individuals from the background data set. In 13 of the 539 source-recipient pairs, available dates of last negative tests indicated that only the opposite transmission direction was possible and in these cases the inferred direction of transmission was set to the opposite direction. The resulting pairs included 501 unique recipient partners, and for reach we retained the most likely transmitting partner. To identify pairs capturing transmission events within the RCCS inland communities, we restricted analysis initially to 236 heterosexual source-recipient pairs in whom both individuals were ever resident in the 36 survey communities. Of these, 142 pairs were from men to women and 94 from women to men. Infection times were estimated for all sampled individuals and refined for the recipient partners in the 236 heterosexual source-recipient pairs. For 4 recipient partners, the phyloTSI estimates were ignored as they were incompatible with inferred transmission direction and survey data, and was based on sero-history data only. The phylogenetically most likely location of both individuals at time of transmission was estimated as their location at the RCCS visit date that was closest to the posterior median infection time estimate. Using this location estimate, 233 of the 236 heterosexual source-recipient pairs were estimated to capture transmission events in RCCS inland communities and were retained for further analysis. A further 6 recipient partners had posterior median infection time estimates outside the observation period from September 2003 to May 2018 and were excluded, leaving for analysis 227 heterosexual source-recipient pairs that captured transmission events in RCCS inland communities during the observation period. This excluded 88 potential source-recipient pairs from our study due to ethical considerations and prior analyses suggesting these pairs most likely represent partially sampled transmission chains (i.e., “false positives”) ²².

Statistical framework

We next estimated the sources of the inferred population-level HIV incidence dynamics from the dated, source-recipient pairs in the viral phylogenetic transmission cohort. Overall, inference was done in a Bayesian frame-work using a semi-parametric Poisson flow model similar to Xi, X. et al. ²⁹, that was fitted to observed counts of transmission flows with transmission direction g ⟶ h (male-to-female or female-to-male), time period p (R10-R15 and R16-R18) in which the recipient was likely infected, and 1-year age bands i, j of the source and recipient populations respectively, where The target quantity of the model is the expected number of HIV transmissions in the study population in transmission direction g ⟶ h (male-to-female or female-to-male), survey round r (R10 to R18) in which infection occurred, and 1-year age bands i, j of the source and recipient populations respectively, which we denote by . We considered that the expected number of HIV transmissions in the study population is characterized by transmission risk and modulated by the number of infectious and susceptible individuals, which prompted us to express in the form of a standard discrete-time susceptible-infected (SI) model, where is the transmission rate exerted by one infected, virally un-suppressed individual of gender g and age i on one person in the uninfected (“susceptible”) population of the opposite gender h and age j in a standardized unit of time in round r. With model (2), we express expected transmission flows with a population-level mechanism of how transmission rates from individuals with unsuppressed HIV act on the susceptible population, and we preferred model (2) over a purely phenomenological model of the for the generalizing insights it provides. The main simplifying approximations in (2) are that all quantities on the right-hand side of (2) are in discrete time and constant in each round, meaning we approximate over changes in population size, HIV prevalence, and viral suppression at a temporally finer scale, and assume further that one generation of transmissions occurs from individuals with unsuppressed HIV in each round. Importantly, in this framework, we can then relate the expected transmission flows to the HIV incidence dynamics and the data from the longitudinal incidence cohort by summing in (2) over the sources of infections, where is the incidence rate per census-eligible, susceptible person of gender h and age j in round and per unit time . Estimates of were calculated in units of 100 person-years as described above and shown in Fig. 1c, and we will constrain the semi-parametric Poisson flow model using these estimates. From the model output, we are primarily interested in the transmission flows and transmission sources during each round as quantities out of 100%, defined respectively by In words, (4b) quantifies the sources of infection in individuals of gender h and age j in round r such that the sum of over i equals one, and (4c) quantifies the sources of infection in the entire population in round r that originate from the group of individuals of gender g and age i such that the sum of over g and i equals one. The width of the boxplots in Fig. 2b shows (4b) and Fig. 2a, c show (4c).

Specification of susceptible and infected individuals

The number of the susceptible population of gender h and age j was calculated by multiplying the smoothed estimate of the census-eligible population of gender h and age j (shown in Extended Data Fig. 1a-b) with 1 minus the posterior median estimate of HIV prevalence in census-eligible individuals of gender h and age j of round r (calculated as described further above). To specify the number of individuals with unsuppressed HIV of gender g and age i, we multiplied the smoothed estimate of the census-eligible population of gender g and age i of round r (shown in Extended Data Fig. 1a-b) with the posterior median estimate of HIV prevalence in the census-eligible population of gender g and age with 1 minus the posterior median estimate of the proportion of census-eligible individuals of gender g and age i in round r that have suppressed HIV (calculated as described further above and shown in Extended Data Fig. 8d). The start and end times of each survey round, and were set as shown in Fig. 1b and specified in units of years, so that the transmission intensity is also expressed in units of years.

Bayesian model

We first present the likelihood of the observed counts of transmission flows under the semi-parametric Poisson flow model that is parameterized in terms of (2). The phylogenetically reconstructed source-recipient pairs capture only a subset of incidence events, and so it is important to characterize the sampling frame. As in Xi, X. et al. ²⁹, we consider the unknown transmission events in round r and assume these are sampled at random within each strata with probabilities that factorise into sampling probabilities of sources of age i and gender g and sampling probabilities of recipients of age j and gender . Using (4a), we also let , where Z_r is the total number of infection events in round r.

Because we have data from both the transmission and incidence cohorts, we are able to constrain the sampling problem the detection probabilities of incidence events. Specifically, setting we let and set the detection probability to the proportion of the expected number of incident cases of gender h and age j that could be phylogenetically reconstructed in time period p, for all rounds r in the two time periods R10-R15 and R16-R18. We focused in (5) on time periods due to the limited phylogenetic count data. The advantage in constraining the transmission model with the detection probabilities (5) is that the estimates of the transmission model will be consistent with the incidence dynamics that we already estimated with data from the incidence cohort. Re-arranging terms between Binomial and Multinomial models, we obtain which shows that the sampling probabilities of recipients terms of the detection probability of infection events, weighted by the relative contribution and sampling of source-specific transmission events to the same incidence group. We still need to specify to complete the sampling model. Here, we approximated the sampling probability of sources with the proportion of individuals of age i and gender h with unsuppressed virus in round r that were ever deep-sequenced. Note that the sampling model (6) will alter the posterior mean transmission flows only when the sampling probabilities and differ between age and gender strata in the same round. Extended Data Fig. 5 visualizes our specifications of and shows that the sampling differences between age and gender groups are relatively modest in any given round, which suggests that the adjustments on the inferred transmission flows based on our modelled sampling probabilities will be modest.

In the semi-parametric Poisson flow model of Xi, X. et al. ²⁹, the sampling model (6) can be analytically integrated out based on standard thinning properties, which in turn allows us to express the likelihood of observing the phylogenetic data with where is the posterior median estimate of the log rate of sexual contacts within communities in one year between one person of age i and gender g and one person of age j and gender h that we estimated from the sexual behavior data, and the remaining terms quantify the transmission probability per sexual contact on the log scale. The model is designed in such a way that the log sexual contact rates describe a fixed age-specific non-zero mean surface, and the remaining parameters describe age-specific random deviations around the mean surface. With this approach, any inferred deviations in transmission rates relative to sexual contact rates are informed by the phylogenetic data and robust to prior specifications on the random deviations. Specifically, γ₀ is the baseline parameter characterizing overall transmission risk per sexual contact, γ_g is a gender-specific offset which is set to zero in the female-to-male direction and a real value in male-to-female direction, γ_r a round-specific offset which is set to zero for the first survey round 10, and γ_p is a time period specific off-set which is set to zero for the first time period. We assume the age-specific structure of transmission rates in terms of the transmitting partners (denoted by i) and recipients (denoted by j) are similar across similar ages, and so we can exploit regularising prior densities ²⁹ to learn smooth, latent transmission rate surfaces from the sparse data shown in Extended Data Fig. 3. In detail, we modeled the age-specific structure of transmission rates non-parametrically with 2 time-invariant random functions with two-dimensional inputs on the domain [15, 50] × [15, 50] that characterize age-age interactions in transmission risk for each gender, 2 × 8 random functions with one-dimensional inputs that characterize time trends in the age of recipients for each gender for survey rounds after round 10, and 2 random functions with one-dimensional inputs that characterize time trends in the age of transmitting partners for each gender for the second time period. We attach to each of these random functions computationally efficient B-splines projected Gaussian process (GP) priors ⁸⁹, which we constructed by describing the random functions with cubic B-splines over equidistant knots and modeling the prior relationship of the B-splines parameters with GPs with squared exponential kernels with variance and length-scale hyper-parameters, denoted respectively by s² and 𝓁. The prior densities of our Bayesian model are where the 2 × 8 recipient-specific time-varying 1D B-splines GPs each have squared exponential kernels with hyper-parameters , the 2 source-specific time-varying 1D B-splines GPs each have squared exponential kernels with hyper-parameters , and the 2 time-invariant 2D B-splines GPs each have squared exponential kernels with hyper-parameters and decomposed as follows, We constrain the model further with a pseudo-likelihood term so that the model’s implied incidence rate in (3b) is around the MLE incidence rate estimate obtained from the incidence cohort. We took this approach in lieu of fitting the model to both the source-recipient and individual-level incidence exposure data to bypass extreme computational runtimes ¹², and in the context that the source-recipient data are not informative of incidence dynamics ⁹⁰. Specifically, we fitted log-normal distributions to the 1, 000 × 50 Monte Carlo replicate rate estimates for individuals of gender h and age j in round r (see above) using the lognorm R package version 0.1.6⁹¹, and then set where and denote respectively the parameters of the fitted log-normal distributions, and the left-hand side is calculated from (7b) and matches the model’s incidence rate in (3b).

Computational inference

Model (7-10) was fitted with Rstan version 2.21.0, using Stan’s adaptive HMC sampler ⁶³ with 4 chains for 3,500 iterations including 500 warm-up iterations. Convergence and mixing were good, with highest Rhat value of 1.0027 and lowest effective sample sizes of 1444. The model presented the data well, with 99.63% data point inside 95% posterior predictive intervals and the fitted model was consistent with all the available data (Extended Data Fig. 6), indicating that the data met the assumptions of the statistical model. There were no divergent transitions, suggesting non-pathological posterior topologies.

Closing half the gap in viral suppression rates in men relative to women

In this scenario, we considered the impact of reducing by half the gap in the proportion of men with unsuppressed HIV compared to the same proportion in women. To this end, we first calculated for each 1-year age band the average of the estimated proportion of census-eligible infected men in round 18 with suppressed virus and the same proportion in women, . Next, we set to the smoothed estimate of census-eligible men of age i in round R18 multiplied with the posterior median estimate of HIV prevalence in census-eligible men of age i, and with .

Closing the gap in viral suppression rates in men relative to women

In this scenario, we considered the impact of achieving the same proportions of men with unsuppressed HIV as in women. To this end, we set to the smoothed estimate of census-eligible men of age i in round R18 multiplied with the posterior median estimate of HIV prevalence in census-eligible men of age i, and with .

95-95-95 in men

In this scenario, we considered the impact of achieving viral suppression in 85.7% (0.95 × 0.95 × 0.95) in each 1-year age group of men with HIV. The number of remaining men with unsuppressed HIV in round 18, , was calculated by multiplying the smoothed estimate of the census-eligible men of age i in round R18 with the posterior median estimate of HIV prevalence in the census-eligible men of age i, and with 1 − 0.857.

Sensitivity in incidence rate estimates to the GAM incidence model specification

The longitudinal age-specific HIV incidence rates of the central analysis were estimated with a log-link generalized additive effects Poisson regression model with a linear predictor comprising relatively simple main and interaction effects by age and survey round, fitted to individual-level 0/1 incidence outcomes and exposure times specified as offset on the log scale. To assess sensitivity against the relatively simple linear predictor, we considered a more complex mean specification comprising in-dependent LOESS smoothers to capture age-specific incidence trends in each survey round, and fitted this mean model for computational reasons to crude HIV incidence rates. Specifically, we fitted LOESS regressions as implemented in the R package stats version 3.6.2 with span argument set to 0.7 to the crude age-, gender- and round-specific HIV incidence rates in all 50 imputation data sets, and weighted by the corresponding, group-level aggregated exposure times. The HIV incidence rate estimates under the LOESS model had as expected a smaller mean absolute error against the crude estimates as compared against the GAM model (0.0048 [0.0046-0.0051] versus 0.0053 [0.0051-0.0056]) (Supplementary Fig. S3). Overall, the contribution of men to incidence was more variable across rounds while the shifts in the median age at infection were similar in the central and this sensitivity analysis (Supplementary Table S10).

Sensitivity in incidence rate and transmission flow estimates to limited communities

Over time some communities were added and others left the RCCS (see Supplementary Table S2). We repeated our analysis on the subset of 28 consecutively surveyed communities. We found similar incidence rates with slightly faster declines in male new infections and larger gender disparities (Supplementary Fig. S4). All other primary findings remained insensitive (Supplementary Table S10).

Sensitivity in estimating transmission flows to uncertainty in infection time estimates

In the central analysis, phyloTSI infection time estimates associated to source-recipient pairs were refined using the inferred transmission direction, age, and serohistory data. To assess sensitivity to the infection time estimates used, we inferred transmission flows on the basis of the raw phyloTSI infection time estimates as long as they were compatible with the inferred transmission direction, and otherwise on the basis of the refined estimates. Overall, we found source-recipient pairs were potentially allocated to earlier or later time periods reflecting the wide uncertainty in infection time estimates, though across the sample the age distribution of sources and recipients was remarkably stable (Extended Data Fig. 4). All primary findings were insensitive to using the raw infection time estimates (Supplementary Table S10).

Sensitivity in time since infection estimates to higher transmissibility during acute infection

In the central analysis, transmission flows were estimated using the centre of gravity of the uncertainty region associated with the refined infection time estimates. To account for higher transmission rates during acute infection of the transmitting partner ⁸⁸, we assumed that the transmission hazard was 5 times higher in the first 2 months after infection of the transmitting partner as compared to the following period, and obtained the resulting mean infection time estimate under this assumption by generalizing our Monte Carlo approach used in the central analysis to an importance sampling approach under piecewise linear transmission hazards. The primary results were insensitive to these changes as less than 5% of source-recipient pairs were attributed to different survey rounds (Supplementary Table S10).

Sensitivity in estimating transmission flows to right censoring of likely transmission pairs

The RCCS transmission cohort was defined retrospectively and so it is possible that some transmission events, especially in later rounds, remain as of yet unseen because the corresponding individuals are not yet in the survey or do not yet have virus deep-sequenced. To assess sensitivity to right censoring, we excluded from analysis those source-recipient pairs for which virus of the source or the recipient was deep-sequenced only after rounds 17, 16 and 15. The primary findings were insensitive to these analyses because the probabilities of detecting infection events in the phylogenetic data changed accordingly (Supplementary Table S10 and Supplementary Fig. S5).

Sensitivity in estimated transmission flows to limited sample size of likely transmission pairs

The number of observed infection events in the incidence cohort was ≈4 times larger than the number of reconstructed transmission events, prompting us to explore the effect of sampling uncertainty on the transmission flow estimates. We bootstrap sampled source-recipient pairs at random with replacement three times, and repeated inferences on these bootstrap samples. Our primary findings remained insensitive (Supplementary Table S10).

Sensitivity in estimated transmission flows to modelled sampling estimates

The sampling adjustments in (6) require assumptions including that sampling is independent of infection and transmission, independent between source and recipient, at random within strata, and well approximated by approximating sources with individuals with unsuppressed virus. We repeated flow inferences without any adjustments and without adjustments for potentially unequal sampling of sources. Our primary findings were insensitive across these analyses (Supplementary Table S10).

Sensitivity in transmission flow estimates to the phylo-SI model specification

In the central analysis, the log transmission rates that underpin the estimated transmission flows were estimated using the linear predictor in (7c), and this model specification was associated with overall smallest mean absolute error and posterior predictive coverage as shown in Supplementary Table S6 against the following alternative models, Models specifying transmission rates without a round-specific random function on the age of infected individuals, (12a)-(12c), did not fit the data well (Supplementary Table S6). The remaining models, (12d)-(12f) performed as well as the model used in the central analysis (Supplementary Table S6) and our primary findings remained insensitive (Supplementary Table S10).

Sensitivity in counterfactual intervention impacts to assumptions on viral suppression levels in non-participants

Infection and viremia in the non-participant census-eligible population remained unknown and in the central analysis, we considered as proxy of virus suppression levels among non-participants data from first-time participants. We performed two sensitivity analyses, assuming first that all non-participants with HIV were also viremic across all rounds, and assuming second that virus suppression was identical among non-participants and participants of the same age, gender and survey round. Together, the two scenarios likely encompass the true, unknown viral suppression levels in non-participants. These scenarios were implemented by updating the number of individuals with viremia in (2), and refitting the model. The sensitivity analysis assuming all non-participants with HIV were viremic resulted in larger predicted incidence reductions in women around 75%, while the sensitivity analysis assuming virus suppression was the same among non-participants as among participants of the same age, gender and survey round resulted in similar predicted incidence reductions in women than in the central analysis (Supplementary Table S10).

Sensitivity in counterfactual intervention impacts to potentially higher HIV prevalence in non-participants

In the central analysis, we assumed that HIV prevalence was the same in participants and non-participants of the same age, gender and survey round. We considered three sensitivity analyses, assuming first that prevalence was 25% higher in male non-participants compared to male participants of the same age, gender and survey round, assuming second that prevalence was 25% higher in female non-participants compared to female participants of the same age, gender and survey round, and assuming third that prevalence was 25% higher in female and male non-participants compared to female and male participants of the same age, gender and survey round respectively. These scenarios were implemented by updating the number of virally unsuppressed individuals in (2), and refitting the model. Our primary findings remained insensitive (Supplementary Table S10).

Sensitivity in counterfactual intervention impacts to lower viral suppression thresholds

Different definitions of HIV suppression are currently operational, and we considered the effect of lower thresholds to define viral suppression (<200 copies/mL) than in the central analysis (<1,000 copies/mL). This scenario was implemented by re-estimating the age- and gender-specific proportions of individuals with HIV in the study population who had suppressed virus at the lower threshold, re-calculating gaps in viral suppression levels in men relative to women, and re-calculating the additional number of men needed to reach and maintain viral suppression in the counterfactual intervention scenarios. We found slightly smaller gender gaps in viral suppression at the lower threshold and the predicted incidence reduction in women in the counterfactual that assessed closing the suppression gap in men was around 45%, and all other findings remained insensitive (Supplementary Table S10).