🧬 biology

Simulation and inference methods for non-Markovian stochastic biochemical reaction networks

本文通过将下一反应算法和 $\tau$ -leaping 算法推广以处理任意事件间隔分布，并引入一种能够通过多层蒙特卡洛和多保真度方法实现显著计算增益的耦合方案，为非马尔可夫生化反应网络开发了高效的随机模拟与推断方法。

原作者： Thomas P. Steele, David J. Warne

发布于 2026-08-20

📖 1 分钟阅读☕ 轻松阅读

原作者： Thomas P. Steele, David J. Warne

原始论文采用 CC BY 4.0 许可（http://creativecommons.org/licenses/by/4.0/）。 ⚕️ 这是一篇未经同行评审的预印本的AI生成解释。这不是医疗建议。请勿根据此内容做出健康决定。阅读完整免责声明

细胞内部是一个混乱、嘈杂的地方。分子在不断的随机风暴中碰撞、反应并发生转化。几十年来，科学家们一直使用数学模型来理解这种混沌，将这些分子相互作用视为一系列瞬时的跳跃。在这些传统的模型中，一旦反应开始，它就会立即完成，系统仅根据当前所处的位置向前推进，而忘记了之前发生的一切。这种被称为马尔可夫过程（Markovian process）的方法在理解许多生物系统方面非常有用。然而，它无法捕捉像基因转录这样复杂过程的现实，在转录过程中，从 DNA 指令到生成最终蛋白质的过程需要一段可衡量的时长。在这些情况下，反应的历史至关重要；系统会“记住”它在某项任务上已经投入了多久，而这种记忆会改变任务完成的可能性。

当科学家试图模拟这些“非马尔可夫”系统（即时间和历史起关键作用的系统）时，数学变得异常困难。用于模拟这些网络的标准工具运行速度过慢，难以投入实际应用，而且从实验数据中推断出系统的隐藏规则也几乎变得不可能。这为试图理解从细胞分化到疾病传播等现象的研究人员制造了瓶颈。昆士兰科技大学的一个研究小组现在开发了一套新工具来打破这一障碍。他们创建了更快速的方法来模拟这些复杂的、依赖记忆的系统，并开发了一种巧妙的方法，能够比以前更高效地推测这些模型的隐藏参数。

问题的核心在于生物反应是如何计时的。在一个简单的模型中，反应发生在条件成熟的那一刻。但在现实世界中，像基因转录成信使 RNA 这样的反应涉及一系列需要时间的步骤。研究人员意识到，为了准确建模，他们必须将反应视为一个具有起点和终点的过程，其中两者之间的时间间隔不是固定的，而是取决于细胞的当前状态以及该过程已经运行了多久。为了处理这个问题，他们改进了一种经典的模拟技术，称为“下一反应法”（next reaction method）。这种方法通常追踪系统中下一个化学事件何时发生。研究人员对其进行了修改，使其能够追踪每一个正在进行的反应，记录下每个反应的开始时间和已运行的时长。这使得计算机能够精确计算出延迟反应何时完成，即使完成过程的规则会随着细胞的变化而改变。

虽然这种精确方法很准确，但由于它必须单独计算每一个反应的每一个步骤，因此仍然非常缓慢。为了提高速度，该团队开发了第二种近似方法。这种方法是“ $\tau$ -跳跃”（tau-leaping）技术的延伸，它通过观察系统在微小时间段内的表现，而不是逐个瞬间进行观察。它不再追踪每一个分子的碰撞，而是估算在那个短时间窗口内会发生多少次反应。这种近似法速度极快，但也引入了少量的误差。突破点在于研究人员找到了将这两种方法联系起来的方法。他们创建了一个“耦合方案”（coupling scheme），一种让慢速、精确的模拟与快速、近似的模拟并排运行，并使用相同的随机事件的方法。因为它们运行在相同的随机路径上，所以快速模拟中的误差与精确模拟高度相关。这种相关性是实现高效的关键。

通过利用这种联系，研究人员可以应用一种强大的统计技巧，称为“多层级蒙特卡洛法”（multilevel Monte Carlo）。他们不再需要运行数百万次昂贵的精确模拟来获得精确答案，而是运行少量的精确模拟和大量的廉价近似模拟。快速模拟提供了大部分数据，而慢速模拟则提供了一个消除偏差的修正项。这使得他们能以极小的计算成本获得慢速方法的精度。他们在两个截然不同的生物场景中测试了这种方法。第一个是基因调节网络，其中一种蛋白质会抑制自身的产生，这一过程自然地产生了振荡或节律。第二个是疾病传播模型，其中感染者康复所需的时间不是固定的，而是遵循特定的分布。在这两种情况下，他们的新方法产生的结果都与精确的慢速模拟相匹配，但速度要快得多。

