枯草芽孢桿菌中色噪聲誘導的ComK基因表達動力學

陳雷艷，劉培江，王浩華,3

（1. 海南大學 理學院，海口 570228; 2. 廣東財經大學 統計與數學學院，廣州 510320;3. 海南大學 熱帶特色林木花卉遺傳與種質創新教育部重點實驗室，海口 570228）

復雜生命過程中包含許多的生物化學反應，這些化學反應都具有潛在的隨機性，在物理學中，通常用噪聲來刻畫這些隨機性。經典的隨機動力系統往往用高斯白噪聲來刻畫這些隨機性。事實上，雖然相比系統自身的時間尺度而言，噪聲的自關聯時間很小，但是其始終恒大于零，這對系統的相變具有本質性的影響。因此，考慮噪聲源的有限關聯時間誘導系統的非平衡相變是合理的，考察ComK基因表達系統的動力學瞬態性質和穩態性質具有一定的探索和研究意義。

ComK基因作為枯草芽孢桿菌種群中的看家基因，包含自激活的正反饋環以及一個通過ComS基因介導的負反饋環，耦合調節ComK基因的表達。然而，正負反饋的耦合環路時刻受到外界噪聲的刺激，形成內外共振的表達模式，誘導ComK基因表達過程具有自關聯時間。據報道，噪聲關聯強度能夠誘導基因表達系統的兩種穩定狀態在時間延遲下發生切換，增強其穩定性[1]；噪聲和交叉關聯時間對隨機系統的平穩特性有不同的影響[2]；高斯有色噪聲及交叉關聯時間等隨機效應能影響腫瘤免疫反應對化療系統的穩定性，提高治療率[3]。文獻[4]研究了由乘性色噪聲和加性色噪聲共同驅動下具有強Allee效應和弱Allee效應的種群模型，對穩定和不穩定狀態進行了生物學分析和解釋。目前，相關文獻均僅涉及色噪聲及其交叉關聯時間，而對色噪聲和自關聯時間協同作用下對系統的動力學特性方面的研究尚未見報道。因此，本研究觀察ComK基因調控系統在色噪聲強度和自關聯時間誘導下的非平衡相變現象，著重考察色噪聲誘導下基因調節系統的瞬態和穩態變化等的動力學性質，旨在揭示不同源噪聲的時間關聯性對系統穩態概率分布和平均首通時間的影響。

1 基因轉錄調節系統

在枯草芽孢桿菌（Bacillus subtilis）中，轉錄因子ComK會激活編碼DNA攝取和重組系統的基因，會激活自己的正反饋環，多蛋白集合物MecA會降解ComK基因的表達水平,外界刺激會激活ComS縮氨酸競爭的抑制ComK的降解，而ComK的過量表達會抑制ComS的表達，因此，形成正負反饋環路[5?6]。基于此，KARMAKAR和 BOSE[7]提出了ComK基因表達形成ComK蛋白的基因表達模型，該模型中包含ComK轉錄因子及其復合物、ComK和ComS的相互反饋作用和非線性作用力，其基本機制如圖1所示。

圖1 枯草芽孢桿菌調控環路示意圖Fig. 1 Schematic diagram of Bacillus subtilis regulation loop

根據上述模型所示的生化反應，x和S可以表示為：

式中，a為 ComK的蛋白質合成速率，b為其降解速率，k0為 自激活速度，k1為ComS激活ComK的速率，n為hill系數，x表示ComK濃度，s表示ComS壓制ComK的速率，S表示ComS濃度。

考慮到ComS激酶相對ComK具有較快的時間尺度，因此，令=0，反帶入方程（1），則ComK的蛋白質濃度x隨時間演化的偏微分方程可表示為[8]：

圖2 隨蛋白質濃度 x(t)變化的函數Fig. 2 as a function of protein concentration x(t)

圖3 V(x)隨 蛋白質濃度x(t)變化的函數Fig. 3 V(x) as a function of protein concentration x(t)

2 穩態概率分布

