Viral phylodynamics and the search for an 'effective number of infections'

Abstract

Information on the dynamics of the effective population size over time can be obtained from the analysis of phylogenies, through the application of time-varying coalescent models. This approach has been used to study the dynamics of many different viruses, and has demonstrated a wide variety of patterns, which have been interpreted in the context of changes over time in the 'effective number of infections', a quantity proportional to the number of infected individuals. However, for infectious diseases, the rate of coalescence is driven primarily by new transmissions i.e. the incidence, and only indirectly by the number of infected individuals through sampling effects. Using commonly used epidemiological models, we show that the coalescence rate may indeed reflect the number of infected individuals during the initial phase of exponential growth when time is scaled by infectivity, but in general, a single change in time scale cannot be used to estimate the number of infected individuals. This has important implications when integrating phylogenetic data in the context of other epidemiological data.

Figure 1.

Schematic of different phylodynamic patterns for the relative size function, ν over time. (a) constant; (b) exponential; (c) piecewise expansion; (d) piecewise logistic; (e) constant–expansion–constant; (f) oscillatory.

Figure 2.

Phylodynamics of a simple susceptible-infected model in an open population (equations (2.7) and (2.8) in the main text). (a) Dynamics of the number of infected individuals, I over time in years. The vertical lines denote sampling times, and the number of lineages over time (b) during exponential growth (red), (c) following the peak of infected individuals (blue) and (d) at equilibrium (green). The grey lines represent stochastic simulations; in order to generate a fair comparison between the deterministic model and the stochastic simulations, time was shifted for each simulation such that the peak prevalence occurred at the same time as in the deterministic model. Parameter values are as follows (with time in years); βc = 52, γ = 1/10, μ = 1/70, Λ = 10000/70. The initial conditions were: X(0) = 9999, Y(0) = 1. Sampling times were set at 900/52, 2000/52 and 15000/52 years, and a sample size of 100 was assumed, i.e. A = 100. Numerical simulations were performed in R (R Development Core Team 2009) using the simecol library (Petzoldt & Rinke 2007). Stochastic simulations were performed with SimPy (http://simpy.sourceforge.net). All code is available from S.D.W.F. on request.

Figure 3.

The product of the transmission probability, contact rate and number of infected individuals, βcY, at different stages of the epidemic depicted in figure 2 (smooth black line) obtained from numerically solving equations (2.7) and (2.8), along with numerical estimates of ‘effective population size’ estimated using generalized skyline plots fitted to stochastic simulations (grey lines) on the same scale. During the exponential growth period, the skyline generates good estimates of βcY. Parameter values and initial conditions are as described in figure 2. Skyline plots were generated using the APE library (Paradis et al. 2004) in R (R Development Core Team 2009). (a) Exponential growth; (b) after peak; (c) at equilibrium.

Figure 4.

(a) Dynamics of the number of infected individuals, Y, and the transmission rate f_XY for a susceptible-infected-recovered model with seasonal forcing, given by equations (2.20)–(2.22) in the main text. Parameter values are as follows (with time in days): β₀ = 10/7, β₁ = 0.05, ω = 2π/365, γ = 1/7, μ = 1/25550. The population size, N, was assumed to be 10⁶. Initial conditions: S = 100029.946, I = 142.978, R = 899827.076. (b) The time of sampling for the high prevalence scenario was t = 3465, when I = 384.477 (38.4 per 100 000) (c) and the time of sampling for the low prevalence scenario was t = 3649, when I = 140.7068 (14.0 per 100000).