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

E-mail: duncan.palmer@stats.ox.ac.uk


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 [18], though these results have been called into question [19], 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 [20] 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 [21] (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].

Figure 1: The inference problem. The cartoon in panel displays the dynamic processes which may occur along a branch within the transmission tree. Time increases from left to right. Susceptible individuals are shown in blue. An individual infected with a virus which is wild-type at the epitope under consideration is green, and individuals with a viral strain with the escape mutation are red. Lighter and darker colours indicate the HLA status of the host. Darker indicates an HLA matched host, and lighter an HLA mismatched host. From left to right, transmission occurs within the population. Transmission of a strain which is wild-type at the epitope under consideration may escape. Transmission of an escaped viral strain to an individual who is HLA mismatched may result in reversion of the virus to wild-type at the epitope. Viruses exist within these environments over their evolutionary history. Thus, associated to a collection of individuals sampled at the present (shown at the tips in ) is a colour coded transmission tree, illustrated in . A transmission event is associated with each coalescence, but due to incomplete sampling, unseen transmission events also occur. These are shown by black crosses. This sampled transmission tree is embedded in the full transmission tree, shown in the second panel of . We have sequence data and colourings at the tips of a sampled transmission tree. Using these data we hope to reconstruct the embedded tree in , and use this reconstruction to make inferences about the unknown full transmission process (shown in grey in ).

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 [21]. 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 [21]. 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 [25], 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 [21]. 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 [27] to this HLA typed cross-sectional sequence data, we create a DNA multiple sequence alignment [28] and perform some data trimming. The alignment is then passed to the program BEAST [29] 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 [20]. 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 [30]. 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 [21]. 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 [31] 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 [21]. 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 [32] (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.

Figure 2: Simulation results. Panel shows two instances of simulations under the integrated approach and the ODE model, with the truth set at . 10,000 bootstraps of the data are applied and estimates under the ODE model shown as dots, and coloured according to Mahalanobis distance (the furthest and in Mahalanobis distance are coloured red and orange respectively. The remainder is coloured green - see methods). The ODE point estimate of the tip data is the blue dot. and credible regions of our integrated method are shown in black. The MAP is shown as a black dot. The MAP assuming we know the true underlying tree is the circled green dot. The truth is shown as a purple triangle. Panel shows 1000 simulations assuming the true underlying tree is known, under the ODE model and our integrated approach. Histograms of the ODE point estimates and MAP estimates are shown in blue and green repectively. The truth is overlaid in red. The first and second columns are histograms of estimates of the underlying escape and reversion rates respectively. Rates and are grouped in the histograms. Panel shows the results of 100 simulations of our full model applied to data generated from underlying rates , shown as a white dot. The 100 MAP estimates are shown in red, and contours are coloured using a 2D kernel density estimate [33, 34].

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 [35]. 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 [21] 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 [29] 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.

Figure 3: Panel displays the maximum clade credibility tree for cross-sectional data for epitopes and , with tips coloured as in Figure 1. Credible regions for the four epitopes under consideration, and corresponding consensus trees. Panels and show and credible regions under our integrated method. Panel is a zoom in of the rectangle in panel . Axes are linear on and on the log scale for values , and colourings correspond to the upper right figure legend of . Coloured diamonds are our MAP estimates of . Estimates obtained through the ODE method are displayed as filled circles.
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
TAFTIPSI pol (128-135) B51 I8T 0.126 0.85 5.17 17.8 35.6
RPMTYKAAV nef (77-85) B7 T4S, Y5F 0.166 69.4 37.8 44.3 26.0
Table 1: Summary table of escape mutations analysed. Location is defined by the HXB2 B-clade reference sequence. is the HLA prevalence in Caucasians [citeulike:7931760]. Time to escape and reversion is measured in .

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 [18] of viral sequences. This was followed by partial consideration of the viral genealogy [19] 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 [21]. 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 [21], 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 [41]. 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 [45]. 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 [46], 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 [47]. 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 [48] 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 [49] 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.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description