系统生物学中动态模型的贝叶斯参数估计

系统生物学中动态模型的贝叶斯参数估计

Bayesian Parameter Estimation for Dynamical Models in Systems Biology

摘要
动力系统建模，特别是通过常微分方程组，已被有效用于捕捉信号转导网络中不同生化组分的时间动态行为。尽管近年来实验测量技术取得了显著进展，包括传感器开发和各类“组学”（-omics）研究，这些进展极大地丰富了蛋白质-蛋白质相互作用网络的细节，但系统生物学建模仍缺乏系统性的方法来估计动力学参数并量化相关不确定性。这源于多方面原因，包括实验测量数据稀疏且含噪声、反应背后缺乏详细的分子机制信息，以及缺失的生化相互作用。此外，微分方程系统在状态变量和参数方面固有的非线性进一步加剧了参数估计的挑战。在本研究中，我们提出一个全面的贝叶斯参数估计框架，以完整量化数据和模型中不确定性的效应。我们将该方法应用于一系列数学复杂度逐步增加的信号转导模型。对这些动力系统的系统性分析表明，参数估计依赖于数据稀疏性、噪声水平以及模型结构（包括多重稳态的存在）。这些结果强调，有针对性的不确定性量化能够丰富系统生物学建模，并为参数估计提供额外的定量分析能力。

关键词：参数估计，系统生物学，不确定性量化，贝叶斯方法

作者摘要
生物信号转导网络的数学模型已被广泛用于捕捉此类系统的时间动态行为。将这些模型与日益丰富的实验数据进行校准，对于确保模型准确反映生物现象至关重要。然而，测量噪声、无法测量系统中所有生化物种，以及对所有反应缺乏详尽了解，使得模型校准十分困难，并可能引入误差。在本研究中，我们提出一个原则性强且完整的计算框架，以应对这些挑战下的模型校准问题，并量化校准后模型中的各类不确定性（潜在误差）。我们将该框架应用于一系列示例模型，这些模型展示了生物学中常见的多种动态行为（如极限环、稳态），以说明我们的方法如何为基于模型的研究提供额外的背景信息和深入见解。

1. 引言

   数学建模是系统生物学不可或缺的组成部分；事实上，动力系统分析方法的应用促使我们对生化信号转导的理解发生了范式转变，并使得人们能够识别信号网络的涌现特性 [1, 2, 3, 4]。此外，数学模型使我们能够研究生物系统的动态行为，其范围超出了实验所能达到的限度 [5, 6, 7, 8]。对信号转导动态行为进行建模的一种经典方法是使用常微分方程组（ODEs）[9, 10]。这些方程通常包含非线性函数，以利用米氏动力学（Michaelis-Menten kinetics）和希尔函数（Hill functions）描述协同结合等复杂的生化相互作用 [11]。然而，构建和约束具有预测能力的信号转导模型所面临的一个持续挑战，是如何估计和识别与这些反应相关的动力学参数，并量化与之相关的不确定性 [12, 13, 14]。尽管在更广泛的计算科学领域中，不确定性量化（Uncertainty Quantification, UQ）领域早已广泛采用严格的定量方法来估计动力学参数及其不确定性 [19, 20, 21, 22]，但在系统生物学中，这类方法的应用仍处于早期阶段 [15, 12, 16, 17, 18]。

在信号转导的动力系统建模中，存在多种不确定性来源，包括模型结构本身、模型参数的取值，以及用于模型校准的数据质量。模型方程中的不确定性，即所谓的模型形式或拓扑不确定性（model form or topological uncertainty）[23, 24, 25, 26]，通常在模型开发过程中产生。然而，许多生化反应的反应通量（即ODE的表达形式）可以根据经典速率方程建立 [9, 10, 11]。更具挑战性的问题是为各种通量项确定合适的模型参数，包括动力学速率常数 [16]。这些参数的直接测量通常在孤立的反应体系中进行，无法反映模型所代表的整个反应网络的复杂性。要估计这些生物参数并识别剩余的不确定性，需要选择一个统计模型，然后从可用数据中学习这些参数的分布 [19, 27]。从这一视角来看，生物参数是随机变量，其分布可以是参数化的，也可以是非参数化的。然而，参数估计因系统生物学中数据的噪声性、稀疏性（时间点较少）和不完整性（由于实验限制，只能观测到部分或特定输出）而变得复杂，这些因素都会在生物参数中引入不确定性 [28, 29, 14, 30, 31]。面对这些复杂因素，亟需对参数进行统计建模，以实现不确定性量化。

一个全面的参数估计与不确定性量化（UQ）框架应考虑结构可识别性（structural identifiability）和参数敏感性（parameter sensitivity）的影响 [16, 29, 32, 33, 34]。结构可识别性分析可揭示在给定特定动力系统模型和一组可观测输出的情况下，哪些参数可以被估计 [28, 29, 35, 36]。若某个参数的每个取值都对应唯一的一组模型输出，则该参数被称为全局结构可识别（globally structurally identifiable）[29]。不符合该标准的参数则被视为结构不可识别（structurally nonidentifiable），无法从指定的模型输出中成功估计。结构不可识别性可能源于复杂的非线性方程以及仅能测量系统部分状态的不完整实验数据。此外，参数敏感性分析 [19, 37, 38] 可量化模型输出对模型参数变化的敏感程度。Gutenkunst 等人 [28] 发现，大多数系统生物学模型包含具有广泛敏感度范围的参数，他们将这类参数称为“松散”（sloppy）参数。尽管存在这一挑战，敏感性分析仍有助于根据参数对特定模型输出的贡献对可识别参数进行排序 [39]，例如 Mortlock 等人 [40] 在催乳素介导的 JAK-STAT 信号通路研究中所做的工作。该分析使他们能够从总共 60 个参数中筛选出 33 个对模型输出变化有显著贡献的参数。

系统生物学模型中常用的参数估计方法包括频率学派方法（frequentist）[17] 和贝叶斯方法（Bayesian）[18]。在频率学派框架下，参数估计被表述为一个优化问题，其解是一组能够最好地复现实验数据的参数 [19]。此外，频率学派方法通过围绕最优参数估计置信区间来量化不确定性 [41, 42, 19]。相比之下，在贝叶斯视角下，参数被视为随机变量，其未知的概率分布（称为后验分布）用于量化参数空间中任一取值的概率 [19, 21, 27]。贝叶斯方法的优势在于能够刻画整个后验分布，并通过可信区间（credible intervals）量化参数估计中的不确定性 [27]。已有多种贝叶斯参数估计方法被提出 [18, 43, 44, 45, 46, 47, 48, 49]，它们均旨在利用贝叶斯定理 [19, 27] 来刻画后验分布。例如，Mortlock 等人 [40] 成功运用贝叶斯估计研究了模型预测中的不确定性，并评估了其建模结果的统计显著性。

