The Spline Clock: When Molecular Evolution Doesn't Keep Perfect Time

On March 11, 2020, the World Health Organization declared COVID-19 a pandemic. By that point, researchers had already been sequencing SARS-CoV-2 genomes for months, watching the virus spread from China across the globe. But when they tried to figure out how fast the virus was evolving — a critical piece of information for understanding its pandemic potential — they ran into a puzzling problem. The rate they calculated depended entirely on when they started looking. Sequences collected over weeks told one evolutionary story; sequences collected over months told another. This time-dependent rate phenomenon had been documented before, in foamy viruses and lentiviruses and other pathogens that accumulate genetic changes rapidly. But standard phylogenetic methods assumed the evolutionary rate was constant, like a metronome clicking at the same tempo throughout a virus's history. It wasn't. And that mismatch was distorting everything scientists thought they knew about how these viruses evolved.

A new modeling framework published in 2026 by Datta, Lemey, and Suchard tackles this problem directly. Their "spline clock model" treats evolutionary rates as what they actually are: smooth, time-varying functions that can rise and fall as populations expand, immune pressures mount, or transmission dynamics shift. The approach applies sophisticated mathematical machinery — cubic B-splines, Gaussian Markov random fields, Gauss-Legendre quadrature — to reconstruct how fast viruses have been evolving at every point in time, not just an averaged constant that erases all the interesting variation. In simulation studies, the spline clock model recovered true time-varying rates more accurately and with tighter uncertainty bounds than competing approaches. Applied to real viral data, it found strong time-varying signals in both the long-term evolution of foamy virus over 100 million years and the spatial diffusion of SARS-CoV-2 across Europe during 2020. The work isn't just a methodological advance; it's a recognition that the molecular clock, long treated as a steady tick, sometimes behaves more like a heart rate — accelerating and decelerating in response to physiological demands.

The Science

Reconstructing History from Sequence Data

When biologists want to understand how a pathogen evolved — when it jumped to humans, which lineages spread where, how quickly its genome changed — they typically turn to phylogenetic trees. These branching diagrams represent the evolutionary history shared by a collection of genetic sequences, with the branching structure showing how different lineages are related and the branch lengths reflecting the accumulation of genetic substitutions over time. A phylogenetic tree for SARS-CoV-2, for instance, would show the relationships among thousands of viral genomes sampled from patients around the world, with branches representing transmission chains and substitutions accumulating as the virus replicated in new hosts.

Calibrating these trees — converting branch lengths measured in expected genetic changes into calendar time — is one of the central challenges in phylogenetics. One approach relies on fossil records to date ancient divergence events; another uses the sampling dates of sequences collected throughout an outbreak as direct time markers. The latter approach has become standard for RNA viruses like SARS-CoV-2, where genomes can be sequenced rapidly and repeatedly throughout an epidemic.

The conversion from genetic distance to time typically proceeds through a molecular clock model, which assigns an evolutionary rate to each branch. The rate is defined as the expected number of substitutions per site per unit time — a measure of how fast the genome is changing. The simplest clock model, known as the strict clock, assumes all branches share a single common rate, like a metronome keeping perfect time. More flexible "relaxed clock" models allow rates to vary across branches, but typically treat this variation as random and unstructured — a branch here might happen to evolve faster, another slower, but there's no systematic pattern over time.

Datta and colleagues identified a fundamental limitation in these approaches: they can't accommodate genuine temporal variation in evolutionary rates. When a virus first jumps into a new host population, it may experience explosive growth and rapid evolution. As herd immunity builds or interventions take effect, transmission dynamics change, and so might the selective pressures shaping the genome. Rates that are genuinely time-dependent — that rise and fall systematically over the course of an epidemic or over millions of years of viral diversification — can't be captured by models that assume constant rates or unstructured branch-to-branch variation.

The Time-Dependent Rate Problem

The phenomenon isn't new. Multiple studies have documented that evolutionary rate estimates depend on the timescale of sampling — a result with profound implications for how we interpret phylogenetic trees. If you estimate the evolutionary rate of a virus using sequences collected over a few weeks, you get one answer. Use sequences spanning months, and you get another. Use sequences spanning years, and you get yet another. This isn't measurement error; it reflects genuine biology. Rates really do change over time.

