Viral phylodynamics is the study of how epidemiological, immunological, and evolutionary processes act and potentially interact to shape viral phylogenies. Since the term was coined in 2004, research on viral phylodynamics has focused on transmission dynamics in an effort to shed light on how these dynamics impact viral genetic variation. Transmission dynamics can be considered at the level of cells within an infected host, individual hosts within a population, or entire populations of hosts. Many viruses, especially RNA viruses, rapidly accumulate genetic variation because of short generation times and high mutation rates. Patterns of viral genetic variation are therefore heavily influenced by how quickly transmission occurs and by which entities transmit to one another. Patterns of viral genetic variation will also be affected by selection acting on viral phenotypes. Although viruses can differ with respect to many phenotypes, phylodynamic studies have to date tended to focus on a limited number of viral phenotypes. These include virulence phenotypes, phenotypes associated with viral transmissibility, cell or tissue tropism phenotypes, and antigenic phenotypes that can facilitate escape from host immunity. Due to the impact that transmission dynamics and selection can have on viral genetic variation, viral phylogenies can therefore be used to investigate important epidemiological, immunological, and evolutionary processes, such as epidemic spread, spatio-temporal dynamics including metapopulation dynamics, zoonotic transmission, tissue tropism, and antigenic drift. The quantitative investigation of these processes through the consideration of viral phylogenies is the central aim of viral phylodynamics.
Sources of phylodynamic variation
In coining the term phylodynamics, Grenfell and coauthors postulated that viral phylogenies "... are determined by a combination of immune selection, changes in viral population size, and spatial dynamics". Their study showcased three features of viral phylogenies, which may serve as rules of thumb for identifying important epidemiological, immunological, and evolutionary processes influencing patterns of viral genetic variation.
The relative lengths of internal versus external branches will be affected by changes in viral population size over time Rapid expansion of a virus in a population will be reflected by a "star-like" tree, in which external branches are long relative to internal branches. Star-like trees arise because viruses are more likely to share a recent common ancestor when the population is small, and a growing population has an increasingly smaller population size towards the past. Compared to a phylogeny of an expanding virus, a phylogeny of a viral population that stays constant in size will have external branches that are shorter relative to branches on the interior of the tree. The phylogeny of HIV provides a good example of a star-like tree, as the prevalence of HIV infection rose rapidly throughout the 1980s (exponential growth). The phylogeny of hepatitis B virus instead reflects a viral population that has remained roughly consistent (constant size). Similarly, trees reconstructed from viral sequences isolated from chronically infected individuals can be used to gauge changes in viral population sizes within a host. The clustering of taxa on a viral phylogeny will be affected by host population structure Viruses within similar hosts, such as hosts that reside in the same geographic region, are expected to be more closely related genetically if transmission occurs more commonly between them. The phylogenies of measles and rabies virus illustrate viruses with spatially structured host population. These phylogenies stand in contrast to the phylogeny of human influenza, which does not appear to exhibit strong spatial structure over extended periods of time. Clustering of taxa, when it occurs, is not necessarily observed at all scales, and a population that appears structured at some scale may appear panmictic at another scale, for example at a smaller spatial scale. While spatial structure is the most commonly observed population structure in phylodynamic analyses, viruses may also have nonrandom admixture by attributes such as the age, race, and risk behavior. This is because viral transmission can preferentially occur between hosts sharing any of these attributes. Tree balance will be affected by selection, most notably immune escape The effect of directional selection on the shape of a viral phylogeny is exemplified by contrasting the trees of influenza virus and HIV's surface proteins. The ladder-like phylogeny of influenza virus A/H3N2's hemagglutinin protein bears the hallmarks of strong directional selection, driven by immune escape (imbalanced tree). In contrast, a more balanced phylogeny may occur when a virus is not subject to strong immune selection or other source of directional selection. An example of this is the phylogeny of the HIV envelope protein inferred from sequences isolated from different individuals in a population (balanced tree). Phylogenies of the HIV envelope protein from chronically infected hosts resemble influenza's ladder-like tree. This highlights that the processes affecting viral genetic variation can differ across scales. Indeed, contrasting patterns of viral genetic variation within and between hosts has been an active topic in phylodynamic research since the field's inception. Although these three phylogenetic features are useful rules of thumb to identify epidemiological, immunological, and evolutionary processes that might be impacting viral genetic variation, there is growing recognition that the mapping between process and phylogenetic pattern can be many-to-one. For instance, although ladder-like trees could reflect the presence of directional selection, ladder-like trees could also reflect sequential genetic bottlenecks that might occur with rapid spatial spread, as in the case of rabies virus. Because of this many-to-one mapping between process and phylogenetic pattern, research in the field of viral phylodynamics has sought to develop and apply quantitative methods to effectively infer process from reconstructed viral phylogenies (see Methods). The consideration of other data sources (e.g., incidence patterns) may aid in distinguishing between competing phylodynamic hypotheses. Combining disparate sources of data for phylodynamic analysis remains a major challenge in the field and is an active area of research.
Applications
Viral origins Phylodynamic models may aid in dating epidemic and pandemic origins. The rapid rate of evolution in viruses allows molecular clock models to be estimated from genetic sequences, thus providing a per-year rate of evolution of the virus. With the rate of evolution measured in real units of time, it is possible to infer the date of the most recent common ancestor (MRCA) for a set of viral sequences. The age of the MRCA of these isolates is a lower bound; the common ancestor of the entire virus population must have existed earlier than the MRCA of the virus sample. In April 2009, genetic analysis of 11 sequences of swine-origin H1N1 influenza suggested that the common ancestor existed at or before 12 January 2009. This finding aided in making an early estimate of the basic reproduction number R 0 {\displaystyle R_{0}} of the pandemic. Similarly, genetic analysis of sequences isolated from within an individual can be used to determine the individual's infection time.
Viral spread Phylodynamic models may provide insight into epidemiological parameters that are difficult to assess through traditional surveillance means. For example, assessment of R 0 {\displaystyle R_{0}} from surveillance data requires careful control of the variation of the reporting rate and the intensity of surveillance. Inferring the demographic history of the virus population from genetic data may help to avoid these difficulties and can provide a separate avenue for inference of R 0 {\displaystyle R_{0}} . Such approaches have been used to estimate R 0 {\displaystyle R_{0}} in hepatitis C virus and HIV. Additionally, differential transmission between groups, be they geographic-, age-, or risk-related, is very difficult to assess from surveillance data alone. Phylogeographic models have the possibility of more directly revealing these otherwise hidden transmission patterns. Phylodynamic approaches have mapped the geographic movement of the human influenza virus and quantified the epidemic spread of rabies virus in North American raccoons. However, nonrepresentative sampling may bias inferences of both R 0 {\displaystyle R_{0}} and migration patterns. Phylodynamic approaches have also been used to better understand viral transmission dynamics and spread within infected hosts. For example, phylodynamic studies have been used to infer the rate of viral growth within infected hosts and to argue for the occurrence of viral compartmentalization in hepatitis C infection.
Viral control efforts Phylodynamic approaches can also be useful in ascertaining the effectiveness of viral control efforts, particularly for diseases with low reporting rates. For example, the genetic diversity of the DNA-based hepatitis B virus declined in the Netherlands in the late 1990s, following the initiation of a vaccination program. This correlation was used to argue that vaccination was effective at reducing the prevalence of infection, although alternative explanations are possible. Viral control efforts can also impact the rate at which virus populations evolve, thereby influencing phylogenetic patterns. Phylodynamic approaches that quantify how evolutionary rates change over time can therefore provide insight into the effectiveness of control strategies. For example, an application to HIV sequences within infected hosts showed that viral substitution rates dropped to effectively zero following the initiation of antiretroviral drug therapy. This decrease in substitution rates was interpreted as an effective cessation of viral replication following the commencement of treatment, and would be expected to lead to lower viral loads. This finding is especially encouraging because lower substitution rates are associated with slower progression to AIDS in treatment-naive patients. Antiviral treatment also creates selective pressure for the evolution of drug resistance in virus populations, and can thereby affect patterns of genetic diversity. Commonly, there is a fitness trade-off between faster replication of susceptible strains in the absence of antiviral treatment and faster replication of resistant strains in the presence of antivirals. Thus, ascertaining the level of antiviral pressure necessary to shift evolutionary outcomes is of public health importance. Phylodynamic approaches have been used to examine the spread of oseltamivir resistance in influenza A/H1N1.
Methods Most often, the goal of phylodynamic analyses is to make inferences of epidemiological processes from viral phylogenies. Thus, most phylodynamic analyses begin with the reconstruction of a phylogenetic tree. Genetic sequences are often sampled at multiple time points, which allows the estimation of substitution rates and the time of the MRCA using a molecular clock model. For viruses, Bayesian phylogenetic methods are popular because of the ability to fit complex demographic scenarios while integrating out phylogenetic uncertainty. Traditional evolutionary approaches directly utilize methods from computational phylogenetics and population genetics to assess hypotheses of selection and population structure without direct regard for epidemiological models. For example,
the magnitude of selection can be measured by comparing the rate of nonsynonymous substitution to the rate of synonymous substitution (dN/dS); the population structure of the host population may be examined by calculation of F-statistics; and hypotheses concerning panmixis and selective neutrality of the virus may be tested with statistics such as Tajima's D. However, such analyses were not designed with epidemiological inference in mind and it may be difficult to extrapolate from standard statistics to desired epidemiological quantities. In an effort to bridge the gap between traditional evolutionary approaches and epidemiological models, several analytical methods have been developed to specifically address problems related to phylodynamics. These methods are based on coalescent theory, birth-death models, and simulation, and are used to more directly relate epidemiological parameters to observed viral sequences.
Coalescent theory and phylodynamics
Effective population size The coalescent is a mathematical model that describes the ancestry of a sample of nonrecombining gene copies. In modeling the coalescent process, time is usually considered to flow backwards from the present. In a selectively neutral population of constant size N {\displaystyle N} and nonoverlapping generations (the Wright Fisher model), the expected time for a sample of two gene copies to coalesce (i.e., find a common ancestor) is N {\displaystyle N} generations. More generally, the waiting time for two members of a sample of n {\displaystyle n} gene copies to share a common ancestor is exponentially distributed, with rate
λ n = ( n 2 ) 1 N {\displaystyle \lambda _{n}={n \choose 2}{\frac {1}{N}}} . This time interval is labeled T n {\displaystyle T_{n}} , and at its end there are n − 1 {\displaystyle n-1} extant lineages remaining. These remaining lineages will coalesce at the rate λ n − 1 ⋯ λ 2 {\displaystyle \lambda _{n-1}\cdots \lambda _{2}} after intervals T n − 1 ⋯ T 2 {\displaystyle T_{n-1}\cdots T_{2}} . This process can be simulated by drawing exponential random variables with rates { λ n − i } i = 0 , ⋯ , n − 2 {\displaystyle \{\lambda _{n-i}\}_{i=0,\cdots ,n-2}} until there is only a single lineage remaining (the MRCA of the sample). In the absence of selection and population structure, the tree topology may be simulated by picking two lineages uniformly at random after each coalescent interval T i {\displaystyle T_{i}} .
The expected waiting time to find the MRCA of the sample is the sum of the expected values of the internode intervals,
E [ T M R C A ] = E [ T n ] + E [ T n − 1 ] + ⋯ + E [ T 2 ] = 1 / λ n + 1 / λ n − 1 + ⋯ + 1 / λ 2 = 2 N ( 1 − 1 n ) . {\displaystyle {\begin{aligned}\mathrm {E} [\mathrm {TMRCA} ]&=\mathrm {E} [T_{n}]+\mathrm {E} [T_{n-1}]+\cdots +\mathrm {E} [T_{2}]\\&=1/\lambda _{n}+1/\lambda _{n-1}+\cdots +1/\lambda _{2}\\&=2N(1-{\frac {1}{n}}).\end{aligned}}}
Two corollaries are :
The time to the MRCA (TMRCA) of a sample is not unbounded in the sample size. lim n → ∞ E [ T M R C A ] = 2 N . {\displaystyle \lim _{n\rightarrow \infty }\mathrm {E} [\mathrm {TMRCA} ]=2N.}
Few samples are required for the expected TMRCA of the sample to be close to the theoretical upper bound, as the difference is O ( 1 / n ) {\displaystyle O(1/n)} . Consequently, the TMRCA estimated from a relatively small sample of viral genetic sequences is an asymptotically unbiased estimate for the time that the viral population was founded in the host population. For example, Robbins et al. estimated the TMRCA for 74 HIV-1 subtype-B genetic sequences collected in North America to be 1968. Assuming a constant population size, we expect the time back to 1968 to represent 1 − 1 / 74 = 99 % {\displaystyle 1-1/74=99\%} of the TMRCA of the North American virus population. If the population size N ( t ) {\displaystyle N(t)} changes over time, the coalescent rate λ n ( t ) {\displaystyle \lambda _{n}(t)} will also be a function of time. Donnelley and Tavaré derived this rate for a time-varying population size under the assumption of constant birth rates:
λ n ( t ) = ( n 2 ) 1 N ( t ) {\displaystyle \lambda _{n}(t)={n \choose 2}{\frac {1}{N(t)}}} . Because all topologies are equally likely under the neutral coalescent, this model will have the same properties as the constant-size coalescent under a rescaling of the time variable: t → ∫ τ = 0 t d τ N ( τ ) {\displaystyle t\rightarrow \int _{\tau =0}^{t}{\frac {\mathrm {d} \tau }{N(\tau )}}} . Very early in an epidemic, the virus population may be growing exponentially at rate r {\displaystyle r} , so that t {\displaystyle t} units of time in the past, the population will have size N ( t ) = N 0 e − r t {\displaystyle N(t)=N_{0}e^{-rt}} . In this case, the rate of coalescence becomes
λ n ( t ) = ( n 2 ) 1 N 0 e − r t {\displaystyle \lambda _{n}(t)={n \choose 2}{\frac {1}{N_{0}e^{-rt}}}} . This rate is small close to when the sample was collected ( t = 0 {\displaystyle t=0} ), so that external branches (those without descendants) of a gene genealogy will tend to be long relative to those close to the root of the tree. This is why rapidly growing populations yield trees with long tip branches. If the rate of exponential growth is estimated from a gene genealogy, it may be combined with knowledge of the duration of infection or the serial interval D {\displaystyle D} for a particular pathogen to estimate the basic reproduction number, R 0 {\displaystyle R_{0}} . The two may be linked by the following equation:
r = R 0 − 1 D {\displaystyle r={\frac {R_{0}-1}{D}}} . For example, one of the first estimates of R 0 {\displaystyle R_{0}} was for pandemic H1N1 influenza in 2009 by using a coalescent-based analysis of 11 hemagglutinin sequences in combination with prior data about the infectious period for influenza.
Compartmental models Infectious disease epidemics are often characterized by highly nonlinear and rapid changes in the number of infected individuals and the effective population size of the virus. In such cases, birth rates are highly variable, which can diminish the correspondence between effective population size and the prevalence of infection. Many mathematical models have been developed in the field of mathematical epidemiology to describe the nonlinear time series of prevalence of infection and the number of susceptible hosts. A well studied example is the Susceptible-Infected-Recovered (SIR) system of differential equations, which describes the fractions of the population S ( t ) {\displaystyle S(t)} susceptible, I ( t ) {\displaystyle I(t)} infected, and R ( t ) {\displaystyle R(t)} recovered as a function of time:
d S d t = − β S I {\displaystyle {\frac {dS}{dt}}=-\beta SI} ,
d I d t = β S I − γ I {\displaystyle {\frac {dI}{dt}}=\beta SI-\gamma I} , and
d R d t = γ I {\displaystyle {\frac {dR}{dt}}=\gamma I} . Here, β {\displaystyle \beta } is the per capita rate of transmission to susceptible hosts, and γ {\displaystyle \gamma } is the rate at which infected individuals recover, whereupon they are no longer infectious. In this case, the incidence of new infections per unit time is f ( t ) = β S I {\displaystyle f(t)=\beta SI} , which is analogous to the birth rate in classical population genetics models. The general formula for the rate of coalescence is:
λ n ( t ) = ( n 2 ) 2 f ( t ) I ( t ) 2 {\displaystyle \lambda _{n}(t)={n \choose 2}{\frac {2f(t)}{I(t)^{2}}}} . The ratio 2 ( n 2 ) / I ( t ) 2 {\displaystyle 2{n \choose 2}/{I(t)^{2}}} can be understood as arising from the probability that two lineages selected uniformly at random are both ancestral to the sample. This probability is the ratio of the number of ways to pick two lineages without replacement from the set of lineages and from the set of all infections: ( n 2 ) / ( I ( t ) 2 ) ≈ 2 ( n 2 ) / I ( t ) 2 {\displaystyle {n \choose 2}/{I(t) \choose 2}\approx 2{n \choose 2}/{I(t)^{2}}} . Coalescent events will occur with this probability at the rate given by the incidence function f ( t ) {\displaystyle f(t)} . For the simple SIR model, this yields
λ n ( t ) = ( n 2 ) 2 β S ( t ) I ( t ) {\displaystyle \lambda _{n}(t)={n \choose 2}{\frac {2\beta S(t)}{I(t)}}} . This expression is similar to the Kingman coalescent rate, but is damped by the fraction susceptible S ( t ) {\displaystyle S(t)} . Early in an epidemic, S ( 0 ) ≈ 1 {\displaystyle S(0)\approx 1} , so for the SIR model
λ n ( t ) ≈ ( n 2 ) 2 β I ( t ) {\displaystyle \lambda _{n}(t)\approx {n \choose 2}{\frac {2\beta }{I(t)}}} . This has the same mathematical form as the rate in the Kingman coalescent, substituting N e = I ( t ) / 2 β {\displaystyle N_{e}=I(t)/2\beta } . Consequently, estimates of effective population size based on the Kingman coalescent will be proportional to prevalence of infection during the early period of exponential growth of the epidemic. When a disease is no longer exponentially growing but has become endemic, the rate of lineage coalescence can also be derived for the epidemiological model governing the disease's transmission dynamics. This can be done by extending the Wright Fisher model to allow for unequal offspring distributions. With a Wright Fisher generation taking τ {\displaystyle \tau } units of time, the rate of coalescence is given by:
λ n = ( n 2 ) 1 N e τ {\displaystyle \lambda _{n}={n \choose 2}{\frac {1}{N_{e}\tau }}} , where the effective population size N e {\displaystyle N_{e}} is the population size N {\displaystyle N} divided by the variance of the offspring distribution σ 2 {\displaystyle \sigma ^{2}} . The generation time τ {\displaystyle \tau } for an epidemiological model at equilibrium is given by the duration of infection and the population size N {\displaystyle N} is closely related to the equilibrium number of infected individuals. To derive the variance in the offspring distribution σ 2 {\displaystyle \sigma ^{2}} for a given epidemiological model, one can imagine that infected individuals can differ from one another in their infectivities, their contact rates, their durations of infection, or in other characteristics relating to their ability to transmit the virus with which they are infected. These differences can be acknowledged by assuming that the basic reproduction number is a random variable ν {\displaystyle \nu } that varies across individuals in the population and that ν {\displaystyle \nu } follows some continuous probability distribution. The mean and variance of these individual basic reproduction numbers, E [ ν ] {\displaystyle \mathrm {E} [\nu ]} and V a r [ ν ] {\displaystyle \mathrm {Var} [\nu ]} , respectively, can then be used to compute σ 2 {\displaystyle \sigma ^{2}} . The expression relating these quantities is given by:
σ 2 = V a r [ ν ] E [ ν ] 2 + 1 {\displaystyle \sigma ^{2}={\frac {\mathrm {Var} [\nu ]}{\mathrm {E} [\nu ]^{2}}}+1} . For example, for the SIR model above, modified to include births into the population and deaths out of the population, the population size N {\displaystyle N} is given by the equilibrium number of infected individuals, I {\displaystyle I} . The mean basic reproduction number, averaged across all infected individuals, is given by β / γ {\displaystyle \beta /\gamma } , under the assumption that the background mortality rate is negligible compared to the rate of recovery γ {\displaystyle \gamma } . The variance in individuals' basic reproduction rates is given by ( β / γ ) 2 {\displaystyle (\beta /\gamma )^{2}} , because the duration of time individuals remain infected in the SIR model is exponentially distributed. The variance in the offspring distribution σ 2 {\displaystyle \sigma ^{2}} is therefore 2. N e {\displaystyle N_{e}} therefore becomes I 2 {\displaystyle {\frac {I}{2}}} and the rate of coalescence becomes:
λ n = ( n 2 ) 2 γ I {\displaystyle \lambda _{n}={n \choose 2}{\frac {2\gamma }{I}}} . This rate, derived for the SIR model at equilibrium, is equivalent to the rate of coalescence given by the more general formula. Rates of coalescence can similarly be derived for epidemiological models with superspreaders or other transmission heterogeneities, for models with individuals who are exposed but not yet infectious, and for models with variable infectious periods, among others. Given some epidemiological information (such as the duration of infection) and a specification of a mathematical model, viral phylogenies can therefore be used to estimate epidemiological parameters that might otherwise be difficult to quantify.
Phylogeography At the most basic level, the presence of geographic population structure can be revealed by comparing the genetic relatedness of viral isolates to geographic relatedness. A basic question is whether geographic character labels are more clustered on a phylogeny than expected under a simple nonstructured model. This question can be answered by counting the number of geographic transitions on the phylogeny via parsimony, maximum likelihood or through Bayesian inference. If population structure exists, then there will be fewer geographic transitions on the phylogeny than expected in a panmictic model. This hypothesis can be tested by randomly scrambling the character labels on the tips of the phylogeny and counting the number of geographic transitions present in the scrambled data. By repeatedly scrambling the data and calculating transition counts, a null distribution can be constructed and a p-value computed by comparing the observed transition counts to this null distribution. Beyond the presence or absence of population structure, phylodynamic methods can be used to infer the rates of movement of viral lineages between geographic locations and reconstruct the geographic locations of ancestral lineages. Here, geographic location is treated as a phylogenetic character state, similar in spirit to 'A', 'T', 'G', 'C', so that geographic location is encoded as a substitution model. The same phylogenetic machinery that is used to infer models of DNA evolution can thus be used to infer geographic transition matrices. The end result is a rate, me