2025-11-19T01:43:13.983106

Inhomogeneous continuous-time Markov chains to infer flexible time-varying evolutionary rates

Datta, Lemey, Suchard

Reconstructing evolutionary histories and estimating the rate of evolution from molecular sequence data is of central importance in evolutionary biology and infectious disease research. We introduce a flexible Bayesian phylogenetic inference framework that accommodates changing evolutionary rates over time by modeling sequence character substitution processes as inhomogeneous continuous-time Markov chains (ICTMCs) acting along the unknown phylogeny, where the rate remains as an unknown, positive and integrable function of time. The integral of the rate function appears in the finite-time transition probabilities of the ICTMCs that must be efficiently computed for all branches of the phylogeny to evaluate the observed data likelihood. Circumventing computational challenges that arise from a fully nonparametric function, we successfully parameterize the rate function as piecewise constant with a large number of epochs that we call the polyepoch clock model. This makes the transition probability computation relatively inexpensive and continues to flexibly capture rate change over time. We employ a Gaussian Markov random field prior to achieve temporal smoothing of the estimated rate function. Hamiltonian Monte Carlo sampling enabled by scalable gradient evaluation under this model makes our framework computationally efficient. We assess the performance of the polyepoch clock model in recovering the true timescales and rates through simulations under two different evolutionary scenarios. We then apply the polyepoch clock model to examine the rates of West Nile virus, Dengue virus and influenza A/H3N2 evolution, and estimate the time-varying rate of SARS-CoV-2 spread in Europe in 2020.

academic

Inhomogeneous continuous-time Markov chains to infer flexible time-varying evolutionary rates

基本信息

论文ID: 2510.11982
标题: Inhomogeneous continuous-time Markov chains to infer flexible time-varying evolutionary rates
作者: Pratyusa Datta (UCLA), Philippe Lemey (KU Leuven), Marc A. Suchard (UCLA)
分类: stat.ME (Statistics - Methodology), q-bio.PE (Quantitative Biology - Populations and Evolution)
发表时间: 2025年10月13日 (arXiv预印本)
论文链接: https://arxiv.org/abs/2510.11982

摘要

本文提出了一个灵活的贝叶斯系统发育推断框架，通过将序列字符替换过程建模为非齐次连续时间马尔可夫链(ICTMCs)来适应随时间变化的进化速率。该方法将进化速率参数化为具有大量时期的分段常数函数(多时期时钟模型)，使转移概率计算相对廉价且能灵活捕捉速率变化。采用高斯马尔可夫随机场先验实现估计速率函数的时间平滑，并通过可扩展梯度评估的哈密顿蒙特卡洛采样提高计算效率。

研究背景与动机

问题定义

系统发育学中的核心问题是从分子序列数据重建进化历史并估计进化速率。传统方法假设进化速率在时间上保持恒定，但这一假设对于快速进化的病毒等生物体并不成立。

重要性

进化生物学意义：准确估计时变进化速率对理解生物多样化机制至关重要
传染病研究价值：病毒基因组序列在短时间尺度上积累显著遗传变化，需要实时分析能力
时间尺度依赖性：研究表明病毒的进化速率估计严重依赖于采样时间框架

现有方法局限性

齐次CTMC假设：传统方法假设分支上的替换过程为齐次连续时间马尔可夫链
速率变异模式固化：现有松弛时钟模型对速率变异模式做出固定假设
计算复杂性：完全非参数函数方法面临计算挑战

研究动机

开发能够直接将进化速率建模为时间函数的灵活框架，克服齐次CTMC假设的限制，为快速进化病毒等提供更准确的进化速率估计。

核心贡献

理论创新：首次将非齐次连续时间马尔可夫链(ICTMCs)系统性引入系统发育推断
方法突破：提出多时期时钟模型，将速率函数参数化为大量时期的分段常数函数
计算优化：开发线性时间复杂度的梯度评估算法，结合HMC实现高效采样
先验设计：采用适当的高斯马尔可夫随机场先验确保后验分布的适当性
实证验证：在多个病毒数据集上验证方法有效性，包括SARS-CoV-2传播分析

方法详解

任务定义

输入：N个排列的分子序列，采样时间信息输出：系统发育树、时变进化速率轨迹、分歧时间估计约束：速率函数必须为正值且可积

模型架构

1. ICTMC基础框架

对于非齐次CTMC，无穷小生成矩阵为时间函数： $Q(t) = f(t)Q$ ，其中：

$Q$ ：与时间无关的基础无穷小生成矩阵
$f(t)$ ：未知的正值可积速率函数

有限时间转移概率矩阵： $P(t_0, t) = \exp\left[\int_{t_0}^t f(\tau)d\tau \cdot Q\right]$

2. 多时期时钟模型

将速率函数参数化为分段常数： $f(t) = \theta_m, \quad w_m \leq t < w_{m-1}, \quad m = 1,\ldots,M$

其中 $w_M < \cdots < w_1$ 为时间网格点， $\theta = (\theta_1,\ldots,\theta_{M+1})$ 为速率参数向量。

3. 分支长度计算

对于连接节点 $i$ 到 $pa(i)$ 的分支，期望替换数为： $b_i = \theta_{q+1}(w_q - t_{pa(i)}) + \sum_{m=p}^{q-1}\theta_{m+1}(w_m - w_{m+1}) + \theta_p(t_i - w_p)$

4. 贝叶斯推断框架

先验设计：

对 $\zeta_m = \log\theta_m$ 使用高斯马尔可夫随机场先验
一阶差分： $\zeta_{m+1} - \zeta_m | \tau \sim N(0, d_m/\tau)$
适当先验： $P(\zeta|\tau) \propto \tau^{M/2}\exp[-\frac{\tau}{2}\zeta'(D_w - \rho W)\zeta]$

后验采样：使用哈密顿蒙特卡洛方法，利用链式法则计算梯度： $\frac{\partial}{\partial\theta_m}\log P(\theta,\tau,\rho,Q,\alpha,F|Y) = \sum_{i=1}^{2N-2}\frac{\partial\log P}{\partial b_i}\frac{\partial b_i}{\partial\theta_m}$

技术创新点

适当性保证：通过引入参数 $\rho < 1$ 确保GMRF先验的适当性
梯度优化：开发 $O(NCS^2 + NM)$ 复杂度的梯度计算，显著优于传统 $O(N^2CS^2)$ 方法
灵活网格设计：支持等间距或自适应网格点设置
多尺度建模：可处理从周到世纪的不同时间尺度

实验设置

数据集

仿真数据：
- 严格时钟模型仿真
- 对数线性时钟模型仿真( $f(t) = e^{-4.5-0.05t}$ )
真实病毒数据集：
- 西尼罗病毒：104个全基因组(1999-2007)
- 登革病毒3型：352个序列(1972-2010)
- 季节性流感A/H3N2：402个序列(1968-2010)
- SARS-CoV-2：3959个基因组(2020年欧洲)