Risk averse reproduction numbers improve resurgence detection.

The effective reproduction number R is a prominent statistic for inferring the transmissibility of infectious diseases and effectiveness of interventions. R purportedly provides an easy-to-interpret threshold for deducing whether an epidemic will grow (R>1) or decline (R<1). We posit that this interpretation can be misleading and statistically overconfident when applied to infections accumulated from groups featuring heterogeneous dynamics. These groups may be delineated by geography, infectiousness or sociodemographic factors. In these settings, R implicitly weights the dynamics of the groups by their number of circulating infections. We find that this weighting can cause delayed detection of outbreak resurgence and premature signalling of epidemic control because it underrepresents the risks from highly transmissible groups. Applying E-optimal experimental design theory, we develop a weighting algorithm to minimise these issues, yielding the risk averse reproduction number E. Using simulations, analytic approaches and real-world COVID-19 data stratified at the city and district level, we show that E meaningfully summarises transmission dynamics across groups, balancing bias from the averaging underlying R with variance from directly using local group estimates. An E>1generates timely resurgence signals (upweighting risky groups), while an E<1ensures local outbreaks are under control. We propose E as an alternative to R for informing policy and assessing transmissibility at large scales (e.g., state-wide or nationally), where R is commonly computed but well-mixed or homogeneity assumptions break down.

Angular reproduction numbers improve estimates of transmissibility when disease generation times are misspecified or time-varying.

We introduce the angular reproduction number Ω, which measures time-varying changes in epidemic transmissibility resulting from variations in both the effective reproduction number R, and generation time distribution w. Predominant approaches for tracking pathogen spread infer either R or the epidemic growth rate r. However, R is biased by mismatches between the assumed and true w, while r is difficult to interpret in terms of the individual-level branching process underpinning transmission. R and r may also disagree on the relative transmissibility of epidemics or variants (i.e. rA > rB does not imply RA > RB for variants A and B). We find that Ω responds meaningfully to mismatches and time-variations in w while mostly maintaining the interpretability of R. We prove that Ω > 1 implies R > 1 and that Ω agrees with r on the relative transmissibility of pathogens. Estimating Ω is no more difficult than inferring R, uses existing software, and requires no generation time measurements. These advantages come at the expense of selecting one free parameter. We propose Ω as complementary statistic to R and r that improves transmissibility estimates when w is misspecified or time-varying and better reflects the impact of interventions, when those interventions concurrently change R and w or alter the relative risk of co-circulating pathogens.

Accounting for the Potential of Overdispersion in Estimation of the Time-varying Reproduction Number.

BACKGROUND: The time-varying reproduction number, Rt, is commonly used to monitor the transmissibility of an infectious disease during an epidemic, but standard methods for estimating Rt seldom account for the impact of overdispersion on transmission. METHODS: We developed a negative binomial framework to estimate Rt and a time-varying dispersion parameter (kt). We applied the framework to COVID-19 incidence data in Hong Kong in 2020 and 2021. We conducted a simulation study to compare the performance of our model with the conventional Poisson-based approach. RESULTS: Our framework estimated an Rt peaking around 4 (95% credible interval = 3.13, 4.30), similar to that from the Poisson approach but with a better model fit. Our approach further estimated kt <0.5 at the start of both waves, indicating appreciable heterogeneity in transmission. We also found that kt decreased sharply to around 0.4 when a large cluster of infections occurred. CONCLUSIONS: Our proposed approach can contribute to the estimation of Rt and monitoring of the time-varying dispersion parameters to quantify the role of superspreading.

Seasonal dynamics of the wild rodent faecal virome.

