Integrating genealogical and dynamical modelling to infer escape and reversion rates in HIV epitopes
Duncan Palmer, John Frater, Rodney Philips, Angela McLean, Gil McVean
1 Department of Statistics, 2 Nuffield Department of Clinical Medicine, 3 Institute for Emerging Infections, The Oxford Martin School, 4 Department of Zoology, South Parks Road, University of Oxford, Oxford, United Kingdom
The rates of escape and reversion in response to selection pressure arising from the host immune system, notably the cytotoxic T-lymphocyte (CTL) response, are key factors determining the evolution of HIV. Existing methods for estimating these parameters from cross-sectional population data using ordinary differential equations (ODE) ignore information about the genealogy of sampled HIV sequences, which has the potential to cause systematic bias and over-estimate certainty. Here, we describe an integrated approach, validated through extensive simulations, which combines genealogical inference and epidemiological modelling, to estimate rates of CTL escape and reversion in HIV epitopes. We show that there is substantial uncertainty about rates of viral escape and reversion from cross-sectional data, which arises from the inherent stochasticity in the evolutionary process. By application to empirical data, we find that point estimates of rates from a previously published ODE model and the integrated approach presented here are often similar, but can also differ several-fold depending on the structure of the genealogy. The model-based approach we apply provides a framework for the statistical analysis of escape and reversion in population data and highlights the need for longitudinal and denser cross-sectional sampling to enable accurate estimate of these key parameters.
Cytotoxic T-lymphocytes (CTLs) are implicated in the control of human immunodeficiency virus 1 (HIV-1). In fact, they are thought to be the most important mediators in reducing viraemia in individuals able to control HIV infection, showing association with repression of viral replication in long-term non-progressors (LTNPs) [1, 3, 2]. Epitopes are presented to CTLs by human leukocyte antigen (HLA) class I proteins at the surface of almost all nucleated cells in the body. The collection of epitopes which may be presented by the HLA class I molecules is determined by an individual’s combination of alleles at these highly variable loci. Mutations in or close to epitopes in the viral sequence can result in alterations to the binding affinity to the HLA class I, reduce CTL recognition, or abrogate T-cell receptor (TCR) binding. Such mutations are known as escape mutations. Examples of escape mutations have been described in almost all proteins encoded in the HIV-1 genome [4, 5, 6, 7, 8, 9, 10, 11, 12], with the strongest signal of association with host HLA type at the HLA-B locus [15, 14]. After an escape event takes place, escape mutations in the virus may be transmitted between individuals and thus have the potential to spread across the infected population [4, 16], or revert through selection pressure within hosts whose immune responses do not drive escape in a given epitope [13, 17]. Associations between HLA type and putative CTL escapes have been demonstrated statistically in population studies , though these results have been called into question , and it is suggested that the frequency at which escape events take place is lower than previously thought. There is strong interest in understanding the selective pressures applied to the virus at the level of the population as there are clear implications for any putative vaccine. To date, simple ordinary differential equation (ODE) based models have been used to estimate the expected time to escape and reversion by using cross-sectional data across hosts. Such estimates make use of only a small portion of the available data, namely presence or absence of an escape mutation and the HLA type of the sampled hosts (which we denote ), and disregard any remaining sequence information. These methods also make assumptions about the independence of the sampled data which could potentially lead to bias in estimates. Furthermore, deterministic models only provide point estimates and thus cannot provide meaningful confidence regions which account for phylogenetic uncertainty.
Figure 1 illustrates the inference problem that we are addressing. We wish to infer three rates, the rate of viral escape (switching from dark green to dark red in Figure 1), the rate of viral reversion (switching from pink to light green in Figure 1), and the transmission rate. If the underlying transmission tree was known, the problem would be straightforward. However, we only have a collection of tip labellings (sequence data and HLA information) which are the culmination of an embedded subtree of the full process. To make statements about parameters of the full transmission tree, we must reconstruct the subtree together with the dynamic processes occurring along its lineages through time.
We apply dynamic programming  in combination with existing software to combine phylogenetic and statistical approaches with well studied ODE based modelling to integrate available sequence data. By combining these two frameworks we determine more informed estimates and credible regions of population level escape and reversion rates which incorporate the underlying dependency structure present in the viral genealogy. Specifically, we first estimate the underlying genealogy using available HIV sequence data, and then overlay a simple ODE compartment based model of the processes of escape and reversion adapted from  (shown in Figure 1). We may then determine the probability of observing conditional on this genealogy. By integrating over genealogies informed by the viral sequence data, we can then generate credible regions for the parameters of interest. We envisage similar methodologies being applied to a wide range of problems. Our model represents an addition to the highly active area of phylodynamics, in which both stochastic and deterministic approaches are being developed [22, 23, 24].
To test the robustness of our integrated method, we perform a series of simulations and compare the results to those of an existing dynamical model . Notably, our method does not require assumptions about the start time of the epidemic, or transmission rate during the exponential growth phase, as these are estimated by the model. When nucleotide substitution rates are fixed at estimates generated from empirical sequence data, we find that in simulation studies our model successfully estimates escape and reversion rates. By altering the nucleotide substitution rate, we find that a lack of information about the genealogy (through lower substitution rates) can dramatically affect escape and reversion rate estimations using our integrated approach, though we find that the rates of substitution found in HIV are sufficiently large for this effect to be considered negligible.
The integrated approach is then applied to estimate escape and reversion rates in four previously identified epitopes. The four epitopes were chosen in order to explore as much of the space of escape and reversion rates as possible on the basis of previous estimates generated using population level data. Again, we compare our integrated method to the ODE method which generated these estimates . We illustrate the benefit of setting our approach in a model-based framework by demonstrating some simple hypothesis tests.
Our model provides evidence for the hypothesis that rates of escape and reversion within-host are slower than published estimates generated experimentally from individual case studies , and highlights the large amount of uncertainty inherent in estimates that make use of cross-sectional population data.
We wish to estimate escape and reversion rates within host, taking
into account the dependency structure between sampled individuals
arising from the phylogenetic tree. Throughout, we define hosts with
an HLA type known to confer an escape mutation in the virus as HLA
matched, and those without such an HLA type as HLA mismatched. We have
the HLA types and cross-sectional viral sequences from a collection of
hosts, taken from the The Swiss-Spanish intermittent treatment trial
(SSITT) [26, 13]. Epitopes and
associated HLA types are considered one at a time, independent of
other epitopes. We consider four epitopes, selected to explore and
test our method across as wide a range of escape and reversion rate
parameter space as possible, based previous estimates
. The chosen epitopes are shown in table
1 column 1. Throughout, we abbreviate these epitopes
by their first three amino acids (e.g. ). By removing
the epitope under consideration from sequences and determining the
presence or lack of an escape mutation, we consider data from two
processes. The sequence data with epitope removed, , allows us to perform
inference on the genealogy, . The combination of HLA type and
presence or lack of escape, (which we refer to as HLA/escape
information), provides information about the dynamical processes shown
in Figure 1, which occur along the
lineages of over time. By assuming escape information is
uninformative about , we may consider these two processes
separately, with the second conditional on the first. Details are
provided in supplementary text S1. Adding a collection of timestamped
data taken from the Los Alamos HIV sequence database 
to this HLA typed cross-sectional sequence data, we create a DNA
multiple sequence alignment  and perform some
data trimming. The alignment is then passed to the program BEAST
 to obtain samples from the posterior density
