马尔可夫基因表达建模的神经网络矩闭合方法

关键词：基因表达建模；神经网络；矩闭合方法；随机模拟；最大熵原理

基因表达作为理解生物学现象的核心焦点之一，在生命科学领域的快速发展中占据着重要地位。基因表达是基因通过转录和翻译实现基因功能产物合成的过程，是生物体内调控和执行生命过程的关键步骤，通过细胞调控和执行基因的功能维持生物体的正常功能。了解基因表达机制不仅有助于理解生命的本质，还对揭示疾病发生、发展以及药物研发具有重要意义[1-3]。在生物学中概率主宰生物学，概率在噪声塑造生物系统行为方面起着至关重要的作用[4-7]。这里所述的“噪音”源自于活细胞内分子浓度的固有波动，主要是由生化反应的随机性引起的，尤其在低分子数量的生化反应中表现得尤为明显。因此，对于生化反应网络的建模主要倾向于对单个反应进行模拟来表现反应中分子数的随机波动[8]。基于这一见解，以化学主方程（ChemicalMasterEquation，CME）为基础的马尔可夫模型等低分子随机建模方法迅速流行[9-10]。同时，随机模拟算法（StochasticSimulationAlgorithm，SSA）也成为了解和获取基因表达动态过程的重要工具[10]。然而SSA的计算成本很高，适用性受到严重限制，难以应用于大型系统。

矩闭合近似方法（MomentClosureApproximations，MMA）在研究基因表达网络的稳态和极限行为方面取得了重要成就。大多数矩闭合方法主要用于估计分布的矩[11-16]，从而得到关于所有阶及以下联合分布矩的时间演化近似解[11-13]。进一步可以利用最终稳态时刻的近似矩值，使用最大熵原理重构相应的边缘概率分布[17-18]。通过分析系统的矩集合，这类方法能够从全局角度理解基因调控网络的整体行为。然而，传统矩闭合方法仍存在一些挑战和局限性。首先，传统矩闭合方法通常基于线性稳态分析，其在非线性系统中的适用性受限，而许多基因调控网络是非线性的，导致传统矩闭合方法无法充分捕捉基因表达网络潜在过程相互作用的复杂性。其次，传统方法往往局限于特定类型的生化反应网络，难以灵活处理不同生物体和细胞类型之间基因调控机制的差异。此外，在某些复杂生化反应网络中，即使采用更高阶的矩闭合方案，其精度也可能受到闭合方案和生化反应网络潜在物理过程复杂性的限制。总的来说，传统矩闭合方法在适用范围和精度方面存在不足。

人工神经网络（ArtificialNeuralNetwork，ANN）通常都是对自然界某种算法或者函数的逼近，也可能是对一种逻辑策略的表达[19]，近年来人工神经网络与其他学科领域联系日益紧密，在各个领域得到广泛应用，通过对神经网络层结构的探索和改进来解决不同领域的问题[20]。受此启发，本文提出了一种基于神经网络的矩闭合方法，称为神经网络矩闭合（Neuralnetworkmomentclosure）方法。该方法利用人工神经网络学习基因调控网络模型的矩方程组中高阶矩的低阶表示，将未闭合的矩方程组闭合，再通过线性常微分方程组求解获得估计的矩值。与传统矩闭合方法相比，神经网络矩闭合方法无需对系统进行额外分布假设，更能充分利用生化反应网络模型中的未知潜在特性，捕捉背后复杂的物理相互作用。一旦神经网络学习到这种潜在相互作用，说明所提出的方法能够学习到生化反应模型中的物理行为，使矩闭合结果更加真实可信和准确。神经网络矩闭合方法不仅提供了一种获取矩闭合方法的新途径，而且弥补了传统方法在生化反应网络系统模型近似中的不足。本研究有望推动基因表达建模领域的发展，为深入理解基因调控网络的动态行为提供新的视角和方法。

1预备知识

1.1随机模拟算法