Viral discovery studies in wild animals often rely on cross-sectional surveys at a single time point. As a result, our understanding of the temporal stability of wild animal viromes remains poorly resolved. While studies of single host-virus systems indicate that host and environmental factors influence seasonal virus transmission dynamics, comparable insights for whole viral communities in multiple hosts are lacking. Utilizing noninvasive faecal samples from a long-term wild rodent study, we characterized viral communities of three common European rodent species (Apodemus sylvaticus, A. flavicollis and Myodes glareolus) living in temperate woodland over a single year. Our findings indicate that a substantial fraction of the rodent virome is seasonally transient and associated with vertebrate or bacteria hosts. Further analyses of one of the most common virus families, Picornaviridae, show pronounced temporal changes in viral richness and evenness, which were associated with concurrent and up to ~3-month lags in host density, ambient temperature, rainfall and humidity, suggesting complex feedbacks from the host and environmental factors on virus transmission and shedding in seasonal habitats. Overall, this study emphasizes the importance of understanding the seasonal dynamics of wild animal viromes in order to better predict and mitigate zoonotic risks.

A Bayesian nonparametric method for detecting rapid changes in disease transmission.

Whether an outbreak of infectious disease is likely to grow or dissipate is determined through the time-varying reproduction number, Rt. Real-time or retrospective identification of changes in Rt following the imposition or relaxation of interventions can thus contribute important evidence about disease transmission dynamics which can inform policymaking. Here, we present a method for estimating shifts in Rt within a renewal model framework. Our method, which we call EpiCluster, is a Bayesian nonparametric model based on the Pitman-Yor process. We assume that Rt is piecewise-constant, and the incidence data and priors determine when or whether Rt should change and how many times it should do so throughout the series. We also introduce a prior which induces sparsity over the number of changepoints. Being Bayesian, our approach yields a measure of uncertainty in Rt and its changepoints. EpiCluster is fast, straightforward to use, and we demonstrate that it provides automated detection of rapid changes in transmission, either in real-time or retrospectively, for synthetic data series where the Rt profile is known. We illustrate the practical utility of our method by fitting it to case data of outbreaks of COVID-19 in Australia and Hong Kong, where it finds changepoints coinciding with the imposition of non-pharmaceutical interventions. Bayesian nonparametric methods, such as ours, allow the volume and complexity of the data to dictate the number of parameters required to approximate the process and should find wide application in epidemiology. This manuscript was submitted as part of a theme issue on "Modelling COVID-19 and Preparedness for Future Pandemics".

Fundamental limits on inferring epidemic resurgence in real time using effective reproduction numbers.

We find that epidemic resurgence, defined as an upswing in the effective reproduction number (R) of the contagion from subcritical to supercritical values, is fundamentally difficult to detect in real time. Inherent latencies in pathogen transmission, coupled with smaller and intrinsically noisier case incidence across periods of subcritical spread, mean that resurgence cannot be reliably detected without significant delays of the order of the generation time of the disease, even when case reporting is perfect. In contrast, epidemic suppression (where R falls from supercritical to subcritical values) may be ascertained 5-10 times faster due to the naturally larger incidence at which control actions are generally applied. We prove that these innate limits on detecting resurgence only worsen when spatial or demographic heterogeneities are incorporated. Consequently, we argue that resurgence is more effectively handled proactively, potentially at the expense of false alarms. Timely responses to recrudescent infections or emerging variants of concern are more likely to be possible when policy is informed by a greater quality and diversity of surveillance data than by further optimisation of the statistical models used to process routine outbreak data.

Genomic Epidemiology of Early SARS-CoV-2 Transmission Dynamics, Gujarat, India.

Limited genomic sampling in many high-incidence countries has impeded studies of severe respiratory syndrome coronavirus 2 (SARS-CoV-2) genomic epidemiology. Consequently, critical questions remain about the generation and global distribution of virus genetic diversity. We investigated SARS-CoV-2 transmission dynamics in Gujarat, India, during the state's first epidemic wave to shed light on spread of the virus in one of the regions hardest hit by the pandemic. By integrating case data and 434 whole-genome sequences sampled across 20 districts, we reconstructed the epidemic dynamics and spatial spread of SARS-CoV-2 in Gujarat. Our findings indicate global and regional connectivity and population density were major drivers of the Gujarat outbreak. We detected >100 virus lineage introductions, most of which appear to be associated with international travel. Within Gujarat, virus dissemination occurred predominantly from densely populated regions to geographically proximate locations that had low population density, suggesting that urban centers contributed disproportionately to virus spread.

Are Skyline Plot-Based Demographic Estimates Overly Dependent on Smoothing Prior Assumptions?