尽管贝叶斯参数估计在系统生物学中已取得一定成功 [12, 15, 16, 18]，但若未能考虑模型中所有不确定性来源，将严重阻碍参数估计和不确定性量化 [16, 26, 28]。因此，系统生物学中一个全面的 UQ 框架应严格考虑模型结构、不可识别参数、混合参数敏感性，以及噪声、稀疏或不完整的实验数据所带来的不确定性。虽然可识别性和敏感性分析通常在参数估计之前进行 [16, 29]，但要处理模型形式不确定性（model form uncertainty），则需要采用随机模型而非确定性模型 [23, 26, 50, 30]。一种有前景的处理模型形式不确定性的方法是无迹卡尔曼滤波马尔可夫链蒙特卡洛（UKF-MCMC）方法 [26, 51]。该方法同时包含对噪声数据和模型形式不确定性的统计建模；然而，它尚未被适配以应对系统生物学中的独特挑战。类似地，文献 [30] 中提出的参数估计与模型选择方法通过扩展卡尔曼滤波考虑了模型形式不确定性，但由于采用频率学派方法进行参数估计，未能提供完整的不确定性估计。因此，亟需一个综合框架，将结构可识别性分析、全局敏感性分析、针对数据和模型形式不确定性的统计模型，以及贝叶斯参数估计相结合，以实现系统生物学中动力学模型的全面不确定性量化。

本工作的创新之处在于我们开发了一套全面的工作流程，推动了系统生物学中不确定性量化领域的当前技术水平。此外，所提出的框架通过扩展贝叶斯框架和 UKF-MCMC 方法 [26, 51]，同时考虑了数据和模型结构中的不确定性。我们提出的工作流程首先进行结构可识别性分析 [36, 52] 和全局敏感性分析 [37, 53]，以剔除那些无法从测量数据中学习到或对模型输出无影响的参数，然后通过引入约束区间无迹卡尔曼滤波（Constrained Interval Unscented Kalman Filter, CIUKF）[54]，对 UKF-MCMC 方法进行系统生物学领域的扩展。总体而言，上述每一步都定量地应对了模型开发与校准过程中遇到的不确定性，从而提升系统生物学中的预测建模能力。

本文后续部分详细阐述了我们在系统生物学中用于不确定性量化的综合工作流程，并通过若干示例突出该分析方法。首先，在第 2 节中，我们介绍所提出的框架，并提供理解和应用该方法所需的数学细节。接着，在第 3 节中，我们将该框架应用于三个复杂度逐步递增的系统生物学模型，包括一个简单的双态模型 [55]、一个核心丝裂原活化蛋白激酶（MAPK）信号通路模型 [56]，以及一个用于描述长时程增强/抑制（long-term potentiation/depression）的突触可塑性现象学模型 [57]。我们发现，即使在简单模型中，参数估计也依赖于数据噪声水平和数据稀疏程度。最后，该框架能够对包含非线性和多稳态（multistability）的模型结构进行不确定性量化。在所有这些案例中，我们都利用可识别性和敏感性分析缩小待估计参数的子集，然后采用贝叶斯估计方法确定模型结构在参数估计中的作用。这些结果确立了一种以不确定性量化为核心的系统生物学方法，可实现严格的参数估计与分析。最后，在第 4 节中，我们结合上述三个示例讨论这些发现，探讨应用贝叶斯方法所面临的挑战，并为系统生物学中未来不确定性量化研究指明方向。

1. 方法

   本节描述了第 2.1 节中提出的不确定性量化综合框架的技术细节（见图 1）。随后，第 2.2 节介绍动力系统生物学模型及参数估计问题。第 2.3 节和第 2.4 节分别概述结构可识别性分析和全局敏感性分析，以降低参数空间的维度。接着，第 2.5 节介绍贝叶斯估计，第 2.6 节概述 CIUKF-MCMC 算法，第 2.7 节详细描述约束区间无迹卡尔曼滤波（Constrained Interval Unscented Kalman Filter, CIUKF）。此后，第 2.8 节讨论如何构建先验分布，第 2.9 节详述马尔可夫链蒙特卡洛（MCMC）采样。最后，第 2.10 节讨论通过集合建模（ensemble modeling）进行输出不确定性传播，第 2.11 节强调如何为参数选择点估计量，第 2.12 节说明合成数据生成方法，第 2.13 节概述极限环分析。

2.1 系统生物学中动力学模型的全面不确定性量化框架

本节预览了所提出的用于系统生物学中动力学模型参数估计与不确定性量化的综合框架。图 1 概述了该框架及其组成部分，后续各节将对这些部分进行更详细的描述。所提出的框架包含三个主要步骤。

首先，我们假设细胞内信号转导的动力学模型（图 1.A1）采用经典生化速率定律，例如质量作用动力学、米氏动力学和希尔函数 [9, 10, 11]（图 1.A2）。应用这些模型的关键挑战在于：如何根据可用的实验数据（图 1.A4）估计相关参数，例如图 1.A2 中的速率常数 k和 Vmax、平衡常数 Km和 KA，以及希尔系数 n。该综合框架采用贝叶斯推断，基于一组含噪声的测量数据和特定的模型形式，为模型参数估计一个统计模型（即概率分布；见图 1.B）。

我们认为，在参数估计之前（图 1.B1），必须进行可识别性分析和敏感性分析。为消除因不可识别参数带来的不确定性，我们使用结构可识别性分析器（Structural Identifiability Analyzer, SIAN）[36, 52] 进行全局结构可识别性分析（详见第 2.3 节）。不可识别的参数被固定为其文献中的标称值或基于其生理范围的合理值。接下来，进行基于方差的全局敏感性分析 [19, 37]，以根据各可识别参数对模型输出方差的贡献程度对其进行排序（详见第 2.4 节）。从中选取敏感性指数最大的可识别参数子集用于参数估计，其余模型参数则以与不可识别参数相同的方式固定为其标称值。