CME所描述的随机过程本质上是一个连续时间马尔可夫过程，其中连续反应事件之间的时间间隔服从指数分布[21]。由于从指数分布中抽样相对简单，因此模拟生化反应的发生非常便捷且直接。SSA算法基于概率分布的数值抽样，可以模拟底层随机过程的精确样本路径，从而提取准确样本，是一种在状态空间中生成随机轨迹集合的动力学蒙特卡罗方法。这使得SSA能够在分子层面上捕获化学反应的随机性质，提供精确的分子轨迹，并且适用于广泛的化学反应网络。

假设一个生化反应网络系统是由个不同的化学反应物和个分别对应反应通道的反应组成。每个反应都有一个倾向函数，反应系统状态用表示，表示反应物在时刻的分子数，表示向量的转置。直接随机模拟算法的模拟过程如下：首先对将要发生的反应所需的时间间隔步长进行采样，然后对反应集合中的某个具体反应进行采样，从而确定是哪个反应在什么时间完成[22]。具体而言，表示下一个反应在时发生的概率，并且该反应在一个无限小的时间间隔内完成；表示下一个反应是反应的概率。这两个概率可以通过相应的计算公式从fr（n）dt获得，如式（1）、（2）所示：

其中，u1和u2为0到1之间的均匀随机数，SIS代表满足公式的最小整数。直接法首先根据式（3）对下一个反应事件的时间点进行采样，然后根据式（4）对发生某一反应进行采样，迭代更新随机模拟过程的状态向量和时间。

由于随机模拟算法模拟系统中的每个化学反应事件都是明确的，即使对于反应物种类较少的系统，随机模拟算法的计算成本也很高。这种高计算成本的情况在分子数波动很大或单位时间内发生大量反应的情况下尤为明显。在第1种情况下，为了获得统计上准确的结果，必须模拟大量样本。而在第2种情况下，由于反应事件之间的时间变得更短，单次模拟的计算成本也变得昂贵。因此，随机模拟算法的适用性受到严重限制，并且很快就无法适用于大型系统。为了克服这些挑战，近几十年来，研究人员投入了大量精力来发展化学主方程的近似方法，并出现了多种不同的方法。其中一种称为Tau跳跃的方法（Tau-leaping）是一种模拟生化反应的近似方法，它的主要目标是提供比SSA更高效的性能[23]。该方法的核心理念在于通过时间上的离散“跳跃”，跨越多个反应事件，从而避免了对每个单独反应事件进行模拟的需要。这允许系统在有限的时间段内经历多个反应，大幅度减少了必须处理的事件总数，加快了模拟的速度。除了Tau跳跃，还有其他近似方法被提出来，这些方法的共同目标是高效地近似CME的解，以此降低计算的复杂性和成本。

1.2近似方法

CME有很多近似方法，其中3种最常见的近似方法分别是化学朗之万方程（ChemicalLangevinEquation，CLE）、系统尺寸展开（SystemSizeExpansion，SSE）和MA[24-25]。这3种方法易于实施，无需对系统有任何预先的了解，而且它们通常能够进行高效计算，并提供精确近似。因此，它们已被成功应用于各种场合[26-30]。然而，这些方法在某些情况下的准确性可能大幅下降，尤其是当某些物种的拷贝数非常低时。如果关注的是过程的矩，CLE通常比SSE和MA更为准确。但是，CLE在计算上的代价更高，因为它需要进行大量的随机模拟并集中平均来获取过程的矩，而其他方法只需求解一组有限的常微分方程。此外，当CLE定义为实值变量时，在零分子数处会遇到边界问题，实值修正又会引入新的不准确性[31]。通过将CLE扩展到复值变量可以解决边界问题，但会降低模拟的效率[32]。因此，如果只对过程的矩感兴趣，使用系统大小扩展或矩闭合近似似乎是更合适的选择。