The classic case study is foamy virus, a genus of retroviruses that infect various mammals including primates. Foamy viruses are ancient — they've been co-evolving with their hosts for tens of millions of years. But when researchers compare short-term and long-term rate estimates, they find dramatic discrepancies. Short-term estimates, based on contemporary sequences, suggest relatively fast evolution. Long-term estimates, calibrated against host divergence times, suggest much slower evolution. Aiewsakun and Katzourakis (2015) demonstrated that a simple power-law decay model could reconcile these estimates: the rate appears to decrease systematically as you look further back in time, following a predictable mathematical relationship.

Previous work by Datta and colleagues had developed a "polyepoch clock model" that accommodated time-varying rates by dividing time into discrete epochs and allowing the rate to vary across epochs. But this discretization imposed artificial discontinuities — sudden jumps where one epoch ended and another began. In regions of dense taxon sampling, these discontinuities could create spurious oscillations in the estimated rate trajectory, mixing noise with genuine signal.

The Spline Clock Solution

The new framework takes a different approach. Rather than discretizing time into epochs, it models the rate as a smooth, continuous function of time using cubic B-splines. B-splines are piecewise polynomial basis functions — flexible mathematical objects that can approximate almost any smooth curve while maintaining continuity and differentiability at the knot points where pieces join. By representing the log-rate as a cubic B-spline function, the model ensures that the estimated rate trajectory varies continuously and smoothly across the entire temporal domain, without artificial jumps or oscillations.

The mathematical structure requires some care. The likelihood function — the probability of observing the sequence data given the model parameters — involves integrals of the rate function over each branch of the phylogeny. Under the assumption that the time-varying rate scales all elements of the substitution matrix equally, these integrals determine the transition probabilities that connect parent and child nodes in the tree. For a piecewise constant rate (the polyepoch approach), these integrals have closed-form expressions. For a smooth spline function, they don't — but they can be approximated efficiently using Gauss-Legendre quadrature, a numerical integration method that selects evaluation points and weights optimally to maximize accuracy.

The researchers exploited a key property of the spline parameterization: between any two adjacent knots, the cubic spline is a polynomial of degree 3. The exponential of a cubic polynomial is smooth and well-approximated by higher-degree polynomials, making it an ideal candidate for Gauss-Legendre quadrature. By applying the quadrature rule locally on each knot interval, they achieved high accuracy with a fixed, predictable computational cost — five function evaluations per overlapping interval, regardless of how wiggly the true rate trajectory might be.

To prevent overfitting, the model imposes a Gaussian Markov random field (GMRF) prior on the spline coefficients. This prior penalizes roughness by shrinking differences between adjacent coefficients toward zero — a Bayesian version of the smoothing approaches commonly used in nonparametric regression. The amount of smoothing is controlled by a precision parameter that receives its own prior and is estimated from the data, allowing the model to adapt to the strength of the temporal signal without requiring researchers to make ad hoc choices.

The inference proceeds in a fully Bayesian framework, sampling from the joint posterior distribution of all parameters — tree topology, branch lengths, substitution model parameters, and spline coefficients — using Markov chain Monte Carlo (MCMC) methods. The implementation builds on the BEAST X platform, a widely used software package for Bayesian phylogenetic analysis, making the method accessible to researchers already familiar with standard phylogenetic workflows.

What They Found

Simulation Study: Accuracy Under Controlled Conditions

To evaluate the spline clock model's performance, the researchers first conducted a simulation study under controlled conditions where the true rate trajectory was known. They simulated a phylogeny of 100 tip nodes from an exponential growth coalescent model, then simulated sequences along the tree under a time-dependent rate that initially increases log-linearly before declining sharply — a two-phase pattern designed to test whether the model could capture qualitatively different behaviors in different time periods.

The simulation revealed stark differences between the three clock models tested. The uncorrelated relaxed clock model, which allows each branch to have its own rate without any temporal structure, failed to recover the true rate trajectory. Its posterior median initially reflected the increasing trend in the early phase, but then modulated sharply and never recovered the subsequent decline. The model lacked the structural assumption needed to accommodate systematic temporal variation.

The polyepoch clock model performed better but showed the roughness problems its developers had anticipated. While the posterior median broadly followed the true rate, the estimate exhibited noticeable oscillations — the signature of discontinuities introduced by discretizing time into epochs. In regions where many branches passed through the same time point, these oscillations were particularly pronounced.