贝叶斯参数估计通过估计一个非参数统计模型，完整刻画模型参数中的不确定性。贝叶斯定理（见图 1.B2）从初始猜测（先验分布）出发，结合实验数据通过似然函数进行更新，从而得到对参数的最佳估计分布，即后验分布。似然函数衡量数据与模型预测之间的不匹配程度，并为那些产生更接近数据输出的参数集赋予更高的概率。我们采用 CIUKF-MCMC 算法 [26, 51] 来近似似然函数，并同时考虑模型形式、数据和参数中的不确定性。通过马尔可夫链蒙特卡洛（MCMC）采样——可选用延迟拒绝自适应 Metropolis 算法 [58] 或仿射不变集合采样器（affine invariant ensemble sampler）[59]——生成一组样本，用以表征后验分布（图 1.B3）。

最后，我们利用后验分布量化模型参数中的不确定性如何影响模型预测中的不确定性（图 1.C）。使用参数样本进行集合模拟（ensemble simulation），生成多组轨迹（见图 1.C1），以捕捉预测动态中的不确定性。计算如图 1.C2 所示的 95% 可信区间（credible intervals）等不确定性区间，可直观展示这种不确定性。值得注意的是，可信区间与置信区间不同：可信区间包含后验样本中指定百分比的取值，而置信区间本身是随机变量，表示在指定概率水平下估计量所处的区域 [27]（图 2 展示了一个可信区间的示例）。此外，对集合模拟结果进行统计分析，能够以类似于通过多次实验重复分析实验结果的方式，对计算建模结果进行定量分析。

在接下来的章节详细讨论实现全面不确定性量化的方法之前，下一节将正式介绍系统生物学中动力学模型的参数估计问题。

2.2. 以部分观测的常微分方程组形式表示的系统生物学模型的参数估计
我们考虑如下形式的非线性常微分方程模型：

我们现在可以定义所提出的参数估计框架的问题设置。这项工作假设模型形式已知，并寻求通过学习参数的概率分布来估计模型参数和相关不确定性。以下问题陈述形式化了参数估计问题。

下一节概述了在参数估计之前为降低参数向量 θ的维度而进行的结构可识别性分析和全局敏感性分析。

2.3. 使用结构可识别性分析器（SIAN）进行全局结构可识别性分析
结构可识别性分析用于判断参数是否能从可用的测量函数中被唯一估计出来 [29]。结构可识别性是模型本身的一种数学性质，不考虑可用实验数据的质量或数量。以下引自文献 [60] 的定义为全局结构可识别性提供了一个直观的判据，并可被证明等价于其他形式的定义 [33, 36]。

全局结构可识别性（Global structural identifiability），如定义1所述，意味着参数 θᵢ 可以从数据中被唯一确定。如果一个参数是全局结构可识别的，则存在唯一的参数取值，能在相同的初始条件下产生每一个观测到的输出值 [36]。反之，若定义1中的条件仅在参数空间 θᵢ 附近的局部邻域 (θᵢ) 内成立，则该参数可能是局部结构可识别的，例如当 θᵢ ∈ (θᵢ) 时 [29, 33, 36]。这一条件意味着，对于一个局部可识别的参数，有限个不同的参数取值可能产生相同的输出值 [36]。最后，如果一个参数是不可识别的，则无穷多个该参数的取值可以产生相同的模型输出。尽管包含局部结构可识别参数和不可识别参数的模型在整个参数空间内仍有效，但这类参数的存在会干扰成功的参数估计。目前已开发出多种计算方法用于评估基于常微分方程（ODE）模型的结构可识别性 [29, 33, 36, 61]。在本研究中，我们采用文献 [36, 52] 中提出的基于微分代数与幂级数的方法，因为该方法专门用于评估全局可识别性。SIAN（结构可识别性分析器，Structural Identifiability ANalyser）软件 [52] 提供了文献 [36] 所提出方法的数值实现。SIAN 用于评估全局结构可识别性的算法结合了符号计算与概率计算 [36, 52]。首先，SIAN 利用模型方程的泰勒展开式获得系统的多项式表示；其次，算法截断多项式系统，生成一个包含所有参数可识别性信息的最小系统；第三，SIAN 针对随机选取的单个参数集求解可识别性问题，确保结果在用户指定的概率水平 p 下正确（参见文献 [36] 中的定理5）；第四，算法利用第三步的结果将参数划分为全局可识别、局部可识别和不可识别三组。SIAN 已在 Maple（Maplesoft, Waterloo, ON）和 Julia（The Julia Project [62]）中实现。有关 SIAN 的更多数学细节，我们建议读者参考文献 [36]。

我们使用 SIAN 算法的 Julia 实现版本，默认正确性概率 p = 0.99（可在 https://github.com/alexeyovchinnikov/SIAN-Julia 获取）。此外，我们将附加参数 p_mod 设置为 2²⁹ − 3，以加快算法运行速度 [63]。任何经可识别性分析后未被判定为全局结构可识别的参数，均根据现有文献固定为其标称值。虽然这些参数可能传递有意义的生物学信息，但不可识别性意味着从现有数据中在数学上不可能识别它们。接下来，我们将使用全局敏感性分析进一步缩小可识别参数的集合。

2.4. 基于方差的全局敏感性分析
参数敏感性分析用于量化模型参数对模型输出变化的贡献 [19, 37]。具体而言，全局敏感性分析旨在量化模型参数在整个参数空间范围内对所关注输出量的影响 [19, 37]，特别适用于研究非线性模型中的参数 [37, 32]。本研究采用 Sobol 方法 [53] 进行基于方差的敏感性分析，因为它能提供一个定量的输出用于对参数进行排序，并且可利用为参数定义的先验分布。Sobol 敏感性分析根据单个参数的贡献以及参数之间的相互作用，对模型输出的方差进行分解 [37, 53]。输出量 f(θ)的总方差为

本研究使用 UQLab 工具箱 [64, 65] 在 MATLAB（MathWorks, Natick, MA）中执行 Sobol 敏感性分析，并使用 DifferentialEquations.jl 软件包 [66] 在 Julia 中进行分析。这两种软件均采用蒙特卡洛估计器，根据参数样本计算 Sobol 敏感性指数（详见 [19, 37] 及其中引用的文献）。在 MATLAB 中，除非另有说明，一阶敏感性指数和总阶敏感性指数均使用 Sobol 伪随机采样（例如设置 SOpts.Sobol.Sampling = 'sobol'）和默认估计器（SOpts.Sobol.Estimator = 'janon'）进行计算，具体细节参见 [65]。此外，采样数量（SOpts.Sobol.SampleSize）在第 3 节中针对每个具体问题分别设定。在 Julia 中，我们使用 QuasiMonteCarlo.jl（https://github.com/SciML/QuasiMonteCarlo.jl）进行类似的 Sobol 采样，采用 SobolSample() 采样器，并使用默认的敏感性指数估计器，例如将 Ei_estimator 设为 :jansen1999。