另一方面，系统尺寸展开是基于小参数的系统扩展，而矩闭合近似是一种特定的近似方法。系统尺寸展开在大系统容量下可以保证准确性，因此在大规模系统下它更具吸引力。对于矩闭合近似，通常不期望能够在所有情况下保持同样的准确度。另外，系统大小扩展不适用于某些确定性具有多稳态的系统，这是矩闭合方法不具有的限制[33]。更进一步地，系统大小扩展仅在均值上高于线性噪声近似两个阶，在协方差上高一个阶[34]，系统大小扩展的高阶矩修正比矩闭合方法更难以推导和实现；而矩闭合近似则可以推广到各种阶数[35-36]。CLE、系统大小扩展和矩闭合近似通常作为基础构建模块，为开发高级建模策略提供了框架。比如，有限状态投影算法（FiniteStateProjectionAlgorithm，FSP）的思想是将状态空间截断为有限子空间，并使用矩阵幂运算求出该子空间上分布的近似值[37]。鉴于这些因素，选择哪种方法更为合适，将取决于具体问题的细节。

在对比CLE和SSE的基础上，本文选择聚焦于MA中的矩闭合技术。矩闭合方法在操作性上提供了广泛的灵活性，近年来，多领域的专家和学者在人工智能技术的研究和应用中取得了突破性进展[38]。

对于线性反应系统，CME方程可以在一定的期望阶数上进行数值求解。然而，对于非线性系统，低阶与高阶方程相互耦合，导致矩方程的无限耦合层次，因此不能直接求解。矩闭合方法通过一种近似的方式截断了这个无限阶方程组，常用的矩闭合近似就是通过将所有高于阶的矩表示为低阶矩的函数来闭合矩方程。为了实现这个目标，一种方法是假设系统分布具有特定的函数形式，比如正态分布。这样的假设将阶矩方程与高阶矩解耦，从而得到一组有限的解耦合的常微分方程组。数值求解这组闭合的方程就可以获得所需的矩估计值。这样的矩闭合方法称为“M阶矩闭合”。

2神经网络获取矩闭合方法过程

本文提出的神经网络矩闭合方法的核心是假设有限数量的矩能够捕捉到所有必要的系统信息，通过神经网络学习到生化反应系统未闭合的矩方程组中高阶矩的低阶矩表示函数，就可以将矩方程组闭合，随后通过解闭合的微分方程组来获取矩估计值。

图1示出了整个实验流程。实验首先要构造所需的特定生化反应模型和输入数据集。虽然流程图中描绘的是一个基因调控网络模型，但方法同样适用于构建更广泛类型的生化反应模型。针对研究需要的生化反应模型，需要生成大量的随机参数组作为模型的输入，其中每个参数组代表生化反应模型的不同倾向函数的反应过程。为了让神经网络能够捕获生化反应模型的底层特性，需要足够数量具有广泛性和代表性的倾向函数随机参数集。这些参数集的数量和范围可能需要根据实验结果进行进一步调整。

利用生成的有效数据集，一方面，需要获取生化反应系统的原始未闭合矩方程组（Rawmomentequations）。这些方程组将运用神经网络学习到的矩闭合方案并求解闭合方程。另一方面，运用SSA随机模拟并进行集中平均，以获得生化反应模型的矩真实值，此值将作为神经网络训练数据集的参考真值。神经网络的输出是高阶矩的低阶矩表示，为了实现这一点，需要针对不同生化反应网络构造不同的向量表达方式。将神经网络的输出代入到原始矩方程组中，即可成功实现方程组的闭合，这为常微分方程组的求解提供了便利，进而获得了矩的估计值。通过将求解得到的矩估计值与SSA得到的矩真实值进行比较，得到模型的偏差，利用偏差对神经网络进行反向传播更新梯度值，直至满足预期的性能指标。

神经网络的训练过程遵循标准的训练算法，如算法2.1所示。

算法2.1 神经网络训练算法

1 加载数据集并归一化处理；

2 设置学习率 =0：1，正则化系数；

3 随机初始化神经网络权重和偏差W;b；

4 repeat

5 训练集样本进行随机排序；

6 forn2trainsetdo

7 正向传播得到神经网络输出v（i）；

8 闭合矩方程组，使用常微分方程求解得到估计矩值，并求出目标函数；

9 反向传播，计算每一层的误差和导数；

10 更新网络参数；

11 endfor

12 until神经网络在测试集上的错误率不再下降