The spline clock model, by contrast, closely tracked the true rate trajectory throughout the entire time range. The posterior median ran tight along the true curve, and the 95% Bayesian credible intervals were substantially narrower than those from either competing approach. The model had recovered the true dynamics accurately, with appropriate uncertainty quantification.

Simulation: Rate Recovery Accuracy by Model

Label	Value
Spline Clock (SCM)	95 estimated rate accuracy
Polyepoch Clock (PCM)	70 estimated rate accuracy
Uncorrelated Relaxed (UCLD)	40 estimated rate accuracy

True rate shown as dashed line. SCM demonstrates superior accuracy in recovering both increasing and decreasing phases of the true rate function.

To assess robustness, the researchers conducted additional simulations under constant rates and log-linearly increasing rates. Under constant rates, all three models produced reasonable estimates, as expected — this is the setting standard methods were designed for. Under log-linear increase, both the polyepoch and spline clock models recovered the true trajectory, with the spline clock again showing tighter credible intervals. The spline model performed particularly well even when the true rate trajectory was more complex, demonstrating that the smoothness constraint doesn't inappropriately smooth away genuine variation.

The researchers also investigated the choice of link function — the mathematical transformation that ensures the rate remains positive. They compared the exponential link (which treats the spline coefficients as coefficients of the log-rate) against a squared link (which treats them as coefficients of the rate itself). While both can work in principle, the exponential link produced better mixing in MCMC chains and avoided identifiability problems that arose from the non-monotonicity of the squaring operation. The choice matters in practice, not just in theory.

Foamy Virus: 100 Million Years of Evolutionary Dynamics

Having demonstrated the method's performance under simulation, the researchers applied it to real viral data. The first application examined the evolutionary dynamics of foamy virus (SFV) over the past 100 million years — a timescale spanning the diversification of primates and the long coevolutionary history between foamy viruses and their hosts.

The analysis revealed strong time-varying signal in the evolutionary rate. The rate appeared to decrease systematically over the deep time horizon, consistent with the time-dependent rate phenomenon documented in previous studies but now estimated with unprecedented smoothness and precision. The spline representation allowed the researchers to trace the rate trajectory continuously, without the artifacts that epoch-based approaches might introduce over such a vast timescale.

The inferred maximum clade credibility (MCC) tree — a summary tree that captures the highest posterior probability clades — showed the branching patterns of foamy virus lineages across primates, with uncertainty quantified at each internal node. The temporal dynamics in the rate estimates aligned with biological expectations: short-term rate estimates, based on recent diversification events, were higher than long-term estimates, reflecting the systematic decline in apparent rates as one looks further back in time. The power-law decay pattern that had been proposed to explain this discrepancy emerged naturally from the spline analysis, without being explicitly imposed.

SARS-CoV-2: Spatial Diffusion Across Europe

The second application turned to a much shorter timescale: the first ten months of the SARS-CoV-2 pandemic in Europe, from December 2019 to October 2020. Here, the question wasn't about the rate of nucleotide substitutions — the standard molecular clock application — but about the rate of spatial diffusion: how quickly the virus was spreading geographically across European countries.

The spline clock model found strong time-varying signal in the diffusion rate. The pattern wasn't monotonic; it reflected the shifting dynamics of the pandemic as it swept across Europe. Early diffusion rates, capturing the initial spread from early epicenters, were elevated. As travel restrictions, lockdowns, and behavioral changes took effect, the rate of cross-border transmission appears to have modulated. The smooth spline representation captured these dynamics without imposing artificial breaks.

The MCC tree from this analysis showed the phylogenetic relationships among European SARS-CoV-2 sequences, with nodes and branches colored by country. The spatial structure was evident: many clades clustered by geography, reflecting the fact that most transmission chains occurred within countries rather than across borders. But there were also clear examples of international spread, documented by the phylogenetic connections between lineages in different countries. The time-varying diffusion rate provided context for understanding when these cross-border movements were most common.

SARS-CoV-2 Spatial Diffusion Rate in Europe (2020)

