Transmission network parameters estimated from HIV sequences for a nationwide epidemic.

BACKGROUND: Many studies of sexual behavior have shown that individuals vary greatly in their number of sexual partners over time, but it has proved difficult to obtain parameter estimates relating to the dynamics of human immunodeficiency virus (HIV) transmission except in small-scale contact tracing studies. Recent developments in molecular phylodynamics have provided new routes to obtain these parameter estimates, and current clinical practice provides suitable data for entire infected populations.
METHODS: A phylodynamic analysis was performed on partial pol gene sequences obtained for routine clinical care from 14,560 individuals, representing approximately 60% of the HIV-positive men who have sex with men (MSM) under care in the United Kingdom.
RESULTS: Among individuals linked to others in the data set, 29% are linked to only 1 individual, 41% are linked to 2-10 individuals, and 29% are linked to ≥10 individuals. The right-skewed degree distribution can be approximated by a power law, but the data are best fitted by a Waring distribution for all time depths. For time depths of 5-7 years, the distribution parameter ρ lies within the range that indicates infinite variance.
CONCLUSIONS: The transmission network among UK MSM is characterized by preferential association such that a randomly distributed intervention would not be expected to stop the epidemic.

Human immunodeficiency virus (HIV) infection, like all sexually transmitted infections, spreads along the sexual contact network, and contact tracing played a key role in establishing the etiology of AIDS [1]. However, more recent studies have struggled to achieve any reconciliation between contact networks revealed by interview data and transmission networks reconstructed from the viral phylogeny [2, 3]. There are several possible reasons for this but 2 important factors are the long period of infectiousness (in the absence of therapy) and the low average risk of transmission per sexual contact [4]. The consequence of the first factor is that although sexual contact networks constructed from interview data usually work within a time frame of 1 year [5, 6], the HIV infection of the subject may have occurred at any time over a period of several years. The consequence of the second factor is reduction of the confidence that any particular contact was associated with infection. Other possible reasons for the lack of agreement that may be associated with particular risk groups include network complexity and anonymous sex. In addition, earlier studies of HIV-infected communities using sequence data were generally based on a small sample of the infected population, reducing the probability of including individuals from the same transmission network [7]. More recently, often working within defined risk groups, a number of studies have identified transmission networks by means of sequence analysis [8-10]. However, this static picture alone does not permit the quantitative description of the network through which HIV has been transmitted [11].

We have introduced an approach to estimating the HIV transmission network structure of a population based on the analysis of HIV sequences generated in the course of routine clinical care to test for antiretroviral resistance [12, 13]. Under UK guidelines from 2003, such tests are now performed for all patients entering therapy, which is recommended for patients with a CD4 cell count of ≤350 cells/μL [14]. Central collection of the assay results from the majority of clinics performing the assays in the United Kingdom allowed a data set of sequences from >26 000 patients to be constructed, including approximately two-thirds of the total number of men who have sex with men (MSM) receiving care for HIV infection in the United Kingdom [15]. We have already demonstrated that phylogenetic analysis of large population-based data sets reveals much more of the transmission pattern than has been seen before [12] and that in combination with relaxed clock models it has shed new light on the dynamics of the HIV infection epidemic in both the MSM and heterosexual contact risk groups.

Here we extend our earlier analysis [12] to the MSM population of the entire United Kingdom. Using viral sequences from almost 60% of MSM diagnosed with HIV infection in 2009, we infer the degree distribution of the transmission networks, based on the time-resolved viral phylogenies, for a range of time depths. We have already shown that that the degree distribution inferred in this way can be approximated by a power law [13], as shown earlier for sexual contact networks inferred from interview data [5]. We here test a series of distributions for fit to the data to determine whether a preferential attachment model is, as has been suggested, the most appropriate model for transmission in this group.

METHODS

Data