这种高斯分布可以在所有 x值处进行解析计算。因此，需要从后验分布中抽取参数样本，以刻画整个参数空间上的分布。马尔可夫链蒙特卡洛（Markov chain Monte Carlo, MCMC）算法能够从任意分布（例如后验分布）中进行采样（详见第 2.9 节）。在进行贝叶斯估计之前，下一节将介绍约束区间无迹卡尔曼滤波马尔可夫链蒙特卡洛（Constrained Interval Unscented Kalman Filter Markov Chain Monte Carlo, CIUKF-MCMC）算法，该算法能够同时考虑模型和数据中的不确定性。

2.6. 约束区间无迹卡尔曼滤波马尔可夫链蒙特卡洛（CIUKF-MCMC）
系统生物学中动力学模型的完整不确定性量化必须考虑模型形式、参数以及含噪声数据中的不确定性。Galioto 和 Gorodetsky 在文献 [26] 中建议，在方程 (1) 中引入一个过程噪声项，以表征系统中的模型形式不确定性。遵循这一建议，方程 (1–2) 中的模型被重新表述为一个离散时间随机过程：

定理 1 中定义的递推关系与贝叶斯滤波器 [67] 非常相似；因此，它可通过卡尔曼滤波算法进行评估 [26]。对于线性模型，可使用标准的线性高斯卡尔曼滤波器来计算该递推关系（参见文献 [26] 中的算法 2）。然而，若模型或测量过程为非线性，则无法获得该递推关系的精确解。因此，必须采用近似方法，例如扩展卡尔曼滤波器（EKF）、无迹卡尔曼滤波器（UKF）或集合卡尔曼滤波器（EnKF）[26, 67]。UKF-MCMC 算法的原始实现使用 UKF [68] 进行这种近似，因为 UKF 通常稳定且能处理非线性模型和测量过程 [26, 51]。然而，UKF 不适用于系统生物学模型，因为它忽略了状态变量上的约束（如化学浓度的非负性）。第 2.7 节将描述本工作中所采用的约束区间无迹卡尔曼滤波器（CIUKF）[54]，用于在滤波过程中强制执行状态约束。因此，我们称文献 [26] 中使用的、基于 CIUKF 的 UKF-MCMC 方法为 CIUKF-MCMC。

2.7. 约束区间无迹卡尔曼滤波器（CIUKF）

2.8 先验分布

先验分布 p(θ)编码了我们在收集数据和执行贝叶斯估计之前对参数的“信念状态”[19, 27]。先验分布的形式允许我们将不同程度的先验知识纳入模型中。如果某参数的取值已知，可使用信息性先验（informative priors）使先验分布向已知值集中 [27]；例如，可以使用对数正态先验分布，使先验围绕该参数的实验测量值进行中心化 [72]。然而，贝叶斯推断文献常提醒：信息性先验应仅与关于参数值的良好信息结合使用 [27]，因为克服一个“不良”先验可能需要大量数据。相比之下，无信息或弱信息先验反映的是对参数缺乏良好先验知识的情况。例如，若我们仅知道某参数的生理范围，则可构造一个均匀先验，表明该范围内任何参数值的概率相等。当参数的先验知识有限时，选择无信息或弱信息先验通常更为稳妥 [27]。

2.9 马尔可夫链蒙特卡洛采样

该接受概率保证了样本的平稳分布（即在样本数量趋于无穷时样本所收敛到的分布）等于目标分布 [27, 75]。上述“提议”和“接受-拒绝”步骤不断重复，直到样本集合的分布收敛至其平稳分布。尽管在样本数量趋于无穷时收敛性可得到理论保证 [27, 75]，但在实际应用中评估收敛性却并非易事。我们选择采用文献 [59] 中提出的一种方法，利用积分自相关时间（integrated autocorrelation time）来评估收敛性，具体如第 2.9.5 节所述（其他评估方法参见 [19, 27, 75]）。

2.9.3 仿射不变集合采样器（AIES）

尽管基于随机游走 Metropolis 的算法（如 DRAM）能够充分采样复杂的后验分布，但当目标分布高度各向异性时，这些方法会表现出极慢的收敛速度 [59]。在系统生物学中，使用 CIUKF-MCMC 得到的后验分布通常是各向异性的，因为模型参数的尺度可能相差几个数量级，且噪声协方差参数的尺度通常与模型参数不同。幸运的是，AIES 提供了一种能有效采样此类各向异性分布的算法 [59]。仿射不变性的动机在于：各向异性分布可通过仿射变换转化为各向同性分布。因此，一个对这类变换保持不变的算法，在采样各向异性分布时将能有效采样对应的各向同性分布 [59]。

2.9.4 马尔可夫链预热期（burn-in）

在 MCMC 中，马尔可夫链（或链的集合）在收敛到其平稳分布之前通常会经历一个初始暂态过程，称为“预热期”（burn-in）[19, 27, 76]。重要的是，这些预热期样本并不服从平稳分布，因此应从最终样本集中剔除。MCMC 文献中的常见做法是直接丢弃这些初始样本，以消除预热期的影响 [19, 27]。预热期长度的选择通常并非易事，最好通过分析马尔可夫链本身来确定 [76]。除非另有说明，本文中预热期样本数量由积分自相关时间（下文将描述）决定。具体而言，我们在收集大量样本后计算积分自相关时间，并将预热期长度设为该计算值的 5–10 倍。

2.9.5 利用积分自相关时间评估收敛性

马尔可夫链蒙特卡洛采样的一个关键挑战是确定应采集多少样本 N。在本研究中，我们使用积分自相关时间 [59, 76] 来判断马尔可夫链何时已近似收敛至其平稳分布。我们概述单条马尔可夫链积分自相关时间背后的理论与动机，并建议读者参阅文献 [59] 了解关于集合方法的讨论。使用积分自相关时间的动机源于 MCMC 采样在计算期望值时的典型应用

2.10 集合模拟与输出不确定性分析

