Modelling tree shape and structure in viral phylodynamics

Abstract

Epidemiological models have highlighted the importance of population structure in the transmission dynamics of infectious diseases. Using HIV-1 as an example of a model evolutionary system, we consider how population structure affects the shape and the structure of a viral phylogeny in the absence of strong selection at the population level. For structured populations, the number of lineages as a function of time is insufficient to describe the shape of the phylogeny. We develop deterministic approximations for the dynamics of tips of the phylogeny over evolutionary time, the number of 'cherries', tips that share a direct common ancestor, and Sackin's index, a commonly used measure of phylogenetic imbalance or asymmetry. We employ cherries both as a measure of asymmetry of the tree as well as a measure of the association between sequences from different groups. We consider heterogeneity in infectiousness associated with different stages of HIV infection, and in contact rates between groups of individuals. In the absence of selection, we find that population structure may have relatively little impact on the overall asymmetry of a tree, especially when only a small fraction of infected individuals is sampled, but may have marked effects on how sequences from different subpopulations cluster and co-cluster.

Figure 1.

Schematic illustrating cherries and the calculation of Sackin's statistic for a symmetric (a) and an asymmetric (b) six-taxon cladogram.

Figure 2.

(a) Dynamics of the number of infected individuals, I (black line), and ν = Iτ/2 (dashed line) over time in weeks based on equations (3.1)–(3.2), as well as estimates of ‘scaled effective population size’ obtained from applying a Bayesian skyride (grey) to simulated data generated from a forwards-time stochastic version of the model, with 100 replicates. (b) Dynamics of the mean generation time, τ. Parameter values and initial conditions are as follows: β = 0.01, c = 1, , , , S(0) = 9999, I(0) = 1, with a simulation time of 40 years. Simulations of the differential equations were performed using the simecol library [29] in R [30], fitting of the skyline plot used the INLA library [31], while the stochastic simulations were performed using SimPy v. 1.9.1 in Python (see [3] for more details). Simulations were conditioned on reaching a quasi-equilibrium state, and registered by aligning the peaks of the simulated number of infected individuals to the peak of infected individuals from the ordinary differential equations. Code to perform simulations is available from http://code.google.com/p/simonfrost.

Figure 3.

Dynamics of (a) the number of lineages, A, (b) the distribution of tip lengths, (c) the mean cluster size, M, (d) the fraction of sequences clustered, 1−P, (e) the number of cherries, C, and (f) Sackin′s index, K, for the simple model of HIV infection given by equations (3.1)–(3.2). Parameter values, initial conditions, and simulations are as in figure 2.

Figure 4.

Asymmetry and clustering assuming a range of sampling fractions, from (red) to (violet) in steps of 0.1, for different values of the relative infectiousness of acute infection (a), and the relative contact rate in the differential risk model, assuming either (b) proportionate mixing (ρ = 0) or (c) preferential mixing (ρ = 0.9). Parameter values for the acute/chronic model are as follows: c = 1, , , , , S(0) = 9999, I₁(0) = 1, I₂(0) = 0. The infectivity parameters β_i were constrained such that and , where is the mean infectiousness (with ), d_i the mean duration of stage i and k the fold increase in infectiousness during acute infection. Parameter values for the differential risk model are β = 0.01, c₂ = 1, , , , , S₁(0) = 999, I₁(0) = 1, S₂(0) = 9000, I₂(0) = 0. The simulation time is 30 years, with weekly timesteps, assuming 52 weeks per year.