HIV sequence data from the protease and partial reverse transcriptase coding regions were obtained from the UK HIV Drug Resistance Database (UKHIVRDB; http://www.hivrdb.org; 2007 download). Before the data set was released for analysis, the records were fully anonymized and delinked from clinical identifiers. In most cases, up to 1200 nucleotides from the start of the protease coding region are available, according to the method used for genotyping. For this study, sequences <850 nucleotides in length were removed, as were sequences with a large number of ambiguities. A check for identical or near-identical sequences (>99.9% identity) was applied to the entire data set. When these were associated with different patient identifications, the sequence groups affected were examined in detail and individual decisions were made whether to discard them on the basis of the available background information. The first available sequence, which was usually obtained before antiretroviral treatment was initiated, was selected from patients with multiple sequences in the database, giving a total of 25 136 sequences each from different patients. In a current substudy with access to sequential viral load data, approximately 12% of the first sequences reported as being obtained before antiretroviral therapy (ART) was initiated appear from analysis of the viral loads to have actually been taken after ART initiation. HIV-1 subtyping was then applied by means of a combination of the Rega [16] method and nucleotide distance comparisons within the data set, which identified 14 560 sequences belonging to subtype B (Figure 1). Because submission of the entire sequence data set to public databases would permit transmission networks to be identified and thus risk breaching patient confidentiality, following Kouyos et al [10] we have submitted a random sample of 10% to GenBank under accession numbers JN100661–JN101948.

Figure 1.

Flowchart of data manipulation and processing. NJ, neighbor-joining; UK HIVRDB, United Kingdom Human Immunodeficiency Virus Drug Resistance Database.

Ethical approval for the use of anonymized data in this work was given by the London Multicentre Research Ethics Committee (MREC/01/2/10; 5 April 2001). Data held in the UKHIVRDB can be accessed for collaborative projects if a project proposal is approved by the steering committee. The proposal form can be downloaded from the database Web site at http://www.hivrdb.org.

Phylogenetic Analysis

The size of the data set required a hierarchical approach to define clusters. For the initial characterization, a genetic distance threshold was calculated following stripping of antiretroviral resistance-associated sites from the sequence data as previously described [13] to avoid clustering artifacts among those sequences that might be affected by ART. The genetic distances were estimated for all pairwise comparisons by means of the HKY+ γ model [17]), and bootstrap testing with PAUP* software [18] was used. PAUP* (version 4.0) output was analyzed using a combination of Python and R scripts. Thresholds of 3.5%, 4.0%, 4.5%, 5%, and 5.5% were examined using a range of bootstrap percentage support values (Supplementary Table 1); 4.5% genetic distance plus 95% bootstrap support was adopted to identify sequences for further analysis as previously described [12]. A total of 5809 sequences with phylogenetic associations that could be resolved using the neighbor-joining method [19] were analyzed using BEAST (version 1.6) software as previously described [13], but with the use of random collections of clusters in 30 pools of approximately 200 sequences per pool for computational feasibility (Figure 1).

Network Creation

A network representing possible HIV transmissions was inferred from the BEAST consensus trees from all the pooled sequence data. The nodes of the network are the HIV sequences, but they represent individuals because only 1 sequence for each individual was used. Two nodes are linked by an edge if the highest posterior density time to the most recent common ancestor (MRCA) of the respective viral sequences in the BEAST consensus tree is less than or equal to some network depth cutoff value. For example, 2 individuals are linked if their viral sequences are predicted to have diverged (ie, the transmission event) <5 years before the most recent sampled sequence. The tree data were processed using a custom R [20] script utilizing the mrca function in the R package ape [21].

Networks of varying time depths (1–16 years) were created by concatenating the subnetworks from the 30 individual consensus trees from the pooled data sets. A network with a time depth of 31 years, representing the maximum depth and maximum number of connections between the sequences, was also created for comparison. Clusters were identified using the clusters function within the R package igraph.

Network Shape Analysis

The degree of a node (number of links per node) approximates the number of possible infected contacts an individual has had within the period of the network. We used the statnet package [22] within R to fit several models to the degree distribution of the nodes within our networks, in addition to fitting a simple power law (not shown). Maximum likelihood parameter estimates, log likelihood values, Akaike information criterion for small sample sizes (AICC) values, and Bayesian information criterion (not shown) values were obtained for the discrete Pareto, negative binomial, Yule, and Waring distributions. For network depths of 3–7 years we also performed 1000 bootstrap replicates to obtain the 95% confidence intervals on the parameters and AICC scores.