13 输出神经网络模型的参数W;b

值得注意的是，经过一轮训练后，根据神经网络学习到的矩闭合效果，可能需要对参数进行调整，或者对网络结构进行优化，以实现更精确的估算结果。

3实验结果分析

3.1基因调控网络模型及数据集介绍

3.1.1基因调控网络模型 本文实验对象采用的是生化反应中极具代表性的基因调控网络（GeneRegulatoryNetwork，GRN）模型。这种反应网络模型是一个用于描述细胞内或一个特定基因组内基因间相互作用的抽象模型，在众多相互作用关系之中，侧重于基因调控机制的相互作用。基因调控网络是生物体内控制基因表达的关键机制，它涉及基因的转录和信使核糖核酸（mRNA）的翻译过程。图2示出了GRN模型示意图[43]。

3.1.2基因调控网络数据集 为了实施图1所描述的基因调控网络模型的神经网络学习矩闭合方法，需要构建数据集，其中是神经网络的输入，即反应方程组的倾向函数组成的向量，是模型经过30000次SSA随机模拟并进行集中平均得到的精确矩值。由于本模型是双变量，所以用分别代表基因和蛋白质阶和阶时的矩值，针对本文的模型将表示成。数据集大小M设置为4000个，然后按照9∶1划分为训练集和测试集。

图3中展现的趋势和分布情况不仅揭示了蛋白质数量随时间的动态演变，而且也体现了在达到稳态时各个状态的概率分布。通过分析，可以确认数据集中的矩闭合值是在稳态条件下计算的，这一点对于验证数据集的精确性至关重要。此外，还可以观察到数据集具有广泛的代表性，这种特性对于保障数据集在模拟各类生化反应网络时的通用性和适用性极为关键，确保了模拟实验结果的稳定性和可重复性。通过选取覆盖多种可能情境的不同参数组合，确保数据集能够覆盖大范围的数据空间，这进一步证明所选数据集在适用性和可靠性方面的优势。

需要注意的是，本文所采用的基因调控网络模型，虽然是一种简化的抽象表达形式，它对于理解更为复杂的生化反应系统的动态行为提供了初始的框架。然而，对于那些对高度复杂生物过程的建模感兴趣的研究者来说，使用生成的模拟数据集之前，对其可信度进行细致的评估是必不可少的。为了确保所生成的模拟数据集能够准确地反映真实世界的数据特性，需要使用一系列细致的量化指标和对比分析方法：

（1）统计一致性：包括对模拟数据集与真实数据集的平均值、中位数、方差等核心描述性统计指标进行比较，并利用Kolmogorov-Smirnov检验和Q-Q图等方式来详细对比数据分布的相似度；（2）时间序列分析：分析模拟数据集和真实数据集分子数量随时间变化的行为模式，确保模拟数据能够精确地再现真实生物系统的动态特性；（3）再现性测试：对于每组参数多次运行模拟过程，并检查结果的再现性和变异性，有助于验证模拟过程的稳定性。

在实际实施中，需要充分考虑到研究目的的具体性和所使用数据集的独特性质，以便选取最适合的评估工具和方法。

3.2神经网络训练结果

本文构建的人工神经网络旨在学习基因调控网络模型中的内在反应特性，因此神经网络设计相对灵活，允许多种修改和实验，只要能够有效捕捉生化反应模型的关键特征即可。具体而言，针对本文的研究对象所构建的神经网络包括：（1）一个由4个神经元组成的输入层；（2）两个隐藏层，每层各含10个神经元；（3）包含7个神经元的输出层。网络中输入层与隐藏层之间采用ReLU函数作为激活函数。在训练过程中，采用ADAM优化器推荐的标准对学习率进行自适应调整。针对不同的反应网络需要构建不同的神经网络的输出层，如下所示：

