基于贝叶斯配置积分方法的常微分方程系统识别及其在基因调控网络中的应用

《Biometrics》：A Bayesian collocation integral method for system identification of ordinary differential equations

本文推荐一种贝叶斯层次配置方法，用于从噪声时间序列数据中识别高维稀疏常微分方程(ODE)系统结构。研究人员通过结合积分形式的ODE约束、数据似然和组间稀疏惩罚，实现了系统识别和轨迹估计的同步进行。该方法采用尖峰-板片先验进行变量选择，通过模拟研究和T细胞基因表达数据验证，在轨迹估计和不确定性量化方面优于现有频率学方法，为基因调控网络(GRN)分析提供了可靠的统计推断工具。

　　

在系统生物学和计算生物学的快速发展中，科学家们迫切需要有效工具来揭示复杂生物系统的动态行为。基因调控网络(GRN)作为理解细胞功能分子机制的关键，其动态特性通常通过常微分方程(ODE)系统进行建模。然而，从噪声时间序列数据中准确识别高维稀疏ODE系统结构面临诸多挑战：观测数据稀疏且带有测量误差，参数估计的不确定性难以量化，现有频率学方法在轨迹置信区间构建方面存在局限。

传统方法多采用两步配置法或基于导数的估计技术，但这些方法对导数估计的准确性高度敏感，特别是在数据稀疏情况下表现不佳。虽然近年来出现了如GRADE（通过加性微分方程的图重构）和JADE（广义稀疏加性ODE的联合估计）等改进方法，但它们仍局限于频率学框架，无法充分量化参数估计的不确定性。

针对这些问题，发表在《Biometrics》上的研究提出了一种创新的贝叶斯层次配置方法。该方法将积分估计技术与贝叶斯推断相结合，通过统一处理数据似然、ODE约束和稀疏惩罚，实现了系统识别和轨迹估计的同步优化。研究团队来自滑铁卢大学统计与精算科学系，他们开发的方法在保持计算效率的同时，提供了可靠的不确定性量化。

研究方法的核心技术包括：采用三次B样条基函数展开逼近ODE解和加性分量函数；利用高斯求积法则数值近似积分项；设计自适应的组间稀疏惩罚（尖峰-板片先验）进行变量选择；基于哈密尔顿蒙特卡洛（NUTS）算法的后验采样。这些技术方法的有机结合，使得该方法在保持计算效率的同时，提供了可靠的不确定性量化。

模拟研究验证性能

通过包含10个变量的ODE系统模拟实验，研究人员在不同信噪比(SNR)和观测密度条件下评估方法性能。当SNR=25时，该方法在轨迹均方误差(MSE)方面与JADE方法相当，但明显优于GRADE方法。随着噪声水平增加（SNR降至4），该方法的优势更加明显，轨迹MSE显著低于对比方法。在加性分量估计方面，该方法在所有噪声水平下对非零分量和零分量的估计误差均最低，且表现出更好的稳定性。网络识别方面，虽然GRADE的真阳性率(TPR)略高，但该方法的假阳性率(FPR)控制更好，在保持较高识别准确性的同时减少了错误发现。

不确定性量化优势

特别值得关注的是该方法在轨迹不确定性量化方面的突出表现。通过计算95%点覆盖概率，研究发现该方法构建的置信区间能够接近名义水平，而频率学方法往往存在覆盖不足的问题。即使在观测稀疏（21个时间点）情况下，该方法仍能保持可靠的统计推断能力，为科学家评估结果可靠性提供了重要依据。

T细胞基因表达数据应用

在真实数据应用中，研究人员分析了人类T细胞对佛波醇肉豆蔻酸乙酸酯(PMA)和离子霉素处理的响应数据。该数据集包含58个T细胞活化基因在10个不等间距时间点的表达测量，共34个重复样本。该方法成功构建了非线性动态GRN，识别出273个非零调控关系（总共3364个可能函数），证实了网络的稀疏特性。

分析结果与已有研究发现高度一致，如FYB（基因45）被确认为关键调控因子，具有最多的外向连接。除了已知的5个FYB调控基因外，该方法还新发现了SIVA（基因6）、IL-16（基因23）、MyD88（基因37）等5个新的调控靶点。这些新发现的调控关系得到了文献支持，如SKAP2与FYB的直接相互作用，以及SKAP2与MyD88的物理关联都为新发现提供了佐证。

图3展示了FYB、IL-2Rγ、凋亡相关基因和IL3Rα四个基因的表达轨迹估计结果。可以看到，尽管原始数据存在较大噪声，该方法仍能提供合理的轨迹拟合，且点wise可信区间能够反映不同时间点的观测密度差异——在观测稀疏的时间区间，置信带自动变宽以体现更大的不确定性。

图4展示了58个基因构成的GRN网络，其中FYB节点的连接用粗箭头突出显示。网络结构清晰显示了基因间的复杂调控关系，每个基因仅受少数其他基因调控，符合生物系统的稀疏特性。

图5进一步揭示了FYB调控作用的非线性特征。以caspase 7（基因52）为例，在低FYB表达水平时调控效应增强，而随着FYB表达水平升高转为负向调控。这种动态变化的调控模式凸显了该方法在捕捉非线性效应方面的优势。

研究还通过敏感性分析验证了方法的稳健性。在不同先验参数设置下，非零函数的选择结果和计算时间保持稳定，表明方法对先验选择不敏感。此外，该方法计算时间随变量数增加呈近似线性增长，具有良好的可扩展性。

该研究的创新性主要体现在三个方面：首先，将积分估计技术引入贝叶斯框架，避免了导数估计的不稳定性；其次，采用尖峰-板片先验实现自动变量选择，提高了模型识别准确性；最后，基于后验采样的不确定性量化为轨迹估计提供了可靠推断工具。

未来研究方向包括将方法扩展到高阶ODE系统，通过格林匹配技术处理更复杂的动态模型。此外，该方法在传染病动力学、神经科学等领域的应用也值得进一步探索。

这项研究为高维动态系统建模提供了新的统计工具，特别是在基因调控网络分析中展现出重要应用价值。通过可靠的系统识别和不确定性量化，该方法有助于科学家从噪声数据中提取有意义的生物学洞见，推动精准医学和系统生物学的发展。