RESULTS

Cluster Size Distribution

The UKHIVRDB holds HIV sequence data from the protease and partial reverse transcriptase coding regions obtained in the course of routine clinical care, submitted by collaborating clinical virology laboratories across the entire United Kingdom. We analyzed the 2007 download containing sequences from 25 136 infected individuals, of which 14 560 sequences were classified as subtype B. Direct information on risk activity associated with transmission was not available for this data set, but in an earlier data set from the UKHIVRDB from 2004 for which 1288 individuals infected with subtype B virus had the associated risk group data, the proportions were recorded as follows: MSM, 82% (unpublished data); heterosexual contact, 10%; and injection drug use, 6%. In the 2010 data set of the UKHIVRDB, of 17 334 individuals infected with subtype B, 83% are recorded as MSM (S. Hue, private communication, 28 April 2011), 12% as heterosexual contact, and 3% as injection drug use. From this we can reasonably assume that a little more than 80% of individuals analyzed in this data set were MSM.

At a genetic distance threshold of 4.5%, 7617 individuals (52%) had a link to 1 or more other individuals (Figure 1). Sequences from 5809 of these individuals could be resolved in bootstrapped neighbor-joining trees. Of these, 1728 (29%) were linked only to 1 other, 2408 (41%) were linked to 2–10 others, and 1673 (29%) were linked to ≥10 others, in 97 clusters. As this indicates, the cluster size distribution (Figure 2) is highly right-skewed, as has been described previously for this risk group [12].

Figure 2.

Cluster size distribution of the human immunodeficiency virus transmission network among men who have sex with men in the United Kingdom. Horizontal scale unit, 2; vertical scale, number of clusters; left axis, cluster size of 2–9 (black); right axis, cluster size of 10–106 (gray).

Transmission Dynamics

For a complete analysis of the dynamics of HIV transmission in this population, we used BEAST to sample the 1524 clusters randomly in 30 pools containing all 5809 individuals and obtain dated phylogenies for all of them. The epidemic dynamics were reconstructed from these dated phylogenies. As discussed elsewhere [12], because each sequence represents a different individual, the estimated internode intervals represent the maximum intervals between transmissions. The median of the internode intervals for the entire population was 22 months, with 16% of intervals being ≤6 months. Equivalent estimates of 14 months and 25%, respectively, were obtained earlier for large clusters (≥10 individuals) identified in a single large London clinic [12]. Restricting the comparison to large clusters from the current study gives a median of 17 months and 20% of intervals being ≤6 months. These figures confirm the very important role for early infection in onward transmission in this risk group. The full distribution of internode intervals is given in Supplementary Figure 1.

Network Structure

The network structure was reconstructed using the time to the most recent common ancestors estimated in BEAST for all individuals at a given time depth. Time depths of 1–16 years were examined and compared with the maximal depth of 31 years to test whether cluster size was affected by time depth. Even within the maximum depth network, 688 separate clusters could be identified, indicating that the process of creating networks from pooled data did not artificially create large clusters. The node degree statistics for networks of time depths 1–16 and 31 are displayed in Supplementary Table 2; Figure 3 shows the structure of a large 42-node cluster as revealed at time depths of 3, 5, and 7 years.

Figure 3.

Estimation of network structure. The reconstruction of the largest network cluster is shown at 3 time depths: 3 years, 5 years and 7 years. Networks were initially reconstructed from the dated trees generated by BEAST at all time depths of 1–16 years. The networks for time depths of 3–7 years were used in the distribution fitting (Table 1; Figure 2).

Table 1.

Parameters of the Waring and Negative Binomial Distributions Estimated by Maximum Likelihood From the HIV Transmission Network for Men Who Have Sex With Men in the United Kingdom

	Distribution, parameter
	Waring		Negative binomial
Year	ρ	α	k	μ
3	3.8	0.4	0.33	0.7
4	3.2	0.5	0.33	1.3
5	2.8	0.7	0.32	1.8
6	2.7	1.0	0.34	2.4
7	2.7	1.2	0.35	2.9