of the genealogy and evolutionary parameters. Taking a sample of
genealogies from the BEAST output, we consider the embedded tree for
which we have HLA and escape information, and apply our pruning
algorithm using the concepts in . Our tree
pruning is based on the transitions shown in Figure
1. Defining , where
denotes , and
denotes in the epitope under investigation. In our model,
transmission between individuals takes place at rate ,
individuals are assumed to be HLA matched with probability
. Transitions between the four states may be described by the
matrix in equation 1, where describes the
transition from state to state . , where increases towards the present. is the time
before the present of the most recent common ancestor (MRCA), where
. is the probability that a lineage at time in the
past does not have any sampled descendents
. Details of the derivation, excluded from the original paper are provided in supplementary text S2. The scaling is required to avoid
double counting of transmission events. We assume the state at the root node does not have the escape mutation. In calculating the likelihood we must account for the fact that each internal node in the genealogy corresponds to exactly one transmission event. By varying the escape and reversion rates, , over parameter space and integrating over the sample from the BEAST output we determine a posterior density for , and define credible regions for our estimates. A description of the data and full details of the method are provided in supplementary text S3.
Implicit in our model is a variety of assumptions as in . Firstly, we assume homogeneous mixing in the infected population. We suppose that the infected population is in exponential growth and assume constant rates of escape and reversion across the population with and without the relevant HLA, and ignore variation within individuals’ viral populations. All escape mutations within a given epitope are assumed to occur at the same rate, and HLA types are considered to two digits. It is assumed that a single individual seeded the epidemic, and we suppose individuals with the corresponding HLA type are always able to make an immune response. Finally, recombination is not considered in our estimation of the genealogy .
Under our model, assuming the population of infected individuals is in exponential growth, the process generating states at the leaves is a birth-death process  together with escape and reversion events. The birth rate is and the death rate, , is equal to the rate of becoming non-infectious (through death or otherwise). Throughout our simulations, , (established from the average of a BEAST run on gag data), . Where required, we sample 500 genealogies from BEAST output.
Testing the integrated method and comparison to existing approach: We test and compare our model to a differential equation approach previously described . It can be shown that the dynamics of our model when match the relative proportions through time of the ODE model during exponential growth, and that this is equivalent to assuming a completely star-like genealogy in which all lineages emanate from the MRCA (see supplementary text S4). We choose five parameter sets, and for each, generate a full birth-death tree with an MRCA of 25 years together with HLA and escape information. We then set such that the expected number of present day tips is 200, and sample extinct tips at rate such that the expected number of historically sampled tips is 50. We then simulate sequence information at the tips, , each 500 nucleotides long, using mutation parameters, (which we set as the average of the required parameters from a BEAST run on gag data) and a model of substitution. Using HLA/escape information at the present day tips, and sequence information , we apply both methods four times for each parameter set. The ODE model requires an estimate of the transmission rate, death rate and the initiation time of the epidemic. We fix these parameters at their true values. We bootstrap 10,000 times to provide an estimate of sampling uncertainty under the ODE model. An estimate of the death rate, , and sampled proportion at the present, , is required under the exponential coalescent tree prior, which we set at the truth. In order to define an approximation to confidence regions for sampling under the ODE method, we use the Mahalanobis distance measure  (see methods) widely used in cluster analysis, which takes covariance between escape and reversion rates into consideration. Two instances of estimating the parameter set are shown in Figure 2. Four simulations of the five parameter sets are shown in Figures S1 - S5. We investigate this further by running 1000 simulations for each parameter set on independent sampled birth-death trees and comparing the maximum a posteriori (MAP) estimate under our model on the true tree to the ODE point estimate, shown in Figures 2 and S6. Finally we examine the ability of our method to estimate true underlying parameters over a large number of simulations. Setting the truth at , we apply our full integrated method 100 times to distinct collections of sequence and HLA/escape data generated as before. Results are shown in Figure 2.
Robustness to tree topology and impact of mutation rate: In order to estimate the impact of uncertainty in the tree topology we performed two tests. Firstly we permute tip labellings in the true tree 250 times and re-estimate using our pruning algorithm, setting the true values at . Secondly we simulated data as above, but multiply each substitution rate by a factor of and , before applying our method. The truth is set at .
Simulations with known underlying tree:
By determining the MAP of 1000 instances of HLA and escape data, supposing we know the true underlying genealogy, we find our integrated method best estimates the true rates when the truth lies in the centre of our range of parameter simulations. MAP estimates within a factor of of the truth were for parameter sets . This makes sense, very high or low escape leads to a lack information to discern from either an infinite rate, or a rate of 0. Such estimates will result by chance under non-zero rates, the extreme example is data in which all individuals show escape. Our integrated approach substantially outperforms the ODE method when escape and reversion rates are slow.
Simulations with unknown underlying tree: We conduct four simulations over the five parameter sets shown in Figures 2, S1 - S5, and a further 100 simulations for shown in Figure 2. We find large variation in the size of credible regions across genealogies, particularly for low underlying rates. In the large simulation shown in Figure 2, MAP estimates lie within the and credible regions 79, 83, 86, 88, 92, 95 and 100 times out of 100.
Comparison to ODE method: We find that in general, the ODE method estimates the truth well, particularly when escape and reversion rates are fast. This makes sense: along branches culminating in tips, fast escape and reversion leads to convergence to the equilibrium distribution, which is independent of the tree. For slower rates however, our integrated method is favourable. This can be seen in single simulation runs in Figure 2 and S3 - S7, and in rate estimates assuming the true tree is known. Here, the integrated approach performs far more favourably, with a tighter distribution about the truth in all parameter sets. However, as would be expected, the signal begins to drop under both models as underlying rates are reduced further. For the underlying parameter sets , the proportion of ODE point estimates and MAP estimates within a factor of were and respectively (the corresponding values within a factor of 10 were and respectively). Knowledge of the transmission tree increases accuracy of rate estimations across the parameter space of , particularly when escape and reversion rates are low.
Robustness to tree topology and impact of mutation rate: By shuffling tip labellings we investigate robustness of estimates to the tree topology. We find, in addition to an increase in variance of estimations, a systematic bias towards higher rate estimates, shown in Figure S7. This makes intuitive sense: reduction in knowledge of tip labellings will act to randomise any clustering (or lack of clustering) present in the true tree of escaped and wild-type strains at the epitope under consideration, leading to a forced increase (decrease) in the lower bound of the number of escape and reversion events in the tree, increasing (reducing) rate estimates. To investigate the effect of mutation rate on estimates, we multiply and divide substitution rates by a factor of 5, this is displayed in Figure S8. As the mutation rate is increased, we see a reduction in variance and increase in accuracy of estimates as we would expect. Going forward, it is important that this observation is considered in pathogens in which mutation rates are far lower and phylodynamic methods are beginning to be applied. Figure S7 and S8 demonstrate that a lack of knowledge of the underlying genealogy can seriously impact any parameter estimations leading to potentially spurious results.
Analysis of real data: The result of applying the integrated method to the available SSITT data is displayed in Figure 3 and summarised in Table 1. The cross-sectional proportion sampled, , is set at 0.003, based on incidence data . The first aspect to notice, which was also present in many of our simulations, is the similarity between the simple ODE method and the MAP from the integrated approach. This is not surprising as the purely dynamical model  can be written as a composite likelihood, which results from the assumption that all lineages are independent and of equal weight. The success of composite likelihoods is reflected in the similarity seen. However, looking more closely at the and epitopes (see Table 1), we see that the estimate of the underlying genealogy is playing a strong role. The maximum clade credibility trees for the BEAST runs of and using TreeAnnotator  are shown in Figure 3. We see that in the case of , the escape rate is times lower in the MAP estimate than the ODE estimate. This is reflected in the clustering seen in the tree. Of the 24 individuals who have a consensus escape, 13 occur in clusters of 2 or more. Moreover, singleton escaped lineages coalesce deep into the tree. These combine to reveal the existence of a lower escape rate than that seen in the ODE approximation, combined with transmission of escapes indicative of a very low reversion rate. Contrast this to the rate approximation in , in which there is an excess of escaped singletons found in the maximum clade credibility tree, leading to increased escape and reversion rate when incorporating the genealogy. Secondly, we notice the large amount of uncertainty arising through the evolutionary process in our estimations. This can be explained through two main factors: the lack of seen transmissions close to the present, and uncertainty in the states deep into the ancestry. Improvement in each will result from increased sampling. An increase at the present will increase the number of recent coalescent events, and increased historical sampling will add confidence to inferred states deep in the genealogy. Throughout, we have considered the exponential coalescent as our prior on the genealogy, as a birth-death prior with historical sampling and sampling at the present is currently untested in BEAST. Use of a birth-death prior would increase uncertainty in the time to the most recent common ancestor due to its stochastic nature. Indeed, the deterministic nature of the coalescent underestimates uncertainty in in an epidemic setting. However, it allows us to sample large portions of the infected population completely correctly, increasing our power to infer dynamical parameters. Unfortunately, the long terminal branches indicative of exponential growth means that it will always be difficult to estimate parameters of such dynamical models.
|Epitope||Gene||HLA||Escape||to escape||to reversion|
|TSTLQEQIGW||gag (108-117)||B57, B58||T3N||0.096||0.35||0.35||2.60||2.54|
|KRWIILGLNK||gag (131-140)||B27||R2K, R2G, R2Q||0.073||5.00||6.68||6.50||9.41|
|RPMTYKAAV||nef (77-85)||B7||T4S, Y5F||0.166||69.4||37.8||44.3||26.0|
Major advantages of the model-based framework are that we obtain meaningful credible intervals for our parameter estimates, and gain a statistical framework in which hypotheses about these parameters can be tested. For example, consider the null hypothesis that escape and reversion rates are common across the four epitopes, with the alternative that they each have distinct rates. Using a likelihood-ratio test, we reject the null hypothesis (). Testing the difference in escape and reversion rates between and , we find the data cannot reject the null hypothesis that there is a common escape and reversion rate across these two epitopes (). Validation of our use of the likelihood-ratio test is given in the supplementary text S5.
We have combined a dynamical modelling approach with cross-sectional sequence data to infer escape and reversion rates at CTL epitopes within hosts whilst taking the underlying genealogy into account. Previous models have seeked to achieve two distinct goals. Firstly, the detection of CTL escape mutations in HIV, and secondly, to estimate rates of viral escape and reversion of these implied CTL mutants. Methods to detect association first looked for signals of association between HLA types and putative escape mutations without considering the genealogical structure  of viral sequences. This was followed by partial consideration of the viral genealogy  and more complex network based methods [36, 37]. Working towards the second goal, escape and reversion rates associated to known escape mutations were estimated, again without any inclusion of dependency between samples . We compared our integrated approach to this ODE based method through simulations and parameter estimates using cross-sectional data from the SSITT study cohort. Our model is set in a statistical framework and makes use of the information present in sequence data for parameter estimates. We gain meaningful credible regions which consider uncertainty in the true underlying genealogy, and our integrated method provides a probabilistic framework in which hypotheses can be tested. Our most striking conclusion is the large amount of uncertainty present in rates estimates utilising cross-sectional sequence data. Great care must therefore be taken before strong conclusions are made on the basis of such estimations.
Under the ODE approach outlined in , sequence information outside the epitope under consideration is redundant. Cross-sectional data are considered to have arisen independently, up to initial conditions. We show that this is equivalent to the assumption of a completely star-like phylogeny in supplementary text S4. Given this major assumption of the ODE method, we expected our integrated method to perform more accurately in simulations. We find that this is indeed the case, the presented model consistently outperforms the simpler ODE based approach, particularly when strong clustering of escape mutations is more likely to be present in the true tree. Despite this improvement, we find that the ODE method recaptures escape and reversion rates relatively well in simulations, inspite of its simplifying assumptions.
In order to incorporate the underlying dependency structure present in the genealogy of our samples, and assuming exponential growth, there are two major processes to decide between: the exponential coalescent [39, 40, 38] (which we chose) and the sampled birth-death process [30, 41]. In the case of the birth-death process, a large number of likelihoods describing the process have been constructed [44, 43, 30, 42], differing only in what the author(s) decide to condition on. We summarise the links between each likelihood in supplementary text S6. Our choice between exponential coalescent and sampled birth-death process determines the prior on the genealogy. One major distinction is that the sampled birth-death process considers a subpopulation within a larger stochastically growing population, and the exponential coalescent assumes that the subpopulation is a tiny subset within a very large deterministically growing population. There are benefits and drawbacks to each. Under a birth-death prior, the sum of seen transmissions informed by the prior and unseen transmissions given by the matrix is the overall transmission rate over time - a desirable property which the coalescent lacks . The birth-death process incorporates early stochasticity, which more accurately represents the truth in an epidemic setting. Additionally, no coalescent assumption is required, making the prior suitable for datasets in which the number of samples is comparable to the total population size. This is becoming increasingly relevant as such datasets are becoming more commonplace . Under the exponential coalescent, inclusion of time-stamped datasets is straightforward. In contrast, under the birth-death process, assumptions about the sampling rate of these historical events must be made , which are often invalid for many datasets. However, both are only priors on tree shape and if the data is strong, the distribution from which the trees are sampled will be near identical. We show an interesting link between the two processes in supplementary text S7.
In our estimates we are fundamentally restricted by coalescence events in the genealogy. Long terminal branches are indicative of exponential growth, yet the greatest power to inform our parameter estimates comes from coalescence events occurring in the recent history of the virus. Thus, obtaining extra information from the genealogy is intrinsically difficult. Greater sampling at the present will increase the occurrence of recent coalescence events, and provides greater power to distinguish high and low rates from infinity and zero respectively. However, the inclusion of more and more sequence data calls the coalescent assumption into question. Dense sampling can also lead to a breakdown of the assumed connection between the genealogy and the transmission tree due to lineage sorting . Including time-stamped data with with tip information would allow estimation of ancestral states with greater confidence, and thus increase our ability to infer escape and reversion rates. If longitudinal and cross-sectional sequence data could be combined, this would dramatically increase power to estimate rates. Unfortunately, current methods cannot support this extension due to the use of the genealogy as a proxy for the transmission tree. It is clear that greatest power to estimate these parameters lies in longitudinal sampling within cohorts of hosts, but here we create another collection of issues: recombination plays a far larger role within host, and we would require longitudinal sequences across a large number of individuals in order to make any meaningful statements about rate estimates across the infected population.
Models which attempt to integrate the underlying genealogy are currently being developed to incorporate epidemiological dynamics outside the exponential growth phase  assumed under this model. Another potential improvement would be to co-estimate escape and reversion rates within the MCMC scheme. Our model also makes many assumptions about the underlying biological processes. For example, overlapping epitopes which are prevalent across the HIV genome  mean that mutations conferred by one HLA type could be incorrectly inferred to be the result of selection due to another HLA type. HLA types are considered to two digits and escape mutations within a given epitope are grouped together due to the relatively small dataset. All individuals with the restricting HLA type are assumed to be capable of making a response which drives selection at the epitope under investigation. With larger datasets, such complications could be included in a similar model in the same framework with more parameters.
We have constructed a model which integrates sequence data and considers the evolutionary history, transmission, and set of dynamical processes together. The model was created using existing techniques, and we use it to address pressing practical questions.
Funding was provided by the EPSRC to D.P. Many thanks to Helen Fryer and Jonas Schlüter for providing thoughtful comments on the manuscript.
- 1. Migueles, S. A. et al. 2000 HLA B*5701 is highly associated with restriction of virus replication in a subgroup of HIV-infected long term nonprogressors. Proceedings of the National Academy of Sciences, 97, 2709-2714. (doi:10.1073/pnas.050567397)
- 2. Poropatich, K & Sullivan, D. J. 2011 Human immunodeficiency virus type 1 long-term non-progressors: the viral, genetic and immunological basis for disease non-progression. The Journal of general virology, 92, 247-268. (doi:10.1099/vir.0.027102-0)
- 3. Rinaldo, C., Huang, X. L., Fan, Z. F., Ding, M., Beltz, L., Logar, A., Panicali, D., Mazzara, G., Liebmann, J., & Cottrill, M. 1995 High levels of anti-human immunodeficiency virus type 1 (HIV-1) memory cytotoxic T-lymphocyte activity and low viral load are associated with lack of disease in HIV-1-infected long-term nonprogressors. Journal of Virology, 69, 5838-5842.
- 4. Allen, T. M. et al. 2004 Selection, transmission, and reversion of an antigen-processing cytotoxic T-lymphocyte escape mutation in human immunodeficiency virus type 1 infection. Journal of virology, 78, 7069-7078. (doi:10.1128/JVI.78.13.7069-7078.2004)
- 5. Borrow, P. et al. 1997 Antiviral pressure exerted by HIV-l-specific cytotoxic T lymphocytes (CTLs) during primary infection demonstrated by rapid selection of CTL escape virus. Nature Medicine, 3, 205-211. (doi:10.1038/nm0297-205)
- 6. Draenert, R. et al. 2004 Immune selection for altered antigen processing leads to cytotoxic T lymphocyte escape in chronic HIV-1 infection. The Journal of experimental medicine, 199, 905-915. (doi:10.1084/jem.20031982)
- 7. Kelleher, A. D. et al. 2001 Clustered mutations in HIV-1 Gag are consistently required for escape from HLA-B27 restricted cytotoxic T lymphocyte responses. The Journal of Experimental Medicine, 193, 375-386. (doi:10.1084/jem.193.3.375)
- 8. Klenerman, P., Meier, U. C., Phillips, R. E., & McMichael, A.J. 1995 The effects of natural altered peptide ligands on the whole blood cytotoxic T lymphocyte response to human immunodeficiency virus. European journal of immunology, 25, 1927-1931. (doi:10.1002/eji.1830250720)
- 9. Klenerman, P., Rowland-Jones, S., McAdam, S., Edwards, J., Daenke, S., Lalloo, D., Köppe, B., Rosenberg, W., Boyd, D., & Edwards, A. 1994 Cytotoxic T-cell activity antagonized by naturally occurring HIV-1 Gag variants. Nature, 369, 403-407. (doi:10.1038/369403a0)
- 10. Meier U. C., Klenerman, P., Griffin P., James, W., Köppe, B., Larder, B., McMichael, A. & Phillips, R. 1995 Cytotoxic T lymphocyte lysis inhibited by viable HIV mutants. Science, 270, 1360-1362. (doi:10.1126/science.270.5240.1360)
- 11. Price, D. A., Goulder, P. J. R., Klenerman, P., Sewell , A. K., Easterbrook, P. J., Troop, M., Bangham, C. R. M. & and Phillips, R. E. 1997 Positive selection of HIV-1 cytotoxic T lymphocyte escape variants during primaryinfection. Proceedings of the National Academy of Sciences, 94, 1890-1895.
- 12. Yokomaku, Y., Miura, H., Tomiyama, H., Kawana-Tachikawa, A., Takiguchi, M., Kojima, A., Nagai, Y., Iwamoto, A., Matsuda, Z. & Ariyoshi, K. 2004 Impaired processing and presentation of cytotoxic-T-lymphocyte (CTL) epitopes are major escape mechanisms from CTL immune pressure in human immunodeficiency virus type 1 infection. J. Virol., 78, 1324-1332. (doi:10.1128/JVI.78.3.1324-1332.2004)
- 13. Frater, A. J. et al. 2006 Passive sexual transmission of human immunodeficiency virus type 1 variants and adaptation in new hosts. J. Virol., 80, 7226-7234. (doi:10.1128/JVI.02014-05)
- 14. Goulder, P. J. R. & Walker, B. D. 2012 HIV and HLA class I: an evolving relationship. Immunity, 37, 426-440. (doi:10.1016/j.immuni.2012.09.005)
- 15. Kiepiela, P. et al. 2004 Dominant influence of HLA-B in mediating the potential co-evolution of HIV and HLA. Nature, 432, 769-775. (doi:10.1038/nature03113)
- 16. Frater, A. J. et al. 2007 Effective T-Cell responses select human immunodeficiency virus mutants and slow disease progression. Journal of Virology, 81, 6742-6751. (doi:10.1128/JVI.00022-07)
- 17. Leslie, A. J. et al. 2004 HIV evolution: CTL escape mutation and reversion after transmission. Nature Medicine, 10, 282-289. (doi:10.1038/nm992)
- 18. Moore, C. B., John, M., James, I. R., Christiansen, F. T., Witt, C. S. & Mallal, S. A. 2002 Evidence of HIV-1 adaptation to HLA-restricted immune responses at a population level. Science, 296, 1439-1443. (doi:10.1126/science.1069660)
- 19. Bhattacharya, T. et al. 2007 Founder effects in the assessment of HIV polymorphisms and HLA allele associations. Science, 315, 1583-1586. (doi:10.1126/science.1131528)
- 20. Felsenstein, J. 1981 Evolutionary trees from DNA sequences: a maximum likelihood approach. Journal of molecular evolution, 17, 368-376.
- 21. Fryer, H. R., Frater, A. J., Duda, A., Roberts, M. G., Phillips, R. E., McLean, A. R., and The SPARTAC Trial Investigators. 2010 Modelling the evolution and spread of HIV immune escape mutants. PLoS Pathog, 6, e1001196+. (doi:10.1371/journal.ppat.1001196)
- 22. Koelle, K, & Rasmussen, D. A., 2011 Rates of coalescence for common epidemiological models at equilibrium. Journal of The Royal Society Interface, 9, 997-1007. (doi:10.1098/rsif.2011.0495)
- 23. Rasmussen, D. A., Ratmann, O. & Koelle, K. 2011 Inference for nonlinear epidemiological models using genealogies and time series. PLoS Comput Biol, 7, e1002136+. (doi:10.1098/rsif.2011.0495)
- 24. Frost, S. D. W & Volz, E. M. 2010 Viral phylodynamics and the search for an ‘effective number of infections’. Philosophical Transactions of the Royal Society B: Biological Sciences, 365, 1879-1890. (doi:10.1098/rstb.2010.0060)
- 25. Goonetilleke, N. et al. 2009 The first T cell response to transmitted/founder virus contributes to the control of acute viremia in HIV-1 infection. The Journal of Experimental Medicine, 206, 1253-1272. (doi:10.1084/jem.20090365)
- 26. Fagard, C. et al. 2003 A prospective trial of structured treatment interruptions in human immunodeficiency virus infection. Archives of internal medicine, 163, 1220-1226. (doi:10.1001/archinte.163.10.1220)
- 27. Los Alamos HIV sequence database. http://www.hiv.lanl.gov/.
- 28. Gaschen, B., Kuiken, C., Korber, B. & Foley, B. 2001 Retrieval and on-the-fly alignment of sequence fragments from the HIV database. Bioinformatics (Oxford, England), 17, 415-418.
- 29. Drummond, A. J., Suchard, M. A., Xie, D. & Rambaut, A. 2012 Bayesian phylogenetics with BEAUti and the BEAST 1.7. Molecular Biology and Evolution, 29, 1969-1973. (doi:10.1093/molbev/mss075)
- 30. Yang, Z. & Rannala, B. 1997 Bayesian phylogenetic inference using DNA sequences: a Markov Chain Monte Carlo Method. Molecular Biology and Evolution, 14, 717-724.
- 31. Kendall, D. G., 1948 On the generalized “birth-and-death” process. The Annals of Mathematical Statistics, 19, 1-15.
- 32. Mahalanobis, P. C. 1936 On the generalised distance in statistics. Proceedings National Institute of Science, India, 2, 49-55.
- 33. Wand, M. P. Fast computation of multivariate kernel estimators. Journal of Computational and Graphical Statistics, 3, 433-445. (doi:10.1080/10618600.1994.10474656)
- 34. Wand, M. P. & Jones, M. C. 1994 Kernel Smoothing (Chapman & Hall/CRC Monographs on Statistics & Applied Probability). Chapman and Hall/CRC, 1st edition.
- 35. Swiss Confederation website. http://www.bag.admin.ch/hiv_aids/.
- 36. Carlson, J. M. et al. 2008 Phylogenetic dependency networks: inferring patterns of CTL escape and codon covariation in HIV-1 Gag. PLoS Comput Biol, 4, e1000225+. (doi:10.1371/journal.pcbi.1000225)
- 37. Carlson, J. M. et al. 2012 Widespread impact of HLA restriction on immune control and escape pathways of HIV-1. Journal of Virology, 86, 5230-5243. (doi:10.1128/JVI.06728-11)
- 38. Griffiths, R. C. & Tavaré, S. 1994 Sampling theory for neutral alleles in a varying environment. Philosophical transactions of the Royal Society of London. Series B, Biological sciences 344, 403-410. (doi:10.2307/56112)
- 39. Kingman, J. F. C. 1982 Exchangeability and the evolution of large populations. In G. Koch and F. Spizzichino, editors, Exchangeability in Probability and Statistics, pages 97-112. North-Holland, Amsterdam.
- 40. Slatkin, M. and Hudson, R. R. 1991 Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations. Genetics, 129, 555-562.
- 41. Stadler, T. 2009 On incomplete sampling under birth-death models and connections to the sampling-based coalescent. Journal of Theoretical Biology, 261, 58-66. (doi:10.1016/j.jtbi.2009.07.018)
- 42. Gernhard, T. 2008 The conditioned reconstructed process. Journal of Theoretical Biology, 253, 769-778. (doi:10.1016/j.jtbi.2008.04.005)
- 43. Nee, S., May, R. M. & Harvey, P. H. 1994 The reconstructed evolutionary process. Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences, 344, 305-311. (doi:10.1098/rstb.1994.0068)
- 44. Thompson, E. A. 1975 Human Evolutionary Trees. Cambridge University Press.
- 45. Jombart, T., Eggo, R. M., Dodd, P. J., & Balloux, F. 2011 Reconstructing disease outbreaks from genetic data: a graph approach. Heredity, 106, 383-390. (doi:10.1038/hdy.2010.78)
- 46. Stadler, T. 2010 Sampling-through-time in birth-death trees. Journal of Theoretical biology, 267, 396-404. (doi:10.1016/j.jtbi.2010.09.010)
- 47. Pybus, O. G. & Rambaut, A. Evolutionary analysis of the dynamics of viral infectious disease. Nat Rev Genet, 10, 540-550. (doi:10.1038/nrg2583)
- 48. Volz, E. M., Kosakovsky Pond, S. L., Ward, M. J., Leigh Brown, A. J., & Frost, S. D. W. 2009 Phylodynamics of infectious disease epidemics. Genetics, 183, 1421-1430. (doi:10.1534/genetics.109.106021)
- 49. Los Alamos HIV immunology database. http://www.hiv.lanl.gov/content/immunology.