In Bayesian phylogenetics, the coalescent process provides an informative framework for inferring changes in the effective size of a population from a phylogeny (or tree) of sequences sampled from that population. Popular coalescent inference approaches such as the Bayesian Skyline Plot, Skyride, and Skygrid all model these population size changes with a discontinuous, piecewise-constant function but then apply a smoothing prior to ensure that their posterior population size estimates transition gradually with time. These prior distributions implicitly encode extra population size information that is not available from the observed coalescent data or tree. Here, we present a novel statistic, $\Omega$, to quantify and disaggregate the relative contributions of the coalescent data and prior assumptions to the resulting posterior estimate precision. Our statistic also measures the additional mutual information introduced by such priors. Using $\Omega$ we show that, because it is surprisingly easy to overparametrize piecewise-constant population models, common smoothing priors can lead to overconfident and potentially misleading inference, even under robust experimental designs. We propose $\Omega$ as a useful tool for detecting when effective population size estimates are overly reliant on prior assumptions and for improving quantification of the uncertainty in those estimates.[Coalescent processes; effective population size; information theory; phylodynamics; prior assumptions; skyline plots.].

Improved estimation of time-varying reproduction numbers at low case incidence and between epidemic waves.

We construct a recursive Bayesian smoother, termed EpiFilter, for estimating the effective reproduction number, R, from the incidence of an infectious disease in real time and retrospectively. Our approach borrows from Kalman filtering theory, is quick and easy to compute, generalisable, deterministic and unlike many current methods, requires no change-point or window size assumptions. We model R as a flexible, hidden Markov state process and exactly solve forward-backward algorithms, to derive R estimates that incorporate all available incidence information. This unifies and extends two popular methods, EpiEstim, which considers past incidence, and the Wallinga-Teunis method, which looks forward in time. We find that this combination of maximising information and minimising assumptions significantly reduces the bias and variance of R estimates. Moreover, these properties make EpiFilter more statistically robust in periods of low incidence, where several existing methods can become destabilised. As a result, EpiFilter offers improved inference of time-varying transmission patterns that are advantageous for assessing the risk of upcoming waves of infection or the influence of interventions, in real time and at various spatial scales.

Jointly Inferring the Dynamics of Population Size and Sampling Intensity from Molecular Sequences.

Estimating past population dynamics from molecular sequences that have been sampled longitudinally through time is an important problem in infectious disease epidemiology, molecular ecology, and macroevolution. Popular solutions, such as the skyline and skygrid methods, infer past effective population sizes from the coalescent event times of phylogenies reconstructed from sampled sequences but assume that sequence sampling times are uninformative about population size changes. Recent work has started to question this assumption by exploring how sampling time information can aid coalescent inference. Here, we develop, investigate, and implement a new skyline method, termed the epoch sampling skyline plot (ESP), to jointly estimate the dynamics of population size and sampling rate through time. The ESP is inspired by real-world data collection practices and comprises a flexible model in which the sequence sampling rate is proportional to the population size within an epoch but can change discontinuously between epochs. We show that the ESP is accurate under several realistic sampling protocols and we prove analytically that it can at least double the best precision achievable by standard approaches. We generalize the ESP to incorporate phylogenetic uncertainty in a new Bayesian package (BESP) in BEAST2. We re-examine two well-studied empirical data sets from virus epidemiology and molecular evolution and find that the BESP improves upon previous coalescent estimators and generates new, biologically useful insights into the sampling protocols underpinning these data sets. Sequence sampling times provide a rich source of information for coalescent inference that will become increasingly important as sequence collection intensifies and becomes more formalized.

Adaptive Estimation for Epidemic Renewal and Phylogenetic Skyline Models.