Model Fitting to Degree Distribution

In our previous analysis of the UK non-B subtype HIV infection epidemic, we tested the degree distribution observed for fit to a power law [13]. A power law is an approximation to a scale-free network, which can arise from a preferential attachment process, whereby individuals are more likely to link to individuals who already have multiple links. There are explicit models that describe the distribution arising from a preferential attachment process, and in order to apply a more rigorous test, we have compared the HIV transmission network of the UK MSM population to negative binomial, discrete Pareto, Yule, and Waring distributions. These distributions reflect different processes; the simplest is the negative binomial distribution, which is generated if each individual acquires partners at a constant rate, drawn from a γ distribution. The discrete Pareto distribution is a power law distribution, whereas the Yule and Waring distributions arise from specific preferential attachment models [6]. Both the Yule and Waring distributions incorporate (1) the probability that a new link is made to a previously inactive individual and (2) the probability that a new link is made to a person with k partners. The Waring distribution additionally allows for a third class of random, nonpreferential attachments [23, 24]. We used maximum likelihood to estimate parameter values and AICC to compare goodness of fit and found that for all time depths the Waring distribution was the preferred distribution (Figure 4). The greatest difference, at all time depths, was with the discrete Pareto distribution, which clearly fitted less well than any other distribution. The Yule and negative binomial distributions provided the next best fits to the Waring distribution, with only slight differences in AICC for time depths of 3–4 years. An increasing difference for 5 and 6 years in the negative binomial distribution made the Yule distribution our second choice, but at a time depth of 7 years the negative binomial distribution was again preferred of these 2 distributions. The values of the parameter estimates for the Waring and negative binomial distributions are given in Table 1.

Figure 4.

Comparison of the fit of different distributions to the observed network structure. The fit of 4 different distributions, including 2 based on explicit preferential attachment models, to the observed data is shown as the difference between the Akaike information criterion for small sample sizes (AICC) value for each distribution at each time depth and the lowest AICC value at that time depth for time depths 3–7 years (see Table 1).

The difference in fit between the Waring and negative binomial distributions was examined by simulation, using parameter values from the observed data. The bootstrap confidence intervals from these simulations (Figure 5) indicate that the difference in fit was in many cases not significant. Although the negative binomial distribution never gave a lower AICC value than that given by a Waring distribution to a simulated Waring distribution (Figure 5), approximately 10% of bootstrap samples from a negative binomial distribution were fitted better by a Waring distribution. This indicates the level of uncertainty that remains over the existence of a preferential attachment effect in this population.

Figure 5.

Test of the difference in fit between the Waring and negative binomial distributions. Populations were simulated using the parameter values estimated from the data under a Waring distribution (A) and a negative binomial distribution (B), and both distributions were tested for goodness of fit using the Akaike information criterion for small sample sizes (AICC). The simulations show that the negative binomial never fits as well to a Waring distribution as the Waring distribution, although the difference is not great. However, the converse does not hold: the Waring distribution can provide a fit that is not perceptively poorer to a negative binomial distribution.

DISCUSSION

The 14 560 HIV subtype B sequences analyzed here represented two-thirds (68%) of the 21 700 MSM in the United Kingdom who accessed health care for HIV infection in 2006 [25]. Among these sequences, just over one-half (52%) had a recent connection with another sequence in the same data set, inferred as a genetic distance of <4.5%. We established the phylogenetic relationships between the sequences of the virus infecting nearly 6000 individuals, from which we estimated the degree distribution of the entire transmission network. With nearly 30% of these sequences found in 97 groups of 10 or more individuals, this distribution was highly right-skewed (Figure 2), as seen in an earlier study based on a single London clinic [12]. Exhaustive comparisons of cluster numbers at different time depths showed that cluster size was not affected by the time depth chosen (Supplementary Table 2), but the retrospective nature of this study imposes a limitation to our conclusions due to censoring of the database at 2009 and the difficulty of locating recently entered patients correctly in the cluster pattern.