这项工作的真正威力体现在尝试从真实数据中推断系统的隐藏规则时。在生物学中，科学家通常拥有在特定时间点存在的分子数量的测量值，但他们并不了解反应发生的精确速率。为了找到这些速率，他们通常必须运行成千上万次模拟，以观察哪组规则产生的实验数据看起来最接近真实世界。这在计算上是非常耗时的。研究人员将他们这种耦合模拟方法应用于这一问题，使用了名为“多保真度近似贝叶斯计算”（multifidelity approximate Bayesian computation）的技术。他们将自己的方法与标准的、缓慢的方法进行了对比。结果令人震惊。他们的新方法在估计隐藏参数方面达到了与标准方法相同的准确度，但其计算时间减少了约五倍。在某些对比中，提速效果更为显著，达到了一个数量级的提升。

研究人员证明，他们的方法不仅适用于简单的模型，也适用于那些反应时间取决于系统状态和流逝时间的复杂系统。他们表明，通过仔细管理快速模拟与慢速模拟之间的关系，可以大幅缩减理解这些系统所需的时间。这为将这些详细、真实的模型应用于更大、更复杂的问题打开了大门，例如模拟整个细胞，或者在具有更现实康复时间的情况下模拟流行病的传播。这项工作并不声称解决了生物建模中的所有问题，但它提供了一个稳健且高效的基础。它证明了通过以聪明的方式结合精确方法与近似方法，科学家现在可以模拟并理解作为生命基础的、依赖历史的非马尔可夫过程，而不受计算能力的限制。

技术摘要：非马尔可夫随机反应网络的模拟与推断方法

问题陈述
随机反应网络模型对于捕捉复杂生物系统（如基因转录和翻译）中的内在噪声至关重要。虽然标准的表述依赖于马尔可夫过程（即未来动态仅取决于当前状态），但许多生物过程表现出具有历史依赖性的动态特性，这需要使用非马尔可夫建模。这些过程通常涉及随机延迟，其中反应启动到完成之间的持续时间取决于系统状态和时间。

主要的挑战在于非马尔可夫系统在计算上的不可行性，包括模拟和参数推断。

模拟： 非马尔可夫系统的精确模拟方法（例如非马尔可夫吉普斯算法）计算成本高昂，特别是当延迟依赖于系统状态时。此外，现有文献缺乏高效的近似模拟方案（类似于 $\tau$ -leaping）来与精确方法相结合。
推断： 对于这些系统，贝叶斯推断受到化学主方程（CME）导致的似然函数不可解性的阻碍。因此，研究人员依赖于无似然的、基于模拟的推断方法（例如近似贝叶斯计算，ABC）。然而，这些方法计算密集，需要大量的模拟。
多保真度差距： 加速策略，如多层蒙特卡洛（MLMC）和多保真度方案（结合了廉价的低保真度近似与昂贵的高保真度模拟），已成功应用于马尔可夫系统。然而，由于缺乏准确的近似模拟器，并且最关键的是，缺乏用于在精确模型和近似模型之间生成正相关样本路径的精确耦合方案，这些方法在非马尔可夫系统中的应用受到了限制。

方法论
作者开发了一个综合框架，用于模拟和推断具有状态和时间依赖性延迟的良构（well-mixed）非马尔可夫随机反应网络。

数学公式化：
作者扩展了库尔茨（Kurtur）随机时间变化表示法，以适应延迟反应。他们定义了一个系统状态向量，该向量不仅包括化学物种计数，还包括正在进行的延迟反应计数（ $X^*_j(t)$ ）。这使得系统能够追踪每个启动延迟的内部时间，从而能够对非指数等待时间分布进行建模。
精确随机模拟 (nM-NRM)：
作者将下一反应法（NRM）推广到非马尔可夫设定（算法 3）。

该方法不再采样全局等待时间，而是为每个反应通道（包括瞬时启动和延迟启动）以及每个活跃的延迟完成过程维护独立的单位速率泊松过程。
它利用依赖于当前状态和进行中延迟内部时间的累积风险函数。
算法在每次反应事件发生后，都会动态更新所有活跃通道的“剩余间隔”，从而实现对依赖于状态和时间的任意事件间时间分布的精确模拟。

近似随机模拟 (nM-TLM)：
为了实现加速，作者推导了一种非马尔可夫 $\tau$ -leaping 方法（算法 4）。