人们通常希望在刻画完整后验分布（图 2 中的黑线）之外，再为每个参数计算一个点估计值。贝叶斯统计中常见的点估计选择包括后验分布的均值、中位数或众数 [27]。图 2 突出了均值和众数（记作 MAP），并显示了一个次要众数，该次要众数可能干扰对点估计量的选择。后验分布的众数是一个强有力的选择，因为它提供了最可能的参数集合，通常被称为最大后验概率（maximum a posteriori, MAP）点。然而，当后验分布为多峰分布（存在多个众数）时，均值或中位数可能提供更好的点估计。

这些点估计值是从通过 MCMC 采样获得的后验样本中计算得出的。虽然像均值和中位数这样的样本统计量可直接从样本中计算，但计算 MAP 点需要先估计概率密度函数的后验分布，并找到其最大值点。直接计算样本众数无法得到 MAP 点，因为参数是连续随机变量，因此预期不会出现两个完全相同的两参数样本。一种估算 MAP 点的方法是计算样本的直方图，并将具有最高关联概率的区间中心作为 MAP 值。然而，我们发现后验样本的直方图往往噪声较大，且对区间宽度敏感，导致 MAP 估计可能不准确。因此，我们选择采用另一种方法：使用核密度估计器 [84] 来近似后验分布，然后据此计算 MAP 点。核密度估计器提供了一种非参数化的后验分布近似，并直观地平滑了后验直方图。我们在 MATLAB 中使用 ksdensity() 函数进行核密度估计。“support”选项设为由先验分布所界定的区域，“BoundaryCorrection”选项设为使用“reflection”方法以处理边界效应。除非另有说明，其余选项均保留 MATLAB 文档中的默认设置。MAP 点即为 ksdensity() 输出中概率最大的点。下一节将从 CIUKF-MCMC 和贝叶斯估计的细节转向概述本文示例中合成数据的生成方法。

2.12 数值实验用合成数据生成

本工作中的合成数据旨在模拟生物实验中常见的含噪数据。我们通过从确定性模型模拟中抽取样本，并向每个样本独立添加同分布（iid）的扰动来生成含噪合成数据。首先，为数据生成选择一组名义上的生物模型参数和初始条件；这些值构成估计问题的“真实值”（ground truth），并在可能的情况下参考现有文献确定。接下来，使用 MATLAB 中的 ode15s() 积分器（默认容差）求解系统数值解，从而获得状态变量的真实轨迹；同时提供雅可比矩阵（指定 ‘Jacobian’ 选项，使用解析雅可比矩阵）以提高计算效率。随后，应用测量函数并对真实轨迹进行子采样，以模拟稀疏采样。最后，通过向每个样本添加一个零均值、正态分布的独立同分布噪声随机过程来污染数据。为满足 CIUKF-MCMC 算法的假设（见第 2.6 节），我们使用均值为零、正态分布的扰动，其协方差矩阵 Γ为对角阵。我们选择 Γ的元素与相应状态变量的方差成比例，例如：

其中 b是一个通常小于 1 的正常数，用于控制噪声水平。
最后一节讨论如何计算极限环振荡的振幅和周期。

2.13 极限环分析
极限环的振幅和周期被用于表征极限环振荡，以进行全局敏感性分析和输出不确定性分析。这些量在细胞内信号转导中具有重要意义，因为信号的强度（振幅）和时序（周期）被认为能够编码不同的输入信息 [85]。极限环振幅 用于量化振荡最大值与最小值之间的差异，其定义为：

3. 结果

我们将贝叶斯参数估计框架应用于三个不同的模型，这些模型代表了生物复杂性和数学复杂性逐步递增的信号转导级联过程。第 3.1 节使用第一个模型——一个简单的动力学方案——来描述一系列计算实验，以说明测量噪声和数据稀疏性对估计不确定性的影响（图 3）。接下来，我们在两个更能体现系统生物学模型中非线性和过参数化特征的模型上测试我们的框架。第 3.2 节分析了一个代表性的丝裂原活化蛋白激酶（MAPK）通路模型 [56]，该模型根据模型参数的选择表现出多稳态特性；附录 A.1 则使用该模型说明结构可识别性和全局敏感性分析的必要性。最后，第 3.3 节分析了一个简化的突触可塑性模型 [57]，该模型说明即使在现象学的生物学表示中，更高参数不确定性所带来的影响。

3.1 测量噪声和数据稀疏性增加简单信号转导模型中的估计不确定性

我们考虑的第一个模型是一个相对简单的两状态模型，源自文献 [55]，如图 3.A 所示。该模型的控制方程为：

备注 1. 数据点与 MCMC 样本不同。
本研究中的实验最多生成 40 个（模拟的）数据点，即用于参数估计的状态变量的含噪测量值。然而，MCMC 算法需要抽取数万至数百万组参数样本，以刻画后验分布，这要求反复评估似然函数，从而必须对模型进行多次模拟。

3.2 MAPK 级联模型的参数估计

进一步检查 30,000 条后验轨迹中的 50 条子集（图 6.C 中的绿色线条）发现，95% 可信区间内包含许多与真实轨迹（虚线黑线）相似的极限环，少量不匹配真实振荡的极限环，以及更少数量最终达到不动点的轨迹（例如图 6.C 中的点状蓝线）。

我们进一步利用后验样本量化预测极限环的变异性。首先，我们对整个集合中的所有轨迹（图 6.D）进行了表征，发现 90.6%（27,182 个样本）正确产生了极限环振荡，而 9.6%（2,818 个样本）收敛至不动点。对这 27,182 条振荡轨迹分别计算其极限环振幅和周期后发现，尽管模型参数的不确定性很小，预测极限环的特征仍存在明显变异性。我们发现，最可能的极限环振幅和周期（见图 6.E 中的直方图）接近真实值（振幅为 26.12，周期为 13.04；由垂直虚线黑线标示）。然而，即使参数估计具有高度确定性，这些量仍表现出一定变异性：具体而言，极限环振幅常偏离真实值最多达 30%（即 10 nM）；周期的偏离幅度较小，最多约 15%（即 2 分钟），但其分布偏向比真实周期更短的值。

基于以上发现，我们得出结论：模型参数的小幅不确定性会导致极限环振幅和周期出现变化，但不会显著影响我们预测振荡动力学的能力。