The analysis we have performed reveals the overall rate and distribution of HIV transmission within this population group. In the population of nearly 6000 MSM that was analyzed, a median interval between transmissions of at most 22 months was observed, with 16% of transmissions occurring within 6 months of infection. In large clusters (≥10 individuals), 20% of transmissions occurred within 6 months; our previous study based on a population at a central London clinic showed this proportion could reach 25% [12]. Overall, these findings suggest that the higher infectiousness associated with acute infection [4] plays a general role in accelerating transmission in the MSM population as a whole but that this could differ between populations.

There have been at least 2 large studies of sexual partnerships among MSM in the United Kingdom that can be compared with the results presented here, the representative (probability sampling) National Survey of Sexual Attitudes and Lifestyles (NATSAL) [26] and a community-based survey [27, 28]. The distribution of partner number in the past year reported by the NATSAL survey was as follows: 0 partners, 22%; 1 partner, 32%; 2 partners, 8%; 3–9 partners, 23%; and ≥10 partners, 15%. Substantially higher proportions with higher numbers of partners were reported in the community survey [29], with 53% reporting >5 partners in the past year as opposed to 19% for all MSM in the NATSAL survey. The distribution reported here from the HIV transmission network lies between those of the NATSAL and community-based surveys, with 29% linked to only 1 and 41% linked to 2–10 other individuals with a time depth of 5 years, but a much higher proportion with ≥10 links (29%) than that seen in the NATSAL survey. It could be expected that the transmission network would have a more right-skewed distribution than that of the total contact network because it is conditioned on HIV infection. For that reason, it may not be surprising that the distribution does not differ greatly from that reported by the community-based survey because that would be expected to be biased toward inclusion of a higher proportion of the more highly connected individuals.

Right-skewing of networks has been noted before as a prominent feature of degree distributions derived from surveys of sexual contacts and other networks, and in several cases it has been shown that the degree distribution fits a power law [5, 28, 30]. This has led to the concept that a “small world” model might underlie their dynamics. Because there are distributions that represent more precisely the preferential attachment process [23, 31, 32], we wished to test the distributions predicted for these models against the observed transmission network reconstructed for the HIV infection epidemic among MSM in the United Kingdom. We observed that both of the distributions specifically associated with a preferential attachment process, the Yule and Waring distributions [33, 34], could be fitted to the data, but at all time depths the Waring distribution fitted best (Figure 4). However, the data were also reasonably well fitted by a negative binomial distribution [35], which would not imply preferential attachment. After investigating this result in more detail, we found that a Waring distribution with parameters estimated from the data fitted a binomial distribution simulated with those parameters on approximately 10% of occasions. We cannot completely rule out the possibility that partner choice in this population may be random, but it is substantially more likely to be preferential.

A critical property of truly scale-free networks is that for values of the parameter ρ between 2 and 3 they exhibit infinite variance [23]. Under these conditions, an infection will spread regardless of its transmissibility and no random intervention is sufficient to halt the epidemic. The estimates of ρ obtained here for a time depth of 5–7 years are within that range, implying that any intervention must be targeted to high-degree individuals to be effective [11]. The significance of this conclusion is heightened by the increasing awareness of the potential offered by high-efficacy therapeutic interventions for preexposure prophylaxis. One study performed among women has reported a reduction in the infection rate of 65% [36], whereas the multinational Preexposure Prophylaxis Initiative (iPrEx) trial among MSM reported a reduction of 45% [37]. In the latter study, inclusion criteria specified ≥4 episodes of unprotected anal sex in the past 6 months; however, partner number was not specified. The impact of a treatment program targeted in that way on the most connected individuals in the population could be estimated from the approach described here by incorporation of the quantitative description of the transmission network for the whole population derived from viral sequence data into structured epidemiological models [11]. This would then permit enhanced accuracy in estimates of the cost per infection saved over the population as a whole [38].

Supplementary Data

Supplementary materials are available at The Journal of Infectious Diseases online (http://www.oxfordjournals.org/our_journals/jid/).  Supplementary materials consist of data provided by the author that are published to benefit the reader. The posted materials are not copyedited. The contents of all supplementary data are the sole responsibility of the authors. Questions or messages regarding errors should be addressed to the author.