Phylodynamic Inference with Kernel ABC and Its Application to HIV Epidemiology

Abstract

The shapes of phylogenetic trees relating virus populations are determined by the adaptation of viruses within each host, and by the transmission of viruses among hosts. Phylodynamic inference attempts to reverse this flow of information, estimating parameters of these processes from the shape of a virus phylogeny reconstructed from a sample of genetic sequences from the epidemic. A key challenge to phylodynamic inference is quantifying the similarity between two trees in an efficient and comprehensive way. In this study, I demonstrate that a new distance measure, based on a subset tree kernel function from computational linguistics, confers a significant improvement over previous measures of tree shape for classifying trees generated under different epidemiological scenarios. Next, I incorporate this kernel-based distance measure into an approximate Bayesian computation (ABC) framework for phylodynamic inference. ABC bypasses the need for an analytical solution of model likelihood, as it only requires the ability to simulate data from the model. I validate this "kernel-ABC" method for phylodynamic inference by estimating parameters from data simulated under a simple epidemiological model. Results indicate that kernel-ABC attained greater accuracy for parameters associated with virus transmission than leading software on the same data sets. Finally, I apply the kernel-ABC framework to study a recent outbreak of a recombinant HIV subtype in China. Kernel-ABC provides a versatile framework for phylodynamic inference because it can fit a broader range of models than methods that rely on the computation of exact likelihoods.

Keywords: approximate Bayesian computation; human immunodeficiency virus; molecular epidemiology; phylodynamics; tree shape; virus evolution.

Fig. 1.

(A) Response of Sackin’s index (I_S) to the varying heterogeneity in contact rates. Each set of box-and-whisker plots summarizes the distribution of I_S for 100 replicate trees simulated under different values of c₁ ($c_1=c_2=1$ is shaded for reference). Varying c₁ away from $c_2=1$ had a more pronounced effect under preferential ($\rho =0.9$) than proportional (ρ = 0) mixing. (B) Projection of simulated trees to kernel space. Sets of 100 replicate trees are each annotated in the plot with their corresponding c₁ values. The proportion of variation explained by the first two principal components, as estimated from the eigenvalues, is reported by the respective axis labels.

Fig. 2.

Comparison of relative errors for posterior estimates of BDSIR model parameters obtained by kernel-ABC and BEAST2. Each box-and-whisker plot summarizes the parameter distribution sampled from the posterior density. The relative error was calculated as $(\widehat{x}-x)/x$ where $\widehat{x}$ is the estimated value and x is the true value; hence, a relative error of −0.9 corresponds to a 10-fold underestimate. These quantities were offset by 1 for rendering on a log-transformed scale. The vertical span of boxes corresponds to the interquartile range and the whiskers extend to the empirical 95% credible intervals. Population size corresponds to the total number of susceptible, infected and removed individuals. Detailed results are provided in supplementary table S1, Supplementary Material online.

Fig. 3.

Histograms summarizing the distribution of c₁ (contact rate, group 1) estimated by a kernel-ABC analysis on trees simulated with $c_1=0.5$ (blue) and $c_1=2.0$ (red, hatched). The true parameter settings and median estimates are indicated by solid and dashed white vertical lines, respectively. Both runs were initialized at $c_1=1$. The simulated annealing tolerance parameter was set to decay from 0.005 to 0.002.

Fig. 4.

Reconstructed dynamics of a recent outbreak of HIV CRF07_BC in China. Numbers of susceptible (gray) and infected (red) individuals over time were obtained by forward-time simulation in MASTER under the birth–death SIR model, using parameter values sampled by kernel-ABC. Shaded areas indicate the interquartile regions of the respective counts. Note that the infected counts include removed individuals from simulated populations. A solid vertical line indicates the median tree height that approximates the date of the most recent sample in the data set (2010). Dashed lines indicate the interquartile range in tree heights. Simulation time zero corresponds to the estimated origin of the epidemic (1993–1994).