接下来，我们研究了用于估计的数据的性质是否会影响我们预测极限环振荡的能力。为了研究这一点，我们使用两个额外的含噪极限环数据集重复了参数估计（如图7.A所示，噪声水平设置为真实轨迹方差的1%）。我们将这些新数据称为等距采样数据和非等距采样数据。等距采样数据包括在0 < t ≤ 60分钟内每两分钟采集的20个样本，涵盖了初始瞬态过程和几个振荡周期（见图7.A）。非等距采样数据同样覆盖前60分钟，但在初始瞬态阶段（前30分钟）仅每五分钟采集5个样本，剩余30分钟内每两分钟采集15个样本（总计20个样本）。我们将此结果与仅包含振荡数据的结果（图6）进行比较，以研究来自初始瞬态阶段的样本如何影响估计。对于等距采样数据，我们使用150个walker进行了22,500步的MCMC（补充图C.22），而对于非等距采样数据，则使用30个AIES walker进行了36,500步（补充图C.23）。为考虑MCMC预热期的影响，我们丢弃了等距数据每个链的前5,669个样本以及非等距数据的前11,065个样本。此外，我们针对这两个新数据集，利用随机后验样本模拟了包含30,000条轨迹的集合；补充图C.15展示了由此产生的动力学后验分布。

对30,000条模拟轨迹中显示振荡的比例进行比较，清楚地表明输入数据影响我们预测极限环的能力（见图7.B）。作为基线，我们将这三个数据集与从先验分布抽取样本的模拟结果进行比较，我们认为这是未提供任何数据的情况。图7.B显示，所有有数据的情况预测出更高比例的极限环轨迹，而先验样本中仅有17.6%的模拟出现振荡。我们发现，对应于等距采样数据的模拟中，仅有45.9%（13,794个样本）显示出持续振荡。然而，随着来自初始瞬态阶段（0至30分钟的数据）的样本数量减少，振荡模拟的比例增加；对应于非等距采样和仅振荡数据的模拟中，分别有63.1%（18,921个样本）和90.6%（27,182个样本）显示出极限环振荡。

我们还发现，动力学性质的不确定性与模型参数的不确定性密切相关。图7.C-E分别显示了对应于等距采样、非等距采样和仅振荡数据的边际后验分布。具体而言，我们观察到，我们总是可以高置信度地预测k₂，但k₄、k₆和α的不确定性在三个数据集中有所不同。例如，对于等距采样数据，k₄的边际分布（图7.C）具有三个高概率的明显峰值（模态）：一个集中在非常小的值，另一个集中在标称值附近，还有一个更宽的模态位于更高的值。对于非等距采样数据，图7.D中的k₄边际分布将较低概率分配给远离标称值的两个模态，但靠近标称值的模态向更高的k₄值偏移。最后，对于仅振荡数据，我们只看到一个突出的k₄模态，该模态非常接近标称值（图7.E）。α的边际分布也呈现出类似的趋势。补充图C.15中的动力学后验分布和补充图C.16中的极限环特性分布表明，这些模型参数的不确定性增加会导致预测动力学的变异性增大。在模型参数不确定性较高的情况下，例如补充图C.15.A所示，x₃(t)的95%可信区间并未严格约束动力学，且极限环特性变化更大（补充图C.16.A）。

最后，我们利用贝叶斯分析的结果来探索参数对与相应预测动力学性质之间的关系。基于观察到k₄和α的不确定性增加会导致错误预测的固定点轨迹比例升高，并结合先前关于模型参数如何相互作用以控制正负反馈之间平衡的知识[56]，我们假设可能会观察到模型参数对与预测动力学之间的相关性。为了验证这一点，我们在补充图C.17、补充图C.18和补充图C.19中绘制了每个数据集的所有参数对的二维散点图，并在图7.F-H中绘制了k₄和α的关系图。如果对应的轨迹产生极限环，则点被着色为蓝色；如果轨迹达到固定点，则着色为红色。我们观察到k₄和α之间最强的相关性，其中对于等距采样数据（见图7.F），两个参数的值都需要高于最小值但又不能过大才能产生振荡。非等距采样数据（图7.G）和仅振荡数据（图7.H）的图表进一步突出了k₄和α之间的关系。根据这些结果，我们得出结论：k₄和α是相关的，它们的值必须处于特定范围内才能产生极限环。我们注意到，我们能够通过贝叶斯分析发现这种相关性；虽然它在等距采样数据中因高估计不确定性而最为明显，但它实际上是模型本身的属性，而非数据的属性。

总之，对MAPK模型的综合不确定性量化（UQ）揭示了多稳态的存在如何为参数估计引入额外的不确定性。具体而言，敏感性分析识别出MAPK模型的两个参数——k₅和α，它们分别仅对双稳态动力学和振荡动力学重要。在这两种情况下，贝叶斯参数估计都能够预测正确的动力学类型，但对动力学的具体特征仍表现出剩余的不确定性。通过改变用于估计的数据，我们发现所提供数据集的性质会影响我们预测极限环的能力。有趣的是，我们发现更多的数据并不总能改善估计效果，而更聚焦于极限环的较少数据反而减少了估计不确定性并给出了更好的后验参数分布。此外，我们发现k₄和α的值需要处于一定范围内才能预测振荡型MAPK动力学。总体而言，所提出的综合UQ框架找到了对每种动力学状态都重要的参数，并直接量化了这些参数的不确定性如何导致动力学的不确定性。

3.3. 长时程增强/抑制现象学模型中的参数估计
文献[57]提出了一种激酶-磷酸酶耦合开关的现象学模型，其中激酶和磷酸酶的活性影响膜结合型AMPAR（α-氨基-3-羟基-5-甲基-4-异噁唑丙酸受体）的水平，以此作为突触可塑性的报告指标，用以捕捉突触可塑性的关键事件。该激酶-磷酸酶模型包含三个状态：pK(t)、P(t) 和 A(t)，分别对应激酶（文献[57]中为CaMKII）、磷酸酶（文献[57]中为PP2A）以及膜结合型AMPAR的活性形式。其微分方程如下：

接下来，我们使用CIUKF-MCMC算法结合AIES来估计缩减参数集的后验分布。参数是从包含36个全状态测量点的含噪合成数据中估计得出的，这些数据对应于LTP响应。数据在0 ≤ t ≤ 9 (秒)的域内均匀分布，噪声水平为真实轨迹方差的1%。由150条马尔可夫链组成的集合每条链运行8,000步，其最大积分自相关时间为621.58，因此我们丢弃了4,351个样本作为预热期。所有参数的马尔可夫链集合轨迹见补充图C.24。图9.A展示了自由模型参数k₂、k₆、P₀、K₀、Ptot、Ktot、Atot的估计边际后验分布。我们观察到不同估计参数之间存在不同程度的不确定性。例如，Atot的边际后验分布显示出非常高的确定性，几乎所有的概率质量都集中在最大后验概率（MAP）点，且该点与标称值对齐。然而，P₀和K₀的边际后验分布表现出更大的不确定性，因为我们观察到后验概率分布在先验支持范围内的广泛扩散。此外，k₂、k₆、Ptot和Ktot的边际后验分布显示出一个远离标称值、围绕MAP点的大模态，以及一个位于标称值处的小模态。