Label	Value
Dec 2019	0.3 relative diffusion rate
Mar 2020	1.2 relative diffusion rate
Apr 2020	0.4 relative diffusion rate
Jun 2020	0.7 relative diffusion rate
Oct 2020	0.9 relative diffusion rate

Posterior median rate with 95% Bayesian credible interval. Strong time-varying signal detected throughout the European epidemic.

Why This Changes Things

Beyond the Metronome Assumption

The molecular clock hypothesis, proposed by Zuckerkandl and Pauling in 1965, has been one of the most fruitful simplifying assumptions in evolutionary biology. If evolution proceeds at a constant rate, then the genetic distance between two sequences is proportional to the time since their most recent common ancestor. This makes it possible to calibrate trees, to estimate divergence times, to reconstruct historical scenarios from contemporary sequence data. The assumption transformed phylogenetics from a qualitative discipline into a quantitative one.

But evolution doesn't always cooperate with our simplifying assumptions. The time-dependent rate phenomenon reveals that constant-rate molecular clocks are sometimes artifacts of analysis rather than reflections of biology. When rates genuinely vary over time, methods that assume constancy produce averaged estimates that may be misleading — not wrong, exactly, but capturing a mean that obscures all the interesting variation. Short-term and long-term estimates diverge because they're measuring different things.

The spline clock model represents a methodological response to this recognition. Rather than averaging over time-varying dynamics, it models them directly. The approach doesn't abandon the molecular clock framework; it enriches it. The rate is still a function of time, but it's a flexible, data-driven function rather than a single constant. This matters for downstream inference: if you want to know when a particular lineage emerged, or how quickly a pandemic accelerated, the answer depends on the rate, and the answer is more accurate when the rate model is more accurate.

Implications for Infectious Disease Research

The applications to foamy virus and SARS-CoV-2 illustrate the broad relevance of the approach. For long-term viral evolution — studying the deep history of retroviruses, endogenous viral elements, and other ancient pathogen-host relationships — the spline clock provides a way to accommodate time-dependent rates without imposing arbitrary epoch boundaries or accepting the averaged estimates produced by constant-rate models. The reconciliation of short-term and long-term rate estimates, documented for foamy virus, suggests that time-varying models may resolve long-standing discrepancies in viral evolution studies.

For contemporary pandemic dynamics — the rapid, acute spread of pathogens like SARS-CoV-2 — the method offers a way to trace how transmission dynamics change over the course of an epidemic. The analysis of SARS-CoV-2 spatial diffusion in Europe didn't just estimate a single average diffusion rate; it traced how that rate varied over time, reflecting the impact of travel restrictions, lockdowns, and behavioral changes. This is information that standard phylogenetic approaches, which assume constant rates, simply can't provide.

The ability to estimate time-varying rates also has implications for real-time outbreak analysis. As an epidemic unfolds, researchers want to know: is the transmission rate increasing or decreasing? Are new lineages emerging faster or slower? Is the virus adapting to its hosts? Standard clock models, calibrated on accumulated sequences, can track how average rates change over months or years, but they struggle with the rapid dynamics of acute outbreaks. A smooth time-varying model can potentially detect shifts in evolutionary dynamics more quickly, as the signal accumulates over weeks rather than requiring months of data.

Connections to Broader Statistical Methodology

The technical framework draws on a rich tradition of semiparametric and nonparametric function estimation. B-splines have long been used in statistics for their ability to approximate complex functional forms while maintaining smoothness constraints. The Bayesian P-splines approach — imposing random walk priors on spline coefficients to enforce smoothness — synthesizes the flexibility of spline-based function estimation with the principled uncertainty quantification of Bayesian inference. Gaussian Markov random fields provide the mathematical structure for these priors, connecting to a broader literature on spatial and spatiotemporal statistics.

The computational approach — using Gauss-Legendre quadrature to approximate branch integrals under a smooth spline parameterization — is also noteworthy. Numerical integration in phylogenetic models is challenging because the integrals appear in the likelihood function, which is evaluated millions of times during MCMC sampling. The deterministic, fixed-cost nature of the quadrature scheme ensures predictable runtime and reliable convergence, avoiding the adaptive integration problems that can arise with more flexible but less predictable numerical methods.

What's Next

Caveats and Limitations