在KARMAKAR和 BOSE[7]提出的ComK基因表達形成ComK蛋白的基因表達模型的基礎上，PAL等[9]在2013年提出了伴隨著隨機漲落的基本表達模型。根據非平衡系統計算物理和隨機過程的相關知識可知，系統的隨機漲落與概率函數密切相關[10]，從理論角度出發，主方程可為任何伴隨生化反應的系統的概率行為提供建模框架，但是對于一些復雜的網絡，特別是非線性系統，要求解主方程的概率函數是很難做到的。因此，基于這一問題，利用概率分布與生成函數之間的關系，重構系統的聯合概率函數，推導出以其對應的福克?普朗克方程[11]。在圖1模型中，隨機漲落在蛋白質的合成率a和 降解率b上均有影響，考慮這些因素，引入Gauss 色噪聲，即a=a+η(t)和b=b+ε(t)，則方程（3）重寫為朗之萬方程：

則方程（5）具有一般隨機動力學機制對應的Langevin方程：

ε(t) η(t)和 分別是Gauss乘性色噪聲和Gauss加性色噪聲，它們具有以下所示的統計性質：

式中，D和α 分 別為有色乘性以及加性噪聲強度，τ1、τ2是分別為其自關聯時間。

因為蛋白質濃度x是一個不可能小于0的值，所以有x≥0，這里p(x,t)用 來表示動態蛋白質濃度x在時刻t的概率分布函數，因此，朗之萬方程（5）對應的p(x,t)的近似福克?普朗克方程可以應用Novikov定理[12]和Fox近似方法[13]推導出，即概率分布隨時間的演化方程為[14]：

這里xs=x2由方程（3）得到，在定態情況下求解方程（8），可以得到其定態概率分布函數為

可化簡為

式中：

N為歸一化常數，

為了研究色噪聲強度和自關聯時間對概率分布函數影響，根據有效勢函數方程（10）所表達的穩態概率分布函數在基因轉錄調節系統中的解析式[15]，對方程（10）進行數值計算，在其他參數固定的情況下，做出了不同噪聲強度和不同自關聯時間的對蛋白質濃度影響的偽三維圖，如圖4（a） 和5（a）所示，其中，紅色區域內均表示該系統處于雙穩狀態。其他區域表示該系統處于單穩狀態，結果表明，系統的穩定性隨噪聲強度和自關聯時間的變化在穩定和不穩定之間變化隨著時間的變化，會出現噪聲削弱穩定性的現象。為方便討論，記基因處于高水平表達狀態為“on” 狀態，反之稱為“off ”狀態。

圖4 （a）：蛋白質濃度在噪聲強度和自關聯時間下的偽三維圖；（b）：概率分布函數P st(x) 作為自關聯時間t 1 的函數Fig. 4 (a): Pseudo - 3D diagram of protein concentration at the noise intensity and autocorrelation time; (b): The probability distribution function P st(x) is used as a function of the autocorrelation time t1

圖4（b）、5（b）給出了概率分布函數不同取值的噪聲強度和自關聯時間的變化曲線，從4（b）可以看出，隨著自關聯時間的增大，乘性噪聲D越大，蛋白質高濃度態逐漸變低，最后消失只剩低濃度態，這說明蛋白質處于失活狀態。從圖5（b）可以看出，隨著自關聯時間的增大，加性噪聲 α越大，蛋白質高濃度態逐漸變低，最后消失只剩低濃度態，即蛋白質處于“off ”的狀態。

圖5 （a）：蛋白質濃度在噪聲強度和自關聯時間下的偽三維圖；（b）：概率分布函數P st(x)作 為自關聯時間t 2 的函數Fig. 5 (a): Pseudo - 3D diagram of protein concentration at the noise intensity and auto correlation time; (b): The probability distribution function P st(x) is used as a function of the autocorrelation time t2