利用后验样本，一组30,000次模拟（见第2.10节）代表了对LTP诱导型输入的预测动力学的后验分布，如图9.B所示。我们观察到，尽管输入是LTP诱导型，但归一化EPSP的95%可信区间同时覆盖了LTP（归一化EPSP > 1）和LTD（归一化EPSP < 1）响应。然而，在MAP点（蓝色虚线）评估的轨迹与真实轨迹（黑色虚线）相匹配，表明大多数轨迹与预期的LTP响应一致。对集合模拟中各单个轨迹（图9.C，顶行）和归一化EPSP稳态值（图9.D，顶行）的检查证实，对LTP诱导型输入既有LTP（蓝色轨迹）也有LTD（黑色轨迹）响应。具体而言，我们发现76.59%（22,972个样本）的响应正确预测了LTP，而23.41%（7,024个样本）的响应错误地预测了LTD（图9.E）。因此，尽管数据质量很高（测量点多且噪声低），我们仍观察到预测的归一化EPSP存在显著的不确定性。

最后，我们研究了在LTP状态下估计的后验分布是否能够预测对LTD诱导型输入的响应。最初的假设是，在响应LTD诱导型输入时，大多数轨迹会预测LTD，但有一部分重要的子集会预测错误的响应——LTP。另外一组30,000次模拟，采用Ca²⁺(t) ≡ 2.2 [μM]（脉冲时间为2 ≤ t ≤ 3秒）的LTD诱导型输入，用于确定动力学的后验分布。对图9.C（底行）中30,000条轨迹中的100条轨迹的可视化再次显示了LTP（蓝色）和LTD（黑色）响应；然而，我们意外地观察到，对LTD诱导型输入的响应中LTP多于LTD。归一化EPSP响应的分布（图9.D底行）证实，大量响应处于LTP状态，且在正确响应附近有一个较小的模态。对这些结果的量化突出表明，仅有31.07%（9,320个样本）的响应处于LTD状态（预期响应），而68.93%（20,680个样本）的响应处于LTP状态（非预期响应）。总之，该模型公式能够正确捕捉一定范围模型参数下的相同LTP行为，但在预测LTD行为方面能力丧失。

尽管LTP和LTD响应对同一组参数敏感（见图8.C-D），但从LTP测量数据估计出的后验分布更倾向于那些能稳健预测LTP而非同时正确预测两种响应的参数集。基于这些结果，我们得出结论：对于该模型，我们需要以高度确定性学习模型参数，才能区分LTP与LTD响应，因为敏感性分析揭示出同一组参数控制着这两种不同的输出。这一发现强调，对于具有多稳态性的系统，仅靠敏感性分析不足以区分参数不确定性，此处概述的综合框架对于揭示此类模型复杂性是必要的。

1. 讨论
   在本研究中，我们开发了一个用于系统生物学动力学模型综合不确定性量化（UQ）的框架（见第2.1节）。所提出的框架结合了可识别性分析和敏感性分析以缩减参数空间（第2.3和2.4节），随后采用CIUKF-MCMC进行贝叶斯参数估计（第2.6节）。我们将该框架应用于三个系统生物学模型，以展示其适用性，并强调关注不确定性如何转变基于建模的研究范式。首先，我们在一个简单的双态模型上进行了两次计算实验，展示了噪声和数据稀疏性如何导致估计不确定性。接着，我们将该框架应用于两个更能代表实际生物读出建模的模型——MAPK模型和突触可塑性模型。利用这些模型，我们阐明了综合不确定性量化如何实现对两种具有生物学意义的动力学行为（极限环振荡和稳态响应）的定量分析。我们还发现，即使数据质量很高，模型结构本身所固有的特性仍可能导致无法消除的不确定性。这些示例为在实践中应用我们的框架以及在不确定性条件下解读系统生物学研究提供了重要的起点。

我们的结果表明，聚焦于不确定性量化能够为基于建模的研究带来新的洞见。例如，在第3.2节中，我们成功获得了能够以高概率预测极限环振荡的参数后验分布。后验分布是在将可用数据纳入统计模型后，对模型参数分布的最佳估计。因此，通过集合模拟（如第2.10节所述）近似得到的动力学后验分布，代表了在已知模型和数据的前提下，对系统动力学的最佳估计。对于MAPK极限环振荡而言，这一最佳估计包含了一系列具有不同极限环特性的动力学行为。对产生极限环振荡的参数进一步分析揭示了k₄与α之间的相关性，并表明它们的取值必须同时落在特定范围内才能产生预期的动力学行为。因此，将不确定性量化纳入建模过程，为理解所研究系统提供了额外的背景信息和更深入的认识。

我们还观察到，在高MAPK稳态（第3.2节）和LTP/LTD响应（第3.3节）的例子中，预测结果并不总能正确捕捉到系统固有的双稳态特性。在这两种情况下，预测动力学集合的95%可信区间同时覆盖了高稳态和低稳态（在突触可塑性例子中即LTP和LTD）。这些结果表明，所估计的模型后验分布中包含了某些参数集，这些参数集在给定特定输入（MAPK模型的初始条件或突触可塑性模型的钙离子水平）下不再表现出双稳态。对这些系统的进一步分析可检验这些参数集究竟是完全丧失了双稳态，还是仅仅使分岔点（即决定系统趋向哪个稳态的输入临界值）发生了偏移。总体而言，我们的结果揭示了模型参数与输入之间复杂的相互作用，这种相互作用可能使多稳态系统的参数估计变得困难。