该方法假设系统状态在一个时间步长 $\tau$ 内保持不变。
它通过在固定间隔 $\tau$ 内对强度函数（propensities）进行积分来近似累积风险函数。
至关重要的是，它在时间步长末端以“组”（批量到达）的形式处理延迟启动，将 $\tau$ 内内的多次启动视为具有相同的启动时间。这种离散化允许生成的近似样本路径比精确方法计算成本更低。

耦合方案：
作者引入了一种新型耦合机制（详见附录 B），用于生成精确（nM-NRM）和近似（nM-TLM）样本路径的正相关对。

两种模拟都由相同的底层不均匀泊松过程（共同随机数）驱动。
这种耦合确保了精确路径与近似路径之间的差异被最小化，这是实现 MLMC 和多保真度推断所需的方差缩减的前提条件。

推断框架：
作者将这些模拟工具应用于多保真度近似贝叶斯计算 (MF-ABC)。

他们使用 nM-TLM 作为低保真度模拟器，使用 nM-NRM 作为高保真度模拟器。
利用耦合方案，他们实现了一种权重策略（公式 11），该策略可以纠正由低保真度近似引入的偏差，从而产生无偏的后验分布估计量。

核心贡献

广义精确模拟： 一种非马尔可夫扩展的下一反应法（nM-NRM），能够处理依赖于时间及系统状态的延迟，推广了以往仅限于固定或仅依赖于时间的延迟的研究。
非马尔可夫 $\tau$ -Leaping： 首个支持任意事件间时间分布的非马尔可夫网络近似随机模拟算法（nM-TLM），提供了一个计算高效的低保真度替代方案。
精确耦合方案： 一种用于生成非马尔可夫系统下正相关精确与近似样本路径的新型方法。这填补了一个关键空白，使得首次将 MLMC 和多保真度方法应用于非马尔可夫设定成为可能。
多保真度推断： 证明了这些模拟工具可以集成到 MF-ABC 中，以显著加速非马尔可夫模型的参数推断。

结果
作者通过两个案例研究验证了其方法：

延迟自抑制基因调节： 一个具有状态依赖性延迟的 mRNA 转录和蛋白质翻译模型。
- 收敛性： 当延迟分布变为指数分布（ $\alpha \to 1$ ）时，nM-NRM 正确地退化为标准的马尔可夫 NRM。
- 近似误差： nM-TLM 表现出预期的弱收敛（ $O(\tau)$ ）和强收敛（ $O(\tau^{1/2})$ ）速率。虽然近似会引入路径误差（例如，较大的 $\tau$ 会高估蛋白质峰值水平），但耦合方案在比较精确路径与近似路径时能有效降低方差。
非马尔可夫 SIS 流行病模型： 一个感染和恢复过程遵循 Weibull 等待时间分布的模型。
- MLMC 效率： 使用无偏 MLMC 估计量，作者实现了与直接蒙特卡洛采样相比的数量级的计算效率提升（MSE vs. 代价）。研究发现，由于更强的耦合效应，较小的步长 $\tau$ 在此设定下更为高效。
推断性能：
- 在基因调节模型中，MF-ABC 方案（使用 nM-TLM 作为低保真度并使用 nM-NRM 进行修正）在达到等效后验均值精度的情况下，比标准 ABC 拒绝采样实现了 5 倍的计算加速。
- MF-ABC 方法成功恢复了未知参数（ $\beta_{m*}, \beta_p, \alpha$ ）的后验边缘概率密度，与金标准 ABC 估计相比，实现了 10 倍的计算加速。

意义与主张
本文声称提供了一个基础框架，消除了在实践中实现非马尔可夫模型（在系统生物学及其他领域）的关键障碍。通过将下一反应法和 $\tau$ -leaping 扩展到非马尔可夫设定，并且最重要的是，建立了一个有效的耦合方案，作者使得在这些复杂系统中应用强大的方差缩减技术（MLMC 和多保真度方法）成为可能。

作者强调，他们的贡献并不局限于所展示的特定 ABC 实现。这些模拟算法和耦合方案是通用的，适用于任何依赖于精确与近似模拟耦合的基于模拟的推断方法（包括顺序蒙特卡洛、神经后验估计或其他 MLMC 变体）。这项工作表明，这些方法可以使全细胞模型以及其他因非马尔可夫随机模拟成本过高而难以处理的计算密集型定量生物学问题变得可行。

技术摘要：非马尔可夫随机反应网络的模拟与推断方法

类似论文