Estimating temporal changes in a target population from phylogenetic or count data is an important problem in ecology and epidemiology. Reliable estimates can provide key insights into the climatic and biological drivers influencing the diversity or structure of that population and evidence hypotheses concerning its future growth or decline. In infectious disease applications, the individuals infected across an epidemic form the target population. The renewal model estimates the effective reproduction number, R, of the epidemic from counts of observed incident cases. The skyline model infers the effective population size, N, underlying a phylogeny of sequences sampled from that epidemic. Practically, R measures ongoing epidemic growth while N informs on historical caseload. While both models solve distinct problems, the reliability of their estimates depends on p-dimensional piecewise-constant functions. If p is misspecified, the model might underfit significant changes or overfit noise and promote a spurious understanding of the epidemic, which might misguide intervention policies or misinform forecasts. Surprisingly, no transparent yet principled approach for optimizing p exists. Usually, p is heuristically set, or obscurely controlled via complex algorithms. We present a computable and interpretable p-selection method based on the minimum description length (MDL) formalism of information theory. Unlike many standard model selection techniques, MDL accounts for the additional statistical complexity induced by how parameters interact. As a result, our method optimizes p so that R and N estimates properly and meaningfully adapt to available data. It also outperforms comparable Akaike and Bayesian information criteria on several classification problems, given minimal knowledge of the parameter space, and exposes statistical similarities among renewal, skyline, and other models in biology. Rigorous and interpretable model selection is necessary if trustworthy and justifiable conclusions are to be drawn from piecewise models. [Coalescent processes; epidemiology; information theory; model selection; phylodynamics; renewal models; skyline plots].

Using information theory to optimise epidemic models for real-time prediction and estimation.

The effective reproduction number, Rt, is a key time-varying prognostic for the growth rate of any infectious disease epidemic. Significant changes in Rt can forewarn about new transmissions within a population or predict the efficacy of interventions. Inferring Rt reliably and in real-time from observed time-series of infected (demographic) data is an important problem in population dynamics. The renewal or branching process model is a popular solution that has been applied to Ebola and Zika virus disease outbreaks, among others, and is currently being used to investigate the ongoing COVID-19 pandemic. This model estimates Rt using a heuristically chosen piecewise function. While this facilitates real-time detection of statistically significant Rt changes, inference is highly sensitive to the function choice. Improperly chosen piecewise models might ignore meaningful changes or over-interpret noise-induced ones, yet produce visually reasonable estimates. No principled piecewise selection scheme exists. We develop a practical yet rigorous scheme using the accumulated prediction error (APE) metric from information theory, which deems the model capable of describing the observed data using the fewest bits as most justified. We derive exact posterior prediction distributions for infected population size and integrate these within an APE framework to obtain an exact and reliable method for identifying the piecewise function best supported by available epidemic data. We find that this choice optimises short-term prediction accuracy and can rapidly detect salient fluctuations in Rt, and hence the infected population growth rate, in real-time over the course of an unfolding epidemic. Moreover, we emphasise the need for formal selection by exposing how common heuristic choices, which seem sensible, can be misleading. Our APE-based method is easily computed and broadly applicable to statistically similar models found in phylogenetics and macroevolution, for example. Our results explore the relationships among estimate precision, forecast reliability and model complexity.

An exact method for quantifying the reliability of end-of-epidemic declarations in real time.

We derive and validate a novel and analytic method for estimating the probability that an epidemic has been eliminated (i.e. that no future local cases will emerge) in real time. When this probability crosses 0.95 an outbreak can be declared over with 95% confidence. Our method is easy to compute, only requires knowledge of the incidence curve and the serial interval distribution, and evaluates the statistical lifetime of the outbreak of interest. Using this approach, we show how the time-varying under-reporting of infected cases will artificially inflate the inferred probability of elimination, leading to premature (false-positive) end-of-epidemic declarations. Contrastingly, we prove that incorrectly identifying imported cases as local will deceptively decrease this probability, resulting in delayed (false-negative) declarations. Failing to sustain intensive surveillance during the later phases of an epidemic can therefore substantially mislead policymakers on when it is safe to remove travel bans or relax quarantine and social distancing advisories. World Health Organisation guidelines recommend fixed (though disease-specific) waiting times for end-of-epidemic declarations that cannot accommodate these variations. Consequently, there is an unequivocal need for more active and specialised metrics for reliably identifying the conclusion of an epidemic.

Key questions for modelling COVID-19 exit strategies.