神经网络的训练使用标准反向传播算法来进行权重更新和训练。为了衡量训练的有效性，本文追踪了损失函数的变化，并通过训练周期的演进来评估模型性能（图4）。如图4所示，损失函数在训练初期迅速下降，表明模型从初始状态迅速学习并调整参数以最小化损失。随着训练的深入，损失函数下降的速度减慢，并最终趋于稳定。定义成功的收敛标准为，若损失函数在连续20个训练周期内保持在一个特定的范围内波动，便认为模型已经稳定学习到了数据的特征。在本实验中，损失函数在后续30个周期内保持稳定，由此可以判断模型已经成功收敛。

3.3结果准确性

由于本文实验采用的基因调控网络模型最终得出6个矩估计值，因此评估结果也集中在这6个矩值上。图5示出了估计值的不同方法箱型图。图5中的箱型图对比了基于神经网络的矩闭合方法、SSA和传统矩闭合方法在所考虑的基因调控网络模型中的准确度表现。图中的SSA方法表示模型经过2000次SSA随机模拟到达稳态后计算出的三阶矩以下矩值，低数量模拟的SSA方法由于其固有的随机性，准确度会受到部分限制。图中的“Normal”和“DM”分布代表传统矩闭合方法，分别对应于第1.2节中的正态分布矩闭合方法和微分匹配矩闭合方法。

从图5中的结果来看，神经网络矩闭合方法在准确性方面明显超越了低数量SSA模拟计算得到的矩估计值。尽管这是基于较少数量的随机模拟得出的结论，但依然能展示神经网络矩闭合方法的相对准确性，从侧面说明了SSA方法在获得精确的矩估计值时需要进行大量的计算平均，而这正是矩闭合方法的价值所在，它显著减少了计算量的需求。从图中还可以看到，神经网络矩闭合方法在、、、这几个矩估计值上表现得优于传统矩闭合方法，直接证明了神经网络矩闭合方法在准确度方面相比于传统的矩闭合方法在基因调控网络模型具有显著优势。

R2是一个统计指标，用于衡量观测数据与拟合模型之间的吻合程度，取值范围从0到1，越接近1表示模型与观测数据的拟合度越高。图6示出了神经网络矩闭合方法得到的矩估计值的拟合图，突显了这些矩值之间的高度相关性，以进一步验证本文方法在基因调控网络模型中的可靠性。从图中可以清晰地看出，每个矩值的拟合值都接近1，表明神经网络矩闭合方法能够有效地捕捉到这些矩之间的紧密关联，进一步说明了神经网络矩闭合方法在揭示基因调控网络模型中生化反应动态内在规律性的能力。

神经网络矩闭合方法在灵活性上优于传统矩闭合技术，特别是在满足精度要求的可调整性方面。研究者不仅可以针对整体模型精度进行优化，还能够对特定参数进行细致的调校，这一切均通过修改训练阶段目标函数（参考式（10））中的权重实现，或者可以在目标函数中添加额外感兴趣的项以进一步细化。

3.4结果快速性

表1所示为神经网络矩闭合方法与其他一些算法单次获得矩闭合估计值所需的平均计算时间对比结果。具体来说，对于数据集中一组数据，SSA方法和Tau-leaping方法的时间消耗包括了随机模拟过程和集合平均获取矩值；传统矩闭合方法时间消耗包括获取矩方程组、利用传统公式闭合矩方程组和求解闭合方程组获得矩估计；FSP方法包括计算系统的概率密度向量和计算矩值；而神经网络矩闭合方法的时间消耗则包括获取矩方程组、训练神经网络、利用神经网络输出闭合矩方程组合求解闭合方程组获得矩估计。平均计算时间基于本文4000组参数的数据集得出，该时间反映了求得最终矩估计值所需的平均时长。SSA方法，使用的是3.2节中选择的10000次模拟并作为真值的数据。Tau-leaping方法和SSA相同，也是进行了10000次模拟并集合平均。对于传统矩闭合方法，表中平均计算时间为正态分布矩闭合和微分匹配矩闭合两种方法的平均计算时间。