No model is perfect, and understanding its limitations is essential for appropriate application. The spline clock model requires choosing the number and placement of knots, which governs the flexibility of the estimated rate trajectory. The researchers used a modest number of equally spaced knots combined with the GMRF prior to control smoothness, but the choice is still a modeling decision that could influence results in borderline cases. More extensive simulation studies, varying the number and placement of knots, would help characterize how sensitive inferences are to these choices.

The approach also requires the rate function to be positive and integrable, a technical assumption that ensures the mathematical machinery of continuous-time Markov chains applies. The exponential link function satisfies this assumption, but it means the spline coefficients are coefficients of the log-rate rather than the rate itself. This is interpretable — the coefficients represent log-rate changes — but it requires some care in interpretation. The squared link, which would model the rate directly, caused identifiability and mixing problems that the researchers chose to avoid.

The GMRF prior with a first-order random walk structure enforces smoothness but also imposes a particular structure on the rate trajectory. It assumes the rate changes gradually rather than jumping dramatically. For many biological applications, this is biologically reasonable — evolutionary rates don't typically change discontinuously. But there may be cases where sudden shifts in rate — reflecting dramatic changes in selection pressure or population dynamics — are real and worth detecting. The spline model might smooth over such shifts, though the tightness of the credible intervals in the simulation study suggests the model can still capture rapid changes when they're supported by the data.

Finally, the Bayesian framework provides principled uncertainty quantification but requires specifying prior distributions for all parameters. The priors used in the paper — normal priors on the intercept, GMRF priors on the spline coefficients, gamma priors on the precision — were chosen based on standard practices and showed good performance in simulation. But prior sensitivity analysis, exploring how results depend on prior choices, would strengthen confidence in the robustness of the approach.

Open Questions and Future Directions

The applications in the paper raise several directions for future research. The foamy virus analysis documented strong time-varying signal over 100 million years, but the biological interpretation of the rate trajectory deserves deeper investigation. What selective pressures or population dynamics drove the apparent rate decline over deep time? How does the spline-based rate estimate compare to paleontological estimates of divergence times? Connecting the rate dynamics to underlying biology, rather than just documenting their existence, is the next step.

The SARS-CoV-2 spatial diffusion analysis applied the spline clock to a phylogeographic model — tracking not just who infected whom but where transmissions occurred. The framework can accommodate any discrete trait evolving along the tree, so the approach could be extended to study the temporal dynamics of host jumps, the geographic spread of antimicrobial resistance, or the diffusion of viral subtypes across tissue types within a single host. The spline clock provides a general framework for any context where trait evolution rates might vary over time.

The method could also be extended to accommodate uncertainty in the substitution model — the researchers used standard HKY models, but more complex substitution processes with across-site rate variation are common in practice. The integration scheme should extend naturally, but numerical performance in more complex settings remains to be characterized. Similarly, extending the framework to multivariate rate functions — where different aspects of evolution (substitution, recombination, host switching) have potentially different time-varying dynamics — would broaden the scope of applications.

The Broader Significance

The spline clock model is part of a broader movement in computational biology and phylogenetics toward more flexible, realistic models that accommodate the full complexity of evolutionary processes. The assumption of constant rates, like all simplifying assumptions in science, was productive — it enabled tractable inference and generated insights — but eventually the field reaches the limits of what constant-rate models can tell us. Time-varying rate models represent the next step, and the spline framework provides a principled, computationally tractable way to take that step.

The implications extend beyond virology and phylogenetics. Any domain where molecular sequences are used to infer historical dynamics — microbial ecology, ancient DNA analysis, cancer genomics, metagenomics — depends on rate calibration. When rates are time-varying, the accuracy of historical inference depends on modeling that variation. The spline clock offers a general framework for doing so, implemented in widely used software, and could become a standard tool in the evolutionary biologist's toolkit.

At the same time, the work highlights a general principle: our assumptions shape what we can discover. The constant-rate molecular clock was an assumption so central that it became nearly invisible — a background condition rather than a hypothesis being tested. Time-varying rate models make the assumption explicit, ask whether it's justified, and provide a framework for answering that question with data. The discovery that rates genuinely vary over time isn't a failure of earlier methods; it's an advance enabled by more sophisticated modeling. As methods continue to improve, we should expect to find more such "artifacts of analysis" — assumptions embedded so deeply that we forgot they were assumptions — and to revise our understanding accordingly.