為了驗證方程(5)近似理論解的正確性，進行數值模擬是非常有必要的，應用歐拉算法模擬了朗之萬方程(5)和方程(7)，圖4（a）和5（a）給出了在不同噪聲強度和自關聯時間下對蛋白質濃度x(t)的時間序列對概率分布函數的影響變化，如圖4（b）和5（b）所示，對比近似的理論解和數值模擬得到模擬解，可以看出兩個方法的結果基本一致，這就意味著經過考慮噪聲強度和噪聲的自關聯時間來計算該模型對應的福克-普朗克方程得到的理論解是可信的。綜合上述，內外噪聲強度和其各自關聯時間增大，會導致一個穩態隨機地失去穩態性，切換到另一個穩定狀態，實現狀態的切換。因此，得出內外噪聲強度和它們的自關聯時間可以誘導基因的切換現象，換句話說，可以將噪聲強度和自關聯時間作為控制基因網絡開關切換的參數。

3 平均首通時間

雙穩系統中噪聲會影響兩穩態間的相互轉化，即系統從一個穩態在外力作用下出發穿越勢壘?V(x)進入另一個穩態的變化，如圖6所示，為得到一個兩態間轉化的確定時間值，通常用統計的方法[16]，直接考察平均首通時間(MFPT)，可把平均首通時間的精確表達式寫成：

圖6 外力作用下系統穩態的變化趨勢Fig. 6 The change of the steady state of the system under the action of external force

但方程（11）特別復雜，用一般方法很難處理，因此，在絕熱消磁近似的條件下，應用最快下降法[17?18]，在加性噪聲強度 α和乘性噪聲強度D遠小于勢壘?V(x) 時，利用MFPT的定義和最速下降法可以得到：

在基因表達系統中，躍遷是動力學特性中的瞬態性質，而這種瞬態性質常用平均首通時間來刻畫，圖7（a）給出了在乘性噪聲D的自關聯時間t1下 平均首通時間T12（off態x1躍 遷到on態x2）的函數圖像，在不同乘性噪聲強度下，MFPT會隨著t1的 增大呈現出單調遞減的趨勢，圖7（b）給出了在加性噪聲α 的自關聯時間t2下平均首通時間T12（（off態x1躍 遷到on態x2）的函數圖像，在不同加性噪聲強度下，MFPT會隨著t2的增大呈現出單調遞減的趨勢，這說明，隨著自關聯時間t1和t2的增大，噪聲強度越大，蛋白質濃度從低濃度態轉向高濃度態所需要的時間越少，即從“off ”狀態到“on”的轉換變得越容易，即加速了狀態之間的轉化，換句話說，自關聯時間會削弱高蛋白質濃度態的穩定性。

圖7 色噪聲自相關時間影響下，系統平均首達時（MFPT）的變化規律Fig. 7 The change of the mean first passage time (MFPT) of the system under the influence of color noise autocorrelation time

4 結 論

本研究應用諾維科夫理論和福克斯近似方法給出ComK基因表達系統中穩態概率分布函數的近似表達式，分析了色噪聲誘導蛋白質濃度轉換的現象以及該系統動力學特性中的瞬態性質，通過分析可以得到以下結論：噪聲強度及其自關聯時間的變化會影響概率分布函數曲線的峰值位置、峰值個數的變化，即噪聲強度及其自關聯時間確定會影響系統發生相變。在不同乘性噪聲和加性噪聲的噪聲強度下，隨著自關聯時間的增大，蛋白質從低濃度態躍遷至高濃度態所需時間減少，加速了狀態之間的轉化，說明噪聲強度和自關聯時間會引起蛋白質的濃度經歷了“開”→“關”的轉換，即自關聯時間在一定程度上會削弱蛋白質濃度態的穩定性，表明不同源噪聲的自關聯時間能誘導ComK基因在各種表型狀態之間切換，提高其生存概率。通過分析影響基因調節系統中蛋白質濃度轉化的因素和噪聲強度對系統的影響，發現蛋白質濃度轉化的一般規律和噪聲強度及自關聯時間的大小對系統產生的影響。本研究為基因方面的藥理學研究提供了一定的理論方法。