{'fraction_non_alphanumeric': 0.023919109054920944, 'fraction_numerical': 0.01189993561157084, 'mean_word_length': 5.236466389054135, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Persistent HIV-1 replication maintains the tissue reservoir during therapy

Lymphoid tissue is a key reservoir established by HIV-1 during acute infection. It is a site of viral production, storage of viral particles in immune complexes, and viral persistence. Whilst combinations of antiretroviral drugs usually suppress viral replication and reduce viral RNA to undetectable levels in blood, it is unclear whether treatment fully suppresses viral replication in lymphoid tissue reservoirs. Here we show that virus evolution and trafficking between tissue compartments continues in patients with undetectable levels of virus in their bloodstream. A spatial dynamic model of persistent viral replication and spread explains why the development of drug resistance is not a foregone conclusion under conditions where drug concentrations are insufficient to completely block virus replication. These data provide fresh insights into the evolutionary and infection dynamics of the virus population within the host, revealing that HIV-1 can continue to replicate and refill the viral reservoir despite potent antiretroviral therapy.

genetic divergence over time, nor the emergence of drug resistance mutations have been convincingly documented 9 , 10 . As earlier studies only examined viral sequences derived from the blood of patients who continued to suppress viral replication in that anatomic compartment 11 , the conclusions are not necessarily generalizable to other compartments in the body, particularly to lymphoid tissue where the frequency of infection per cell is mostly higher 12 and the intracellular drug concentrations are much lower than in blood 13 . Under low drug concentrations, the virus may continue to replicate and evolve in sanctuary sites within the reservoir of cells in lymphoid tissue, and remain undetectable in the bloodstream for a time depending on viral population migration dynamics between the two compartments. Here we use a multi-pronged strategy of deep-sequencing, time-calibrated phylogenetic analysis, and mathematical modeling to characterize the distinct temporal structure and divergence of compartmentally sampled viral sequences. We discover ongoing replication in lymphoid tissue sanctuaries of patients despite undetectable blood levels of virus. Our sampling approach differs fundamentally from those of previous studies 14 -16 , which do not address evolutionary dynamics within lymphoid tissue, and better suits investigation of the dynamic nature of the viral reservoir during treatment with potent antiretroviral drugs.

HIV-1 sequence determination
To investigate the evolution and spatial dispersion of virus with high accuracy, we deep sequenced (using the Roche 454 Sciences' GS-FLX sequencing platform) HIV-1 DNA in cells from blood and inguinal lymph nodes collected from three subjects at three separate times (at day 0, and after 3 and 6 months of treatment) described elsewhere 13 . Previous work established that viral sequences contemporaneously sampled from lymphoid tissue in different locations are genetically homogeneous 17 , consistent with CD4+ T cell homing and trafficking 18 . Consequently, detailed assessment of a portion of a solitary lymph node is no more susceptible to bias than wider anatomical sampling. We also sequenced viral RNA in the plasma (day 0) from these three study subjects. Two subjects (1727 and 1679) had wellsuppressed infections (< 48 copies/mL of plasma); and the third subject (1774) continued to have measureable amounts of viral RNA in plasma after 3, but not 6, months of treatment (Extended Data Fig. 1). Subjects 1727 and 1679 were each infected with HIV-1 for approximately 3 to 4 months and were antiretroviral drug naïve before the study. Subject 1774 was infected with HIV-1 for approximately 17 years and was antiretroviral-therapy experienced, but had not received any treatment for at least 1 year prior to the study.
We aligned individual reads with an average length of 548bp to a consensus sequence using reference-guided assembly, and corrected sequencing errors for potentially inflated estimates of genetic diversity 19 . We then used a previously described approach 20 to reconstruct the minimum number of viral haplotypes needed to adequately explain the observed reads. We calculated the sequencing error rate and set the cut-off for the subsequent analyses using a known internal control sequence. We found no significant evidence for recombinant sequences that could bias the analysis. High coverage enabled us to correct PCR and sequencing errors (leaving an average 25,000 final long-reads per sample) to detect variants present in at least 0.04% of the virus population (see Methods). The sequences from the Pol region of HIV-1 that spanned the genomic region encoding the viral enzymes protease or reverse transcriptase showed no evidence of new mutations that confer resistance to the particular antiretroviral drug used (data not shown).
To avoid uncertainties in the haplotype detection due to sparse sampling and prevent systematically biased evolutionary analyses, we established a higher number of template molecules than the depth of the sequence data by a median 9.6 fold (range 3.2 to 362 fold). The high coverage of ultra deep sequencing ensured reliable detection of low-frequency viral variants (see Methods). In support of this conclusion, we found limited variability across inferred haplotypes in two completely independent technical replicates made from the same RNA or DNA sample at each time point from each subject using the same procedures (average Spearman rank correlations coefficient between single nucleotide polymorphism frequencies across replicates of ρ = 0.832, interquartile range across samples, ρ = 0.820 to 0.851; 93.7% of haplotypes above 1% frequency appeared above 1% frequency in the replicate; 97.7% of haplotypes below 1% frequency appeared below 1% frequency in the replicate), indicating that the haplotype representation is not notably biased from random amplification of some sequences and not others. The high degree of concordance between technical replicates validated our approach for the computational characterization of the viral populations, which is robust with respect to experimental error and stochastic effect 20 .

Phylogenies show temporal structure
The inferred haplotypes corresponding to the Gag or Pol regions of HIV-1 were subjected to maximum-likelihood methods of phylogeny estimation (Extended Data Figs. 2 and 3). We masked out the guanosines in the APOBEC3 trinucleotide contexts of the edited sites from the alignments to avoid distortion and retain the phylogenetic information 20 (Extended Data Figs. 2 and 3). Phylogenetic relationships between the distinct haplotypes showed a temporal structure consistent with the molecular clock (continuing nucleotide substitutions occurring at a constant rate), as evidenced by strong correlation between root-to-tip distance and sampling date in the regression analyses (despite the short branches). Branch support computed using an approximate likelihood ratio test and the proportion of sites that are different (p-distance) verify the divergence of haplotypes between day 0 and after 6 months of treatment in most of the Gag and Pol regions of HIV-1 analyzed (Extended Data Table 1). Consistent with ongoing replication rather than sampling of different virus populations in lymph nodes, viral sequences contemporaneously sampled from lymphoid tissue and blood showed a similar degree of divergence. With removal of the haplotypes found to harbor repetitive inactivating base substitutions of guanosine-to-adenosine (G-to-A) the evolutionary lineage emerged that lead up to the APOBEC3-mediated hypermutation event, and the evolutionary rate estimates (range, 6.24 × 10 −4 to 1.02 × 10 −3 substitutions per site per month; Extended Data Table 2) are consistent with those of intra-host virus estimations (range, 5.22 × 10 −4 to 8.42 × 10 −4 substitutions per site per month) 21 . We therefore conclude that continued virus replication contributes to the viral reservoir.

Genetic differentiation due to migration
The pairwise fixation index (F ST ), a standard measure of genetic differentiation between populations, confirmed significant genetic variation between lymph node and blood at each time point (Extended Data Table 3) 22 , 23 . Because F ST measures of population structure   among sampled haplotypes in spatially distinct compartments can be affected by selection  and migration, we used an unrestricted branch-site random effects model to test whether the  proportion of sites along the branches of the phylogeny significantly differ among lineages,  indicating episodic diversifying selection  24 (see Methods and Extended Data Table 4). We restricted the test to internal branches, which capture at least one and likely multiple rounds of virus replication, to mitigate the biasing effects of neutral or deleterious mutations on the ratio of nonsynonymous-to-synonymous substitution rates (ω) estimates where selection has not yet fully filtered such population level variation 25 , 26 . Except for one study subject (1679), where a small proportion of sites (0.3%) were under strong diversifying positive selection along internal tree branches (likelihood ratio test, P~10 −6 ), we found scarce evidence to suggest the virus is evolving in response to strong selective forces. We concluded that the F ST values are more likely due to migration of haplotypes from one compartment to another.

The phyloanatomic history of HIV-1
To infer evolutionary patterns and population dynamic processes from the time-structured sequence data, we used a Bayesian statistical framework that estimates the substitution rate, divergence time, and demographic history of the sampled viral lineages to structure-rooted, time-resolved phylogenies 27 , 28 (see Methods). This approach resolves evolutionary patterns to infer the timing and direction of the key migrations of the virus within hosts. Figure 1 shows the phyloanatomic history of HIV-1 within the study subjects inferred from the ancestral and descendent haplotypes on a Bayesian maximum clade credibility (MCC) phylogenetic tree. The branching patterns in the trees, which reconstruct the origins and trace the flow of HIV-1 within hosts, and the tissue of origin of the internal nodes, show strong statistical support (highlighted by their posterior probabilities). Branch lengths in these time-structured trees (left panels with all haplotypes; right panels with the putative Gto-A hypermutant sequences removed to avoid distortion) represent posterior median estimates of calendar time. The temporal structure in the trees shows a strong clock-like signal and rates consistent with HIV-1 within-host evolution. The best-fit model included a strict molecular clock and assumed a constant population size, though a model with a relaxed molecular clock gave qualitatively similar results (data not shown).
Based on inferred MCC tree topologies and the compartment assignments to unobserved internal nodes, an underlying conformity and strong correlation exist between genetic and anatomic locations and the direction of the virus's spread in the body. A particular pattern recurs: the haplotypes in lymph nodes are the source of viral lineages that migrate from lymph node to blood (Fig. 1). We deduce that viral lineages in blood are derived from replicating virus in lymph nodes with little or no evidence of an additional source in blood.
The time-scaled trees show a strong and significant result and are robust to different substitution models (see Methods). These data, only revealed through temporarily and spatially resolved sampling, further support the conclusion that the pattern does not result from distinct populations of haplotypes being sampled from different compartments, but rather migration and colonization of haplotypes between lymphoid tissue and blood. A structured coalescent model 29 , less prone to potential bias in spatial inference estimates, shows higher migration rates from lymph node to blood (Extended Data Table 5), confirming that the direction of flow is not due to oversampling of a particular anatomic location that would have increased estimates of traffic into that location 30 .
Our results, which reconstruct the dynamics of HIV-1 spread within the body, imply that in patients with no detectable viral RNA in plasma, the virus reservoir is constantly replenished by low-level virus replication in lymphoid tissue. Distinguishing between low amounts of viral replication and pools of latently infected cells that may persist and rekindle HIV-1 infection is methodologically difficult. A small number of HIV-1 sequences isolated at consecutive time points that persisted without evidence of genetic change might be due to long-lived central memory cells, a fraction of which may have reverted to a resting state, or latently infected transitional memory cells that persist by clonal expansions (driven by homeostatic proliferation) or survivorship for long-lived infected CD4+ T cells that contain replication-competent virus 8 , 15 , 22 , 31 -34 . Two of the subjects (1774 and 1679) showed that some haplotypes persist as a single tree branch through time ( Fig. 1d and f), consistent with proliferation of HIV-1 infected cells or long-term cell survival 33 , 34

A spatial dynamic model
Even though viral genetic diversity accumulated over time, the infected cells did not produce new virus with drug-resistance mutations, conferring a putative fitness advantage that might have led to a systemic viral rebound. To provide a mechanistic explanation of this scenario, we developed a drug concentration-dependent mathematical model of virus replication and spread between spatially distinct compartments to explore the deterministic components of viral dynamics 36 . In this model (see Fig. 2 and Supplementary Information), the virus can occupy two spatially distinct compartments that have limited traffic of virus particles or cells. The main compartment has a larger volume and high effective drug concentration. The other, smaller volume compartment (<0.01% of the size of the main compartment) has a lower drug concentration and represents a sanctuary site within the lymphoid tissue reservoir 37 . The model includes competition between two viral strains, one that is sensitive to a potent multidrug regimen and one that is partially resistant to the full complement of antiretroviral drugs, but less fit than the drug-sensitive strain in the absence of treatment.
Within each compartment the balance between the strains is determined by the difference in their replicative fitness without treatment, the effectiveness of drug therapy to curb new rounds of infection in that region (determined by drug concentration) and the susceptibility of each strain to antiretroviral therapy. Figure 3 illustrates a hypothetical fitness landscape, which portrays the relationship between the fitness of each strain and the evolutionary adaptation across a range of drug concentrations. The fitness landscape illuminates the persistence of the pool of virus that is unaffected by antiretroviral drugs. The effective reproductive number R (the average number of secondary de novo infections of cells produced by one infected cell) is a function of the basic reproductive number R 0 and the effectiveness of treatment (R = R 0 in the absence of treatment). The fitness cost for drug resistance determines the difference between the effective reproductive numbers for drug-sensitive and drug-resistant strains (R S and R R , respectively) at zero effective drug concentration. In a competitive system where both strains are present, the maximum of the two effective reproductive numbers R S and R R equals the effective reproductive number R for the system as a whole at that effective drug concentration.
With two spatial compartments, heterogeneity in the distribution of the drug can lead to heterogeneity in R. In the sanctuary site, where drug-penetrance is low, the drug-selective pressure on the replicating virus population is too low to compensate for the fitness costs associated with resistance; thus, the effective reproductive number for a drug-sensitive strain is greater than that of the partially drug-resistant strain (R S > R R ), enabling de novo infection dominated by the drug-sensitive strain 38 . As drug effectiveness increases, the partially drug-resistant strain gradually becomes more fit relative to the drug-sensitive strain such that at some intermediate drug concentrations, characterized by threshold levels, the effective reproductive number for the resistant strain can be greater than that for the drugsensitive strain (R R > R S ). This would allow for de novo infection dominated by the partially drug-resistant strain. In our model, however, we assume that drug concentration in the main compartment exceeds this threshold. At this high drug concentration, infection is no longer sustainable because the effective reproductive numbers are less than unity for both strains (R R and R S ).
By assuming a relatively simple two-compartment model with drug concentration differences between them, the model predicts that virus, dominated by the drug-sensitive strain, can proliferate in a small sanctuary site where the drug concentration is low (Fig. 4a). Although all single point mutations will be generated sufficiently often to prompt partially drug-resistant strains within individuals, their numbers will remain low through competition. The likelihood is exceedingly low (see Supplementary Information and Supplementary Table  2) of stepwise accumulation of mutations from a partially drug-resistant strain to one that confers resistance to all drugs, each of which would have to come to fixation in the absence of drug selection (Extended Data Fig. 4), or the presence or absence of recombination. The model calculations show that increasing drug effectiveness or penetration across spatial regions can affect evolutionary dynamics, and lead either to the emergence of drug-resistant strains (Fig. 4b) or the elimination of ongoing replication (Fig. 4c). Our model predictions fit the data well (Extended Data Fig. 5) and confirm that both competition between strains and regional spatial heterogeneity in antiretroviral concentrations help capture the observed dynamics of the viral reservoir in these well-suppressed patients.
Though probabilistic models suggest (and this model allows) production of partially drugresistant strains, any which arise cannot populate the sanctuary site because of low drug penetration and competition from drug-sensitive strains. Equally, they cannot repopulate the larger, main compartment where drug concentrations preclude any ongoing replication. In agreement with our phylogenetic inferences, these results suggest that the low-level viral replication in lymphoid tissues where antiretroviral drugs concentration are low could allow drug-sensitive strains to grow and spillover to the blood 39 .

Conclusion
This study reveals how dynamic and spatial processes work together to permit HIV-1 to persist within the infected host and avoid development of resistance despite antiretroviral therapy. From these temporally and compartmentally structured sequence data, we conclude that continued virus production from infected cells in lymphoid tissue sanctuary sites, where drug concentrations are not fully suppressive, can continue to refill the viral reservoir and traffic to blood or lymphoid tissue 18 . We further show that the virus does not inexorably develop resistance to antiretroviral drugs because the lower concentrations of drug in the sanctuaries are not sufficient to confer a competitive advantage upon drug-resistant strains.
Our findings explain the failure of treatment intensification to fully suppress de novo infection and highlight issues surrounding the barriers to delivering antiretroviral drugs at clinically effective concentrations in the infectious viral reservoir. The state-of-the-art sequencing approach, innovative time-calibrated phyloanatomic tree construction, and a novel model of compartmentalized intra-host population dynamics provide a new perspective on HIV-1's seemingly untouchable strongholds in the body. Achieving optimal cellular pharmacokinetics and spatial distribution of antiretroviral drugs in lymphoid tissue to fully suppress viral replication and preserve immune function is, therefore, a prerequisite to the elimination of the viral reservoir and ultimately a cure for HIV-1 infection.

Study subjects
The three study subjects were enrolled into a clinical protocol where treatment was started after we obtained peripheral blood by venipucture, inguinal lymph nodes by excisional biopsy, and ileum and rectal biopsies through colonoscopy 13 . Two subjects were antiretroviral drug naïve (1727 and 1679) and one subject (1774) had not received any drugs for at least 1 year before enrollment. Subjects 1679 and 1774 received emtricitabine (FTC), tenofovir (TFV; as the disoproxil fumarate [TDF]), and atazanavir with ritonavir (ATV/R). Subject 1727 received FTC, TDF and efavirenz (EFV). In all three subjects, conventional typing methods confirmed that the plasma virus was sensitive to the potent antiretroviral regimen that they received. Subjects 1727 and 1679 were well-suppressed patients. Subject 1774 continued to have measureable amounts of HIV-1 RNA in plasma after 3, but not 6 months of treatment. We obtained peripheral blood, inguinal lymph node, and terminal ileum and rectum biopsies once again at 3 and 6 months after starting treatment. Laboratory procedures for tissue management, in situ hybridization, and analytical pharmacology are described elsewhere 13 .

Extraction and quantification of viral nucleic acids
Nucleic acids were extracted from frozen cells obtained from blood or lymphoid tissue using the MasterPure Complete DNA and RNA Purification Kit (Epicentre, Madison, WI). Viral RNA was isolated from plasma using the PureLink Viral RNA/DNA Mini Kit (Life Technologies, Carlsbad, CA). HIV-1 was quantified using a quantitative reverse transcription PCR assay. The relative amount of HIV-1 target DNA was normalized to the quantification cycle for a concentration calibrator by using an external standard curve of serial 10-fold dilutions of reference DNA for the Gag region of HIV-1 derived from the plasmid pNL-43.
All reactions were performed in triplicate on the ABI 7900HT sequence detector (Applied Biosystems, Foster City, CA).

Library preparation for deep sequencing
As the number of viral templates in the sampled material is low, we used an amplicon-based deep sequencing strategy. To minimize biased amplification of the target sequence, primer locations in the Gag and Pol regions of HIV-1 were selected on the basis of the alignment positional entropy in the multiple sequences aligned from the Los Alamos National Laboratory HIV-1 sequence database (http://www.hiv.lanl.gov/). The primers were computationally screened for cross-dimer interactions and the concentration of each primer was optimized for amplification. For each sample, blanks were included to screen for contamination. We used designated, physically separated areas within the laboratory to setup PCR, which avoids contact with potentially contaminating amplicons.
To generate the long read-length PCR amplicons sequenced in this study (range, 509bp to 587bp read-lengths per run depending on the gene region being analyzed), we amplified the  for concentration and size distribution. The concentration of the product DNA was normalized before pooling to achieve sequence uniformity across amplicons. Controls were spiked into the reaction to monitor the library construction process and potential index crosscontamination. The known internal control sequence (clonal sequence of 456 bp) introduced into the reaction was used to calculate the single nucleotide error rate and set the cut-off for the sequence analysis where no control errors could be detected. PCR errors are more common in later cycles of amplification, and being limited to small copy numbers they have little impact on the haplotype distribution. Bidirectional sequencing using the 454 Life Sciences' GS-FLX sequencing platform (Roche, Basel, Switzerland) provided independent confirmation of sequence information. The long read-length sequences were sorted based on index sequences, trimmed to remove residual adapter bases from the ends of the reads, and filtered for length and duplicates prior to alignment. Read depth and coverage estimation that met predetermined coverage thresholds were performed as in 20 . Sequencing quality metrics were calculated for all samples using FastQC and only high-quality sequencing libraries were used in the ensuing analyses.

Sequence clean up and assembly
Experimental precision along with deep coverage allows for accurate estimation of the underpinning diversity of the virus population. We binned the sequence reads by multiplex identifier barcodes, and identified and excluded sequencing errors and misaligned regions from the analysis by computational methods. A small number of reads had a disproportionate number of errors that accounted for most of the inaccuracy in the full dataset. After quality filtering and trimming to a uniform length, we proceeded to build the different haplotypes present in each of the samples. We began by collapsing the identical sequences into haplotypes using reference-guided assembly to avoid the use of uninformative sequence repeats. Haplotypes at prevalence above the error threshold (defined by the internal control spiked in the sequencing runs), which corresponded to variants present above 0.04% of the total existing variants in the collapsed alignments, were used in the analysis.
For the processing of the temporally and spatially linked deep-sequencing data for studying viral diversity, the viral sequences were first aligned against the HXB2 reference sequence (GenBank accession number K03455) using Segminator II (version 0.1.1). We then generated a consensus viral sequence for each patient as a reference for assembly to improve the alignment quality. A statistical model that utilizes platform error rates in conjunction with patterns within nucleotide frequencies derived from data obtained from related samples, that is, different compartments within the host or temporally linked samples, was used to separate low frequency platform error from true variation. This model is termed "probabilistic read error detection across temporally obtained reads" (PREDATOR) and is implemented within Segminator II. This statistical framework maintains the reading frame and corrects for deep sequencing errors 19 . We corroborated the number of haplotypes and the frequency of haplotypes that explain the data using a second reconstruction algorithm based on combinations of multinomial distributions to analyze the k-mer frequency spectrum of the sequencing data implemented with QuRe 40 . We screened the sequence alignments for recombinant sequences using the GARD algorithm implemented in HyPhy 41 , 42 .

Sampling variance
The high coverage of massively parallel sequencing is necessary to ensure reliable detection of low-frequency viral variants. A simple calculation based on the geometric distribution shows that in order to guarantee that a viral template occurring at frequency f is detected with probability p or better, it is necessary to sequence at least log (1-p) / log (1-f) -1 templates. A variant of frequency 0.01 (that is, 1%), for example, would require 450 sequences to ensure its detection at probability 0.99 (that is, 99%) or better. Conversely, a study using 100 single-genome sequences would detect a variant of frequency 0.01 with probability 0.64. A low number of input DNA templates derived from one compartment that catches the spillover from another does not account for the complexities of partial observation and spatial heterogeneity that could lead to measurement error elsewhere 16 . This complication emphasizes the challenge in trying to extrapolate from single template sequencing the magnitude and character of the virus population that comprises the viral reservoir.

Maximum-likelihood tree construction
Maximum-likelihood phylogenies were created with PhyML using the general timereversible model with the proportion of invariant sites and gamma distribution of among-site rate variation (GTR+I+Γ 4 ) nucleotide substitution model 43 applying an approximate likelihood ratio test for branch support 44 . We estimated trees on viral sequence sets from which gaps in the alignment were removed and considered as missing data for the reason that maximum-likelihood tree error may increase with inclusion of unreliable sites. We assessed the temporal structure of the trees by performing linear regression on the root-to-tip distances of samples versus the time of sampling and tested the validity of the timedependency of the evolution rate estimates with the assumption of a strict molecular clock using the program Path-O-Gen v1.4 (http://tree.bio.ed.ac.uk/software/pathogen/). We used the Highlighter sequence visualization tool (www.HIV.lanl.gov) to trace commonality between sequences in an alignment based on individual nucleotide changes.

Compartmentalization
We tested for subdivision of viral sequences into sub-populations in the different compartments at each time point. We calculated genetic distances using the Wright's F ST and S nn test statistics 45 , 46 . We used a bootstrap test to determine the confidence of the estimates and performed a permutation test (1000 iterations) to assess the significance levels of the obtained scores.

Identifying selection
We used a modification of a random effects branch-site model to detect positive selection and test whether the phylogeny diverged over time 47 . A likelihood ratio hypothesis test compared the fit of the model using a 3-bin ω distribution (ω 1 and ω 2 in [0,1], ω 3 unrestricted) to describe the evolution of all branches in the tree, to the fit of the model where all ω 3 is restricted to be in [0,1]. We tested whether or not a proportion of sites along internal branches of the intra-host viral phylogeny have been subject to episodic selection (ω > 1) 24 , restricting the test to internal branches to lessen the biasing effects of neutral or deleterious mutations on ω estimates 48 , 49 and serve as a proxy for population level selection 25 , 26 .

Time-calibrated phylogenetic tree construction
To resolve the phyloanatomy, we reconstructed the temporal and spatial dynamics of the viral haplotype lineages with a Bayesian statistical framework using Markov chain Monte Carlo (MCMC) sampling for evolutionary hypothesis testing, as implemented in BEAST version 2.1.2 27 . This approach was used to sample phylogenies from their joint posterior distribution, in which the viral haplotypes are restricted by their known date of sampling, using a simple substitution model described by Hasegawa-Kishino-Yano (HKY) to avoid over-parameterization 50 . Models differing in assumptions on mutation rate and effective population size were run for 100 million generations each and compared using the Bayes factor as implemented in Tracer version 1.6. We determined that the best-fit model included a strict molecular clock and assumed a constant population size. We used a symmetric transition model with constant rates over time that considered a discretized diffusion process among the different compartments. This was formalized as a continuous time Markov chain model to reconstruct the spatial dynamics between compartments. All chains were run for sufficient length and convergence of the relevant parameters was assessed using Tracer version 1.6, ignoring 10% of the chain as burn-in.
We summarized the connections between virus evolution and anatomical compartment history using an annotated MCC phylogenetic tree estimated with BEAST. The model and its parameters were chosen after computing the posterior probability of several models to obtain the discriminatory Bayes factors. Because population structure, whether due to spatial segregation or limitations to gene flow, may affect evolutionary dynamics, we confirmed that the direction of flow was not due to oversampling of a particular environment, by running a two-deme Bayesian inference under a structured coalescent model with a HKY substitution model assuming a strict molecular clock 29 , which is less susceptible to sampling issues than our trait-based analysis. For completeness, we conducted a search for topologies and divergence times assuming a relaxed molecular clock as well. In the analyses performed, HIV-1 showed a high degree of clock-like evolution and a mean nucleotide substitution rate expected to be within the bounds necessary to obtain meaningful phyloanatomic information from sequence data. Using the location of each of the haplotypes, a discrete trait was included in the inference. We used BEAST to estimate the probabilities of each of the possible states.

Statistical analysis
Standard descriptive statistics were performed with the use of the STATA, GraphPad or R packages.
Maximum-likelihood trees were constructed using gene sequences from the Gag region of HIV-1 from lymph node and blood before and after the guanosines within all possible APOBEC3 trinucleotide sequence context of edited sites were masked in the alignments, regardless of their presence in hypermutant or non-hypermutant sequences, to avoid their distortion in the phylogenetic reconstructions. Branch tips are colored according to compartment sampled: red for plasma; gold for lymph node; and blue for blood. The progressive shading of the colors of the branch tips indicate the points in time sampled. Phylogenetic trees reconstructed from the haplotypes in which the guanosines in the APOBEC3 trinucleotide context of the edited sites are masked in the alignments correct the skewing effect caused by clustering of shared haplotypes that harbor repetitive G-to-A substitutions and longer branch lengths caused by a larger number of these mutations in the hypermutated sequences whilst retaining the phylogenetic information. The horizontal scale indicates the expected number of substitutions per nucleotide site per unit time with haplotypes from later time points having diverged more. The Highlighter plots show the haplotypes from the lymphoid tissue and blood time point clusters aligned to the plasma virus sequence from day 0. The particular nucleotide changes are color coded in the alignment (thymidine, red; adenosine, green; cytosine, blue; and guanosine, orange). Magenta circles represent APOBEC3-induced G-to-A change in a trinucleotide context of the edited sites, which are distinguishable from the more random error-prone viral reverse transcriptase and RNA polymerase II replicating enzyme induced mutations 6 . Gene sequences from the Gag region of HIV-1 from Subject 1774, who continued to have measureable amounts of HIV-1 RNA in plasma on treatment, and Subjects 1727 and 1679 who were well-suppressed on treatment (a, b and c, respectively) before and after the guanosines within the particular APOBEC3 trinucleotide sequence context of edited sites were masked in the entire sequence alignment (left and right panels, respectively).
sampled. The horizontal scale indicates the expected number of substitutions per nucleotide site per unit time with haplotypes from later time points having diverged more. The Highlighter plots show the haplotypes from the lymphoid tissue and blood time point clusters aligned to the plasma virus sequence from day 0. The particular nucleotide changes are color coded in the alignment (thymidine, red; adenosine, green; cytosine, blue; and guanosine, orange). Magenta circles represent APOBEC3-induced G-to-A change in a trinucleotide context of the edited sites. Gene sequences from the Pol region of HIV-1 that spanned the genomic region encoding the viral enzyme reverse transcriptase from Subjects  1774, 1727, and 1679 (a, b and c, respectively) before and after the guanosines within the particular APOBEC3 trinucleotide sequence context of edited sites were masked in the entire sequence alignment (left and right panels, respectively).
Extended Data Figure 4. Alternative drug-dependent fitness landscape plots a, Fitness landscape plot for a partially drug-resistant strain that confers relatively low-level resistance to drugs as compared with the fitness costs imposed by the drug-resistant mutations. The drug-resistant strain (blue line) does not outcompete the drug-sensitive strain (orange line) at any effective treatment concentration where it can grow. There are two phases to the dynamics. At lower effective drug concentrations (left of grey line) the drugsensitive strain thrives. Beyond this threshold, neither strain can continuously replicate. b, Fitness landscape plot for a highly drug-resistant strain. This strain confers a high-level of drug resistance relative to the replicative fitness cost imposed by the resistance mutations. At low effective drug concentrations (left of grey line), the drug-sensitive strain outcompetes the drug-resistant strain. At high effective drug concentrations, the drug-resistant strain outcompetes the drug-sensitive strain and can continuously replicate. We argue that, typically, fully drug resistant mutants of this sort neither exist in the viral population of patients before treatment, nor arise through random mutation during the course of antiretroviral therapy (see Supplementary Information and Supplementary Table 2). Drugresistant strains, which are capable of ongoing replication at high effective drug concentrations are not typically generated in individuals because: they are generated in a single step very rarely; and stepwise generation from partially resistant strains is also rare because partially resistant strains are outcompeted in the sanctuary site which constantly refills the pool. The strain specific effective reproductive numbers for the drug-sensitive R S (orange) and drug-resistant R R (blue) strains are shown. For simplicity, only the impact of changes to the effectiveness of a single drug in a single compartment is shown. Circles demonstrate average data from the 3 patients discussed in detail in this study and an additional 9 patients presented elsewhere 13 . Where the average value was indeterminate because of test sensitivity, the data are fitted below the upper limit of the average log 10 infectious units. b, Demonstrates the dynamics over a longer period. The model predicts the persistent low-level viral RNA in plasma. The diamond symbol represents data relating to the long-term persistent virus, as measured using quantitative reverse transcription PCR (see Methods). The optimal model fit parameters are presented in Supplementary Table 1.
Extended Data Table 1 The genetic distance measured between the haplotypes in the Gag or Pol regions of HIV-1 We used sequence data from the Gag or Pol regions (protease [pol1] and retrotranscriptase [pol2]) regions of HIV-1 (range, 315bp to 487bp) to perform a genetic distance analysis. We masked the guanosines in APOBEC3 trinucleotide contexts of the edited sites from the alignments. We compared the proportion of substitutions per site (p-distance) for haplotypes within lymph node or blood from each subject at 6 months to the most common haplotype present at day 0.

Author Manuscript
Author Manuscript Author Manuscript

Author Manuscript
Extended Data Table 4 Episodic selection estimated across sites along internal branches of the phylogeny A branch-site evolutionary model provided a statistical framework to search for evidence of episodic selection on internal branches in the tree. The ω distribution is given for both internal and external branches. Migration rates (in fraction of emigrants per month) were estimated in the Gag region of HIV-1 for the three study subjects using Bayesian inference under the structured coalescence model, assuming a constant population size and a strict molecular clock. Migration rates from lymph node to blood and vice versa are shown with their standard error of the mean (SEM) after running at least 50 million MCMC steps and reaching high values of estimated sample sizes (ESS) for all parameters.   In the main compartment (i=0; the majority of lymphoid tissue and the blood) drug concentration is high (grey and red). In the sanctuary site (i=1; a small fraction of the lymphoid tissue and localized extracellular fluid) drug concentration is low (pink). There are uninfected cells (X i ), long-lived infected cells (Y i ), and short-lived infected cells (Q i ), as well as virus particles (V i ) that can be bound by few (F i ) or many (G i ) receptors on the follicular dendritic cell network. The dashed lines represent the effect of treatment in blocking infection and production of infectious virus particles. For graphical simplicity, we do not show the emergence of drug resistance, the production of noninfectious virus particles, virus clearance, nor cell death.  In the main compartment, antiretroviral therapy is wholly effective against drug-sensitive virus (z 0 = z0 = 1). In the sanctuary site, a, low treatment effectiveness (z 1 = z1 = 0.3) favors drug-sensitive strains. Partially drug-resistant strains will exist, but at levels below those favoring stepwise evolution towards a fully drug-resistant strain. b, Intermediate treatment effectiveness (z 1 = z1 = 0.6), favors partially drug-resistant virus at levels that may suffice for evolution towards a fully drug-resistant strain. c, high treatment effectiveness (z 1 = z1 = 1) favors the decline of all strains and the cessation of virus replication.