正如第3.1和3.2节所强调的，我们证明了估计不确定性源于可用实验数据的质量、数量和性质。正如预期，在简单双态模型中（图3），我们展示了噪声更大、样本更少的数据会导致更宽的不确定性区间。此外，对于MAPK振荡（图7），我们发现采样策略——即是否包含初始瞬态阶段的数据，或仅聚焦于振荡阶段的数据——会极大地影响估计不确定性。这一发现与“更多数据总是更好”的普遍直觉相悖，因为我们发现使用样本最少的“仅振荡数据”反而获得了最小的估计不确定性。我们认为，这一现象可从似然函数的角度理解：似然函数通常采用类似L²范数的方式度量数据与估计轨迹之间的不匹配程度（例如见公式(12)）。当同时包含初始瞬态和振荡数据时（如第3.2节中的等距和非等距采样情形），所有能较好拟合初始瞬态的轨迹（无论其后续是否振荡或趋于固定点）都会获得相似的似然概率。而仅使用振荡数据则使似然函数聚焦于长时间动力学行为，并对趋于固定点的轨迹赋予更低的概率。尽管基于L²范数的似然函数在贝叶斯方法中被广泛使用，也可以设计其他形式的似然函数，例如明确判断轨迹是否具备与数据相同特征的基于特征的方法（如文献[90]所述）。这些发现表明，精心设计的实验方案可以减少不必要且可能昂贵的数据采集，并提高参数估计和模型预测的确定性。

在所提出的框架中，我们强调全局可识别性分析和敏感性分析对于实际参数估计和不确定性量化至关重要。通过在补充材料附录A.1中运行不含这些分析的框架变体，我们发现每种分析都能有效降低估计不确定性，并使估计问题更易求解。在完整框架中，我们首先筛选出全局可识别的参数，然后基于Sobol敏感性指数设定一个截断阈值（见第2.4节），以选择有显著影响的参数子集。在本文所有案例中，我们都选择了一个合理的敏感性指数阈值，该阈值通常对应敏感性指数的急剧下降点或一个逻辑上合理的数值，例如我们在多个示例问题中采用的Sᵢ ≥ 0.1标准。然而，我们提醒读者，选择有影响的参数子集是一个复杂的问题，需要用户自行决定在不确定性分析中可以接受舍弃多少信息[91, 92]。除阈值法外，其他方法还包括文献[91]中讨论的方法、组敏感性分析[37]或Lasso回归[93]。尽管这些方法旨在选择有影响的参数子集，但它们与阈值法类似，都依赖于确定可从分析中排除的自由度数量。未来的研究应致力于开发完全定量的参数子集选择方法。

尽管本文三个示例模型的计算仍属可行，但必须认识到，完整的不确定性量化会显著增加参数估计的计算成本。文献[26]指出，原始UKF-MCMC算法所需的浮点运算量远高于那些未考虑所有不确定性来源的方法。类似地，我们发现（见补充材料附录A.2），单次CIUKF似然函数评估的耗时是不考虑模型形式不确定性的类似似然函数的81倍。对CIUKF似然函数执行过程的深入分析表明，约60%的运行时间花费在调用ode15s()上，我们使用该函数对微分方程模型进行离散化。未来可通过开发优化的微分方程离散化和似然函数评估代码，降低CIUKF带来的额外计算开销。

贝叶斯方法（如CIUKF-MCMC）能够实现完整的不确定性量化；然而，MCMC采样及其所需的重复似然函数评估会带来更高的计算负担。特别是，随着模型状态变量和参数数量的增加，模拟成本上升、维度增高，以及后验分布可能具有的复杂几何结构，都会对标准MCMC方法构成挑战。来自更广泛不确定性量化领域的创新方法，结合日益普及的高性能计算资源，使得对具有大量状态变量和参数的模型进行贝叶斯分析成为可能[19, 94, 95]。此外，尽管我们已证明可识别性和敏感性分析能显著降低参数空间维度，但在某些应用场景中，缩减后的空间维度可能仍然较高。许多MCMC采样算法通过利用分布几何信息或并行运行马尔可夫链，能够有效采样高维且可能复杂的后验分布。这些方法包括并行MCMC采样器[96, 97]、哈密顿蒙特卡洛（Hamiltonian Monte Carlo）[27]，以及维度无关、似然信息驱动的方法[98]。尽管综合不确定性量化带来了额外的计算成本，但不确定性量化领域的最新进展可在一定程度上缓解这些负担。

在本文工作中，我们对模型方程的可获得性及相关误差的统计特性做出了若干假设。首先，我们假设在参数估计前模型方程是已知的。这一假设反映了标准建模范式，即模型基于生化反应动力学的理论构建。然而，在使用CIUKF-MCMC时，我们通过引入过程噪声在一定程度上弱化了这一假设，以考虑模型形式的不确定性[26]。实际上，所有模型都因对系统做出的假设而存在一定程度的不确定性。因此，考虑模型形式不确定性可对动力学进行正则化，以弥补预测动力学与数据之间的不匹配[26]。然而，在模型结构事先未知的情况下，可能需要同时从数据中估计模型结构和参数。解决模型结构模糊性的一种方法是从数据中直接学习生化反应网络[99, 100]或数学模型本身[101–105]。此外，也可将此类问题置于贝叶斯框架下，同时学习模型形式及其相关不确定性[26, 106]。

此外，在选择测量噪声和过程噪声的统计模型以及模拟生物测量数据时，我们还做出了一些假设。首先，在应用CIUKF-MCMC方法时，我们假设测量噪声和过程噪声均为高斯分布。标准卡尔曼滤波方法通常基于噪声过程服从正态分布的假设[67]。然而，生物系统和实验中的噪声或许更适合用其他概率分布来描述（参见，例如[46]）。若假设测量噪声服从非高斯分布，则需重新构建似然函数（即修改公式(12)中的分布）；而若对过程噪声采用替代分布，则需采用超越卡尔曼滤波的近似方法来计算边缘似然。其次，本文考虑的是线性测量函数（如公式(3)），但CIUKF-MCMC同样适用于非线性测量函数[26]。最后，我们假设仅能获得单条时间序列的测量数据（即单次实验）；然而，大多数生物学实验都会进行多次重复试验。一种直接的方法是使用每个时间点的均值构建时间序列，并基于该均值序列估计参数。此外，也可分别对每条时间序列进行参数估计，然后通过元分析或信息融合原理[107–109]整合多个后验分布，或构建一个能同时处理多条时间序列的统计模型[27]。

本文提出的框架可直接应用于并扩展至系统生物学中遇到的大多数动力学模型，以实现综合不确定性量化。未来的研究应聚焦于：开发严格的参数子集选择方法；将适用于复杂高维分布的MCMC采样器应用于系统生物学；以及整合更符合生物数据特性的统计模型和噪声模型。总之，在不确定性量化与系统生物学建模的交叉领域开展研究，将强化参数估计能力，并推动模型更准确地预测实验观测到的动力学行为。

原文链接：https://arxiv.org/pdf/2204.05415