馬爾可夫基因表達建模的神經網絡矩閉合方法

關鍵詞：基因表達建模；神經網絡；矩閉合方法；隨機模擬；最大熵原理

基因表達作為理解生物學現象的核心焦點之一，在生命科學領域的快速發展中占據著重要地位。基因表達是基因通過轉錄和翻譯實現基因功能產物合成的過程，是生物體內調控和執行生命過程的關鍵步驟，通過細胞調控和執行基因的功能維持生物體的正常功能。了解基因表達機制不僅有助于理解生命的本質，還對揭示疾病發生、發展以及藥物研發具有重要意義[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結束語

在基因調控網絡建模過程中，隨機模擬算法在獲取矩值時需進行大量的隨機模擬并集合平均，導致計算量龐大和復雜性增加。而依賴于簡化假設的傳統矩閉合方法則無法充分描繪真實系統的復雜性，不能有效捕捉大量相互作用的生化反應模型系統的物理細節。因此，本文提出了一種新穎的神經網絡矩閉合方法，它通過在整個生化反應網絡中探索潛在關聯，能夠更全面地捕捉生化反應模型中的動態行為。實驗證明，相較于傳統方法，神經網絡矩閉合方法在對基因表達模型的預測精度和時間效率上都表現出一定的優勢，為基因表達建模研究提供了一種更準確和高效的分析工具。

盡管神經網絡矩閉合方法在生化反應建模方面取得了顯著的進展，但也存在著挑戰和改進的空間。本文的實驗驗證主要局限于特定的基因調控網絡模型，因此該方法在遇到未知情境時的泛化能力可能不足。此外，盡管本文在方法驗證階段使用的是模擬數據集，但與實際生物實驗數據的結合是提升方法可靠性和應用實用性的關鍵。未來的研究應當著重于將神經網絡矩閉合方法應用于更為廣泛的生化反應模型，并提升模型可解釋性，以改善用戶對預測決策的理解。同時，與更多的反應類型的結合也將是增強方法魯棒性和驗證可行性的重要步驟。總而言之，通過解決現有問題并成功地將研究前景轉化為實際成果，神經網絡矩閉合方法有望在生化反應建模領域實現更廣泛的應用。