Combinations of intense non-pharmaceutical interventions (lockdowns) were introduced worldwide to reduce SARS-CoV-2 transmission. Many governments have begun to implement exit strategies that relax restrictions while attempting to control the risk of a surge in cases. Mathematical modelling has played a central role in guiding interventions, but the challenge of designing optimal exit strategies in the face of ongoing transmission is unprecedented. Here, we report discussions from the Isaac Newton Institute 'Models for an exit strategy' workshop (11-15 May 2020). A diverse community of modellers who are providing evidence to governments worldwide were asked to identify the main questions that, if answered, would allow for more accurate predictions of the effects of different exit strategies. Based on these questions, we propose a roadmap to facilitate the development of reliable models to guide exit strategies. This roadmap requires a global collaborative effort from the scientific community and policymakers, and has three parts: (i) improve estimation of key epidemiological parameters; (ii) understand sources of heterogeneity in populations; and (iii) focus on requirements for data collection, particularly in low-to-middle-income countries. This will provide important information for planning exit strategies that balance socio-economic benefits with public health.

Robust Design for Coalescent Model Inference.

The coalescent process describes how changes in the size or structure of a population influence the genealogical patterns of sequences sampled from that population. The estimation of (effective) population size changes from genealogies that are reconstructed from these sampled sequences is an important problem in many biological fields. Often, population size is characterized by a piecewise-constant function, with each piece serving as a population size parameter to be estimated. Estimation quality depends on both the statistical coalescent inference method employed, and on the experimental protocol, which controls variables such as the sampling of sequences through time and space, or the transformation of model parameters. While there is an extensive literature on coalescent inference methodology, there is comparatively little work on experimental design. The research that does exist is largely simulation-based, precluding the development of provable or general design theorems. We examine three key design problems: temporal sampling of sequences under the skyline demographic coalescent model, spatio-temporal sampling under the structured coalescent model, and time discretization for sequentially Markovian coalescent models. In all cases, we prove that 1) working in the logarithm of the parameters to be inferred (e.g., population size) and 2) distributing informative coalescent events uniformly among these log-parameters, is uniquely robust. "Robust" means that the total and maximum uncertainty of our parameter estimates are minimized, and made insensitive to their unknown (true) values. This robust design theorem provides rigorous justification for several existing coalescent experimental design decisions and leads to usable guidelines for future empirical or simulation-based investigations. Given its persistence among models, this theorem may form the basis of an experimental design paradigm for coalescent inference.

Exact Bayesian inference for phylogenetic birth-death models.

Motivation: Inferring the rates of change of a population from a reconstructed phylogeny of genetic sequences is a central problem in macro-evolutionary biology, epidemiology and many other disciplines. A popular solution involves estimating the parameters of a birth-death process (BDP), which links the shape of the phylogeny to its birth and death rates. Modern BDP estimators rely on random Markov chain Monte Carlo (MCMC) sampling to infer these rates. Such methods, while powerful and scalable, cannot be guaranteed to converge, leading to results that may be hard to replicate or difficult to validate. Results: We present a conceptually and computationally different parametric BDP inference approach using flexible and easy to implement Snyder filter (SF) algorithms. This method is deterministic so its results are provable, guaranteed and reproducible. We validate the SF on constant rate BDPs and find that it solves BDP likelihoods known to produce robust estimates. We then examine more complex BDPs with time-varying rates. Our estimates compare well with a recently developed parametric MCMC inference method. Lastly, we perform model selection on an empirical Agamid species phylogeny, obtaining results consistent with the literature. The SF makes no approximations, beyond those required for parameter quantization and numerical integration and directly computes the posterior distribution of model parameters. It is a promising alternative inference algorithm that may serve either as a standalone Bayesian estimator or as a useful diagnostic reference for validating more involved MCMC strategies. Availability and implementation: The Snyder filter is implemented in Matlab and the time-varying BDP models are simulated in R. The source code and data are freely available at https://github.com/kpzoo/snyder-birth-death-code. Supplementary information: Supplementary data are available at Bioinformatics online.

On signalling and estimation limits for molecular birth-processes.

