A key priority in infectious disease research is to understand the ecological and evolutionary drivers of viral diseases from data on disease incidence as well as viral genetic and antigenic variation. We propose using a simulation-based, Bayesian method known as Approximate Bayesian Computation (ABC) to fit and assess phylodynamic models that simulate pathogen evolution and ecology against summaries of these data. We illustrate the versatility of the method by analyzing two spatial models describing the phylodynamics of interpandemic human influenza virus subtype A (H3N2). The first model captures antigenic drift phenomenologically with continuously waning immunity, and the second epochal evolution model describes the replacement of major, relatively long-lived antigenic clusters. Combining features of long-term surveillance data from the Netherlands with features of influenza A (H3N2) hemagglutinin gene sequences sampled in northern Europe, key phylodynamic parameters can be estimated with ABC. Goodness-of-fit analyses reveal that the irregularity in interannual incidence and H3N2' s ladder-like hemagglutinin phylogeny are quantitatively only reproduced under the epochal evolution model within a spatial context. However, the concomitant incidence dynamics result in a very large reproductive number and are not consistent with empirical estimates of H3N2' s population level attack rate. These results demonstrate that the interactions between the evolutionary and ecological processes impose multiple quantitative constraints on the phylodynamic trajectories of influenza A (H3N2), so that sequence and surveillance data can be used synergistically. ABC, one of several data synthesis approaches, can easily interface a broad class of phylodynamic models with various types of data but requires careful calibration of the summaries and tolerance parameters.
Many infectious pathogens, most notably RNA viruses, evolve on the same time scale as their ecological dynamics [1]. One of the perhaps best documented examples are human influenza A viruses, which cause substantial morbidity and mortality as they escape host immunity predominantly through the evolution of their surface antigens [2]. The resulting, dynamical interaction between the ecological and evolutionary processess can be better understood through the formulation and simulation of so-called “phylodynamic” mathematical models, e. g. [3]–[8]. However, while data on disease incidence as well as viral genetic and antigenic variation are increasing for many viruses, e. g. [9]–[13], fitting and assessing phylodynamic models to these data is still not commonly done. Historically, epidemiological time series data have been pervasively used to analyze hypotheses of host-pathogen interactions at the population level [14]–[17]. However, time series data capture the underlying evolutionary processes of pathogens only very indirectly. For flu, this has limited the type of infectious disease models that can be statistically interfaced with time series data, and the number of epidemiological parameters that can be simultaneously estimated [18], [19]. Consequently, the disease behavior of rapidly evolving pathogens is increasingly studied under additional, complementary data sets [1], most typically in ways that attempt to qualitatively reproduce prominent disease attributes [3]–[8]. More recently, coalescent-based statistical methods have been used to elucidate the disease dynamics of RNA viruses from molecular genetic data alone [20]. These methods have been particularly useful to reconstruct epidemiological transmission histories, identifying when and where transmission occurred and how viral populations change over time. For example, coalescent-based analyses have highlighted the importance of the tropics in the complex circulation dynamics of human influenza A (H3N2) virus (in short: H3N2) [9], [21], [22]. However, most coalescent methods estimate past population dynamics within a class of flexible demographic functions including exponential and logistic growth as well as the nonparametric Bayesian skyride [23], [24]; but see also [25]. These demographic functions do not explicitly describe the non-linear population dynamics of RNA viruses. Thus, assessing which ecological interactions underlie observed patterns of sequence diversity, and estimating the respective strength of these interactions, is difficult within this framework. Because of these limitations, we adopt a different statistical approach known as Approximate Bayesian Computation (ABC) to infer the phylodynamics of RNA viruses. ABC allows mechanistic phylodynamic models to be simultaneously fitted against both sequence and surveillance data. This method circumvents explicit likelihood calculations by simulating instead from the stochastic model that defines the likelihood [26]. Recent extensions of ABC allow for model assessment to be carried out at no further computational cost [27]. We further suggest incorporating variable selection procedures to quantify if and to what extent the data provide support for the inclusion of specific model components [28]. To demonstrate the utility of our approach, we consider the phylodynamics of interpandemic H3N2. We obtained weekly reports of H3N2 incidence in the Netherlands from 1994–2009 by combining influenza-like-illness (ILI) surveillance data with detailed records of associated, laboratory-confirmed cases of flu by type and subtype [29], [30], and similarly for France and the USA; see Figure 1 and the supplementary online material (Text S1). In addition, we reconstructed the ladder-like phylogeny of H3N2' s haemagglutinin gene (HA) from dated European sequences collected in 1968–2009 (see Figure 1 and Text S1). To represent H3N2' s global phylodynamics, we focus on a class of spatially structured phylodynamic compartmental models that formalize probabilistically how evolving, antigenic variants interact epidemiologically. These antigenic variants might correspond to the major antigenic clusters that are distinguishable in H3N2 antigenic maps [31], but can in principle also represent a different phenotypic resolution. The evolutionary dynamics of viral genotypes are separately formulated for each antigenic phenotype because genetic distances do not necessarily easily translate into phenotypic relationships [5]. Spatial substructure has been incorporated in several models of H3N2 phylodynamics to reflect the global circulation of the virus [4], [8], [32]. We adopt here a simple source-sink framework, where the sink is thought of as the Netherlands into which viral genetic diversity and antigenic strains are imported on a seasonal scale from a source population where the virus persists [9], [33]. We fit and assess two distinct models to the combined features of sequence and incidence data described in Figure 1 and Table 1. The first model captures H3N2' s antigenic drift phenomenologically through gradual loss of immunity, and the second model describes the antigenic evolution of the virus explicitly with particular assumptions on the tempo of antigenic change.
To perform phylodynamic inference and goodness-of-fit analyses for complex phylodynamic models, we adopt a simulation-based approach that has become known as Approximate Bayesian Computation (ABC) [26]. Our first goal is to estimate the posterior density (1) of epidemiological and evolutionary model parameters under approximations to the likelihood of observed population incidence and phylogenetic data. The prior density can be used to incorporate existing information or limit the range of plausible values of model parameters. Our second goal is to assess fitted phylodynamic models based on a recent extension of ABC [27]. ABC methods circumvent computations of the likelihood by comparing the observed data to simulated data in terms of many, lower-dimensional summary statistics, , such as those in Figure 1. Using a distance function that compares summaries, each simulation is weighted according to the magnitude of the summary error under a weighting scheme, and this value is used in place of the likelihood term in Monte Carlo algorithms. In essence, ABC is a particular auxiliary variable Monte Carlo method, where the summary errors take on the role of auxiliary variables. Integrating these errors out, the ABC likelihood approximation adopted here is (2) where the weighting scheme is typically the Indicator (3) with tolerance parameter or the Exponential (4) with. Intuitively, the summary errors indicate how well a parameterized model reproduces the observed data. Once Monte Carlo algorithms such as the Markov Chain Monte Carlo (MCMC) sampler proposed by Marjoram et al. [34] have converged, the magnitude of the summary errors can be used to diagnose goodness-of-fit with respect to each of the summaries. To use this detailed information on each summary, we prefer using (2) to the Mahanalobis approximation (see [26]). Although uncommon, we typically use the log ratio so that the errors can be uniformly interpreted as fold-deviations. Parameter inference using ABC is approximate in that the ABC target density approaches the posterior density (1) as tends to zero if the summaries are sufficient for [26]. We use a Monte Carlo algorithm that is very similar to the MCMC sampler in Figure 2. A full specification of the algorithm is given in Text S1. It is typically difficult to establish the sufficiency of phylodynamic summaries analytically, and instead a small set of summaries is chosen such that model parameters of interest can be estimated [26]. Table 1 lists basic features of H3N2 epidemiological and phylogenetic data that were primarily considered in this study. Phylodynamic models were fitted and assessed against the features of the Dutch incidence data and the viral phylogeny derived under the Exponential clock model. The differences between these summaries and those derived from the remaining data in Figure 1 were used to set the ABC tolerances large enough so that inference is robust to the choice of phylogenetic reconstruction method and reporting country. Although smaller tolerances can be computationally feasible, these were not supported by the additional data considered. We typically use the Indicator weighting scheme (3) with tolerances that encompass differences in summary values across reporting countries and/or reconstruction methods, see Table 1. When a model never fits a particular summary well, we use (4) to give a mild prior preference to small errors [27]. See Text S1 for further details. A frequent problem in phylodynamic modeling is to determine if a specific model parameter should be included. For example, it can be unclear which types of ecological interactions between antigenic variants underlie pathogen phylodynamics, or if the residual selection parameter in (5c–5d) is required in addition to a given ecological mechanism that induces immune selection. Following existing variable selection procedures [28], we use an additional indicator variable to denote whether a single model parameter is present () or absent () and estimate its posterior probability under equation (2). Here, we use a standard spike-and-slab variable selection procedure [28] to estimate inclusion probabilities of the residual selection parameter.
To illustrate ABC methodology with the summaries in Table 1, we begin with a classical phenomenological model that implicitly accounts for antigenic drift through gradual loss of immunity [37]. H3N2 phylodynamics are represented with a spatial two-tier system of equations that is a special case of (5) when the antigenic emergence rate is set to (6) For simplicity, we will refer to (5) without antigenic variants as the SEIRS model. While several models have been able to simulate phylodynamics that are consistent with some aspects of the observed data, most notably the ladder-like phylogeny of H3N2' s haemagglutinin gene [4], [5], [41], none have been quantitatively fitted and tested against a set of epidemiological and molecular genetic features such as those in Figure 1. Here, we focus on the epochal evolution model as formulated in [7] within the above spatial framework, which is identical to (5) when antigenic variants are interpreted as major antigenic clusters. To fit (5) to the serial replacement of 11 major antigenic clusters within 1968–2002 [31], we define an antigenic cluster as any antigenic unit that survives for at least years and use the summaries in Table 1 as well as the number of antigenic clusters generated in 1968–2002 (nclust). Following [7], the emergence rate is set to increase with age, (7), and the scaling parameter is estimated. For simplicity, we refer to (5) with this antigenic emergence rate and an antigenic resolution that is determined by nclust as the epochal evolution model.
Fitting mechanistic models to infectious disease dynamics of RNA viruses that may escape immunity is notoriously difficult, and key epidemiological parameters such as can be estimated only under tacit assumptions from incidence time series [18], [19], [43]. Currently, alternative statistical synthesis approaches are explored to harness the information in complementary data sources [25], [44], [45]. Considering summaries of interpandemic H3N2 sequence and surveillance data, we show here that ABC can be used to fit and assess complex phylodynamic models which describe how evolutionary and ecological processes of the influenza virus may interact. Key phylodynamic parameters could be estimated under relatively weak assumptions (Table 2), and ABC diagnosed readily if and in which direction the two considered models deviate from all the available data taken together. Phylodynamic parameter inference and goodness-of-fit analyses rely critically on the possibility to combine epidemiological and molecular genetic data. In particular, H3N2 case report data were not sufficient to disentangle the reporting rate from epidemiological parameters, and measures of sequence divergence and diversity were not sufficient to separate the population size from evolutionary parameters. To the extent that other RNA viruses are characterized by different phylodynamic behavior, different sets of summaries must be identified in each case to replace likelihood calculations. ABC relates evolutionary and epidemiological data mechanistically through an evolving dynamic system and thereby allows us to investigate empirical phylodynamic hypotheses more directly than is possible with other statistical data synthesis approaches [44], [45]. Whenever the evolution and ecology of the virus are inseparably linked [1], case report and phylogenetic summaries are co-dependent. In general, this reduces the degrees of freedom of a phylodynamic model in reproducing features of both types of data simultaneously, and may reveal model inconsistencies. For example, the fitted epochal evolution model could not reproduce the TMRCA' s and the population attack rates at the same time (Figure 5L). The reported parameter estimates and summary errors are derived by conditioning only on the phylodynamic summaries and weighting schemes described in Table 1. ABC is sensitive to the chosen summary statistics and the tolerances since they determine how the prior is re-weighted in light of the presented evidence (see for example Table S3 and Figure S4 in Text S1) [26]. Here, we chose broad enough tolerances such that the weighting schemes are robust to differences in surveillance time series from the Netherlands, France and the US. This approach seems appropriate to avoid overfitting in the context of the limitations of syndromic influenza surveillance, but may be less suited in the analysis of other viral infectious diseases. Figures 3 and 5 illustrate that the resulting dimension reduction regularizes the underlying, intractable likelihood into a smooth, yet well-defined surface such that key phylodynamic model parameters are identifieable and goodness-of-fit can be characterized. It remains unclear to what degree the use of sufficient statistics or the full historical data would be desirable. Indeed, when infectious pathogens escape immunity, the likelihood surface can be especially complex [18], [43]. Likelihood-based inference is then sensitive to small changes in the complete historical data [46], [47], which can be problematic when the reported incidence time series or viral phylogeny is itself subject to considerable uncertainty and/or bias [2], [48]. We used ABC to fit mechanistic phylodynamic models of interpandemic influenza A (H3N2) to summaries of surveillance data from the Netherlands and sequence data from Northern Europe. Influenza is a globally circulating virus, and the mechanistic models considered must account for the replenishment of genetic variants from outside Northern Europe in order to reproduce features of influenza' s phylogeny. In contrast, semi- or non-parametric models of population dynamics that are used in coalescent methods do not necessarily require this layer of spatial complexity [9], [20]. Here, the mechanistic structure of Eqns. (5) constrains the set of possible phylodynamic trajectories in such a way that influenza' s global disease dynamics must be explicitly accounted for. Put more generally, the quantitative features of H3N2 sequence and incidence data contain sufficient information to determine at least some basic aspects of phylodynamic process models statistically. The two models we analyzed show clear limitations in their ability to replicate features of H3N2 sequence and surveillance data simultaneously, and the ABC error diagnostics give some indication how these models could be refined (Figures 3 and 5). For example, the phylogenies generated under the SEIRS and the epochal evolution models have, across time, more lineages than the observed HA phylogeny (Table 3). One possible explanation is that localized extinctions may not occur sufficiently often under the re-seeding source-sink framework, suggesting that models with more detailed population structure, either in space or by age, may result in thinner phylogenies. Accounting for these types of population structure can be critical for understanding viral phylodynamics; here we showed that the fit of the epochal evolution model to both sequence and incidence data depends critically on the assumed spatial model structure and the associated φ (see Text S1). The SEIRS model could not generate the irregularity in observed incidence data. In comparison, our analysis of the epochal evolution model demonstrates that epochal evolutionary processes can easily excite irregular between-season dynamics that match observed data (see Figures 5I and S18 in Text S1). Since the virus is known to be under intense immune selection [2], it seems plausible that antigenic evolution is an important co-factor in explaining influenza' s irregular seasonality in temperate regions [49]. Several alternative models have been proposed to reproduce H3N2' s narrow HA phylogeny. Here, we identified an additional, testable constraint for these models on surveillance data, that arises through the phylodynamic interactions in Eqns. (5). The cluster-specific duration of immunity must be sufficiently long to avoid deep phylogenetic branching. If the fitted values of and are correlated, this in turn implies a characteristic range of population level attack rates that can be tested against available data as in Figure 5L. In particular, the duration of immunity can be lower if the time between replacement events is shorter. Thus, while H3N2' s limited standing genetic diversity provides information on the strength of immune interactions between H3N2 antigenic variants, this second constraint may help identify the tempo of antigenic evolution. For the epochal evolution model with source-sink migration dynamics, the average simulated waiting time is years from the emergence of the current antigenic cluster to the next successfully invading offspring antigenic cluster, and this implies phylodynamics that are inconsistent with the molecular genetic and epidemiological summaries in Table 1 taken together. More frequent and more gradual transitions between antigenic variants that are smaller than H3N2' s antigenic clusters would allow for lower estimates of that are more in line with observed population level attack rates, break weaker refractory oscillations in their onset, and might also provide sufficient, continual selection pressures to explain the fast divergence in H3N2' s HA phylogeny [8]. In this case, sequence and surveillance data would point to a finer antigenic resolution than the one suggested through antigenic map analyses [31]. Alternatively, it is also possible that finer population structure, either in space or by age, could increase extinction rates and thereby allow for a narrow HA phylogeny under a broader, more realistic set of epidemiological parameters without accelerating the tempo of antigenic evolution per se. More broadly, both types of data are now increasingly becoming available for RNA viruses [9]–[13]. This study indicates that these data, when considered simultaneously, may drastically constrain parameter space and readily expose model deficiencies, so that ABC appears as a well-suited tool to explore the phylodynamics of RNA viruses.
The infectious disease dynamics of many viral pathogens like influenza, norovirus and coronavirus are inextricably tied to their evolution. This interaction between evolutionary and ecological processes complicates our ability to understand the infectious disease behavior of rapidly evolving pathogens. Most statistical methods for the analysis of these “phylodynamics” require that the likelihood of the data can be explicitly calculated. Currently, this is not possible for many phylodynamic models, so that questions on the interaction between viral variants cannot be well-addressed within this framework. Simulation-based statistical methods circumvent likelihood calculations. Considering interpandemic human influenza A virus subtype H3N2, we here illustrate the effectiveness of these methods to fit and assess complex phylodynamic models against both sequence and surveillance data. We find that combining molecular genetic and epidemiological data is key to estimate phylodynamic parameters reliably. Moreover, the information in the available data taken together is enough to expose quantitative model inconsistencies. Methods such as ABC which can combine sequence and surveillance data appear to be well-suited to fit and assess mechanistic hypotheses on the phylodynamics of RNA viruses.
public health and epidemiology
evolutionary ecology
mathematics
ecology
evolutionary biology
epidemiology
statistics
theoretical ecology
biology
population biology
statistical methods
Phylodynamic Inference and Model Assessment with Approximate Bayesian Computation: Influenza as a Case Study