由结果清楚地显示，神经网络矩闭合方法在计算速度上明显优于SSA方法，并且随着生化反应模型复杂性的提升和模拟规模的扩大，这种速度优势将非常显著。与评估中的其他3种方法相比，神经网络矩闭合方法同样展现出了速度上的优越性。这强调了在进行复杂生化反应模拟时，利用神经网络进行矩闭合近似作为提高计算效率的有力工具，尤其在传统算法难以承受高计算负荷时更显其价值。图中神经网络矩闭合方法虽然在表中仅展示了整体的平均计算速度，但神经网络矩闭合方法中最耗时的环节预计为网络训练过程。后续分析将进一步探究数据量的增加对神经网络训练时间的影响。

图7示出了随着数据集样本量的增加，SSA、传统矩闭合方法和神经网络矩闭合方法在获得矩闭合估计值时所需的平均计算时间的变化。对于SSA和传统矩闭合方法，由于它们在获取矩值时采用了固定的实现途径，因此这两种方法的平均计算时间保持不变，不受数据集规模影响。这一点可以从图中的黑色虚线和浅灰色虚线观察得到。神经网络矩闭合方法的平均计算时间随着数据集样本量的增加而提升，这是因为数据集规模的扩大导致了更长的网络训练时间。值得强调的是，在数据集样本量为1000时，神经网络矩闭合方法已能达到SSA在进行30000次随机模拟后的集合平均矩值精度。从图中可以明显看出，SSA所需的计算时间大约是神经网络矩闭合方法的6倍，而传统矩闭合方法所需时间则约为神经网络方法的两倍半。因此，相较于SSA和传统矩闭合方法，神经网络矩闭合方法在计算效率上具有显著优势。

这种计算效率的显著提升主要归功于神经网络矩闭合方法继承并强化了传统矩闭合方法在近似建模领域的优势，同时规避了SSA在执行大规模随机模拟并集合平均过程中所固有的高计算需求。随着生化反应系统规模的扩张，SSA的计算负担将急剧增加，而矩闭合方法所需的计算资源几乎不受影响。此外，矩闭合技术在求解微分方程组时能够运用先进的时间步长优化技术，根据反应动力学的实际特性动态调整求解步长，由此节约了不必要的计算资源。最关键的是，神经网络矩闭合方法通过神经网络的学习能力，实现了对高阶矩方程组中高阶矩的低阶近似表达，在大量模拟的情况下有效避免了直接计算复杂高阶矩的需求。如果研究者需要对时间效率有极端的要求，迫切需要快速执行大规模模拟时，可以牺牲精度提升时间效率。通过选用较小的数据集合或限制迭代次数，可以大幅缩短神经网络训练所需的时间。尽管这样做可能会影响结果的精细度，但在特定的实验环境中，这种方法仍能有效地满足对快速处理的需求。

4结束语

在基因调控网络建模过程中，随机模拟算法在获取矩值时需进行大量的随机模拟并集合平均，导致计算量庞大和复杂性增加。而依赖于简化假设的传统矩闭合方法则无法充分描绘真实系统的复杂性，不能有效捕捉大量相互作用的生化反应模型系统的物理细节。因此，本文提出了一种新颖的神经网络矩闭合方法，它通过在整个生化反应网络中探索潜在关联，能够更全面地捕捉生化反应模型中的动态行为。实验证明，相较于传统方法，神经网络矩闭合方法在对基因表达模型的预测精度和时间效率上都表现出一定的优势，为基因表达建模研究提供了一种更准确和高效的分析工具。

尽管神经网络矩闭合方法在生化反应建模方面取得了显著的进展，但也存在着挑战和改进的空间。本文的实验验证主要局限于特定的基因调控网络模型，因此该方法在遇到未知情境时的泛化能力可能不足。此外，尽管本文在方法验证阶段使用的是模拟数据集，但与实际生物实验数据的结合是提升方法可靠性和应用实用性的关键。未来的研究应当着重于将神经网络矩闭合方法应用于更为广泛的生化反应模型，并提升模型可解释性，以改善用户对预测决策的理解。同时，与更多的反应类型的结合也将是增强方法鲁棒性和验证可行性的重要步骤。总而言之，通过解决现有问题并成功地将研究前景转化为实际成果，神经网络矩闭合方法有望在生化反应建模领域实现更广泛的应用。