Understanding and uncovering the mechanisms or motifs that molecular networks employ to regulate noise is a key problem in cell biology. As it is often difficult to obtain direct and detailed insight into these mechanisms, many studies instead focus on assessing the best precision attainable on the signalling pathways that compose these networks. Molecules signal one another over such pathways to solve noise regulating estimation and control problems. Quantifying the maximum precision of these solutions delimits what is achievable and allows hypotheses about underlying motifs to be tested without requiring detailed biological knowledge. The pathway capacity, which defines the maximum rate of transmitting information along it, is a widely used proxy for precision. Here it is shown, for estimation problems involving elementary yet biologically relevant birth-process networks, that capacity can be surprisingly misleading. A time-optimal signalling motif, called birth-following, is derived and proven to better the precision expected from the capacity, provided the maximum signalling rate constraint is large and the mean one above a certain threshold. When the maximum constraint is relaxed, perfect estimation is predicted by the capacity. However, the true achievable precision is found highly variable and sensitive to the mean constraint. Since the same capacity can map to different combinations of rate constraints, it can only equivocally measure precision. Deciphering the rate constraints on a signalling pathway may therefore be more important than computing its capacity.

Point process analysis of noise in early invertebrate vision.

Noise is a prevalent and sometimes even dominant aspect of many biological processes. While many natural systems have adapted to attenuate or even usefully integrate noise, the variability it introduces often still delimits the achievable precision across biological functions. This is particularly so for visual phototransduction, the process responsible for converting photons of light into usable electrical signals (quantum bumps). Here, randomness of both the photon inputs (regarded as extrinsic noise) and the conversion process (intrinsic noise) are seen as two distinct, independent and significant limitations on visual reliability. Past research has attempted to quantify the relative effects of these noise sources by using approximate methods that do not fully account for the discrete, point process and time ordered nature of the problem. As a result the conclusions drawn from these different approaches have led to inconsistent expositions of phototransduction noise performance. This paper provides a fresh and complete analysis of the relative impact of intrinsic and extrinsic noise in invertebrate phototransduction using minimum mean squared error reconstruction techniques based on Bayesian point process (Snyder) filters. An integrate-fire based algorithm is developed to reliably estimate photon times from quantum bumps and Snyder filters are then used to causally estimate random light intensities both at the front and back end of the phototransduction cascade. Comparison of these estimates reveals that the dominant noise source transitions from extrinsic to intrinsic as light intensity increases. By extending the filtering techniques to account for delays, it is further found that among the intrinsic noise components, which include bump latency (mean delay and jitter) and shape (amplitude and width) variance, it is the mean delay that is critical to noise performance. As the timeliness of visual information is important for real-time action, this delay could potentially limit the speed at which invertebrates can respond to stimuli. Consequently, if one wants to increase visual fidelity, reducing the photoconversion lag is much more important than improving the regularity of the electrical signal.

Optimal point process filtering and estimation of the coalescent process.

The coalescent process is a widely used approach for inferring the demographic history of a population, from samples of its genetic diversity. Several parametric and non-parametric coalescent inference methods, involving Markov chain Monte Carlo, Gaussian processes, and other algorithms, already exist. However, these techniques are not always easy to adapt and apply, thus creating a need for alternative methodologies. We introduce the Bayesian Snyder filter as an easily implementable and flexible minimum mean square error estimator for parametric demographic functions on fixed genealogies. By reinterpreting the coalescent as a self-exciting Markov process, we show that the Snyder filter can be applied to both isochronously and heterochronously sampled datasets. We analytically solve the filter equations for the constant population size Kingman coalescent, derive expressions for its mean squared estimation error, and estimate its robustness to prior distribution specification. For populations with deterministically time-varying size we numerically solve the Snyder equations, and test this solution on common demographic models. We find that the Snyder filter accurately recovers the true demographic history for these models. We also apply the filter to a well-studied, dataset of hepatitis C virus sequences and show that the filter compares well to a popular phylodynamic inference method. The Snyder filter is an exact (given discretised priors, it does not approximate the posterior) and direct Bayesian estimation method that has the potential to become a useful alternative tool for coalescent inference.