加速RRA的隨機模擬算法

周　文，瞿　佳，柴　甜

(1.安徽師范大學 數學計算機科學學院，安徽 蕪湖　241000；2.無錫市港下實驗小學，江蘇 無錫，214000)

?

加速RRA的隨機模擬算法

周文1，瞿佳1，柴甜2

(1.安徽師范大學 數學計算機科學學院，安徽 蕪湖241000；2.無錫市港下實驗小學，江蘇 無錫，214000)

摘要：文章提出“最后所有可能步進的RRA算法(Final All Possible Step RRA, FAPS RRA)”.該算法改變了RRA算法的數據收集方式，目的是在不失算法精度的前提下，減少模型的運行次數，有效地提高算法的運行速度.仿真實驗結果表明：在相同精度的要求下，FAPS RRA算法比RRA算法的模擬運行速率有顯著提高.

關鍵詞：隨機模擬算法；最后所有可能步進法；RRA算法

在均勻混合的生化反應系統中，分子數目隨著時間的演化通常是通過解常微分方程組得到的.常微分方程組中每個方程表達了一種化學物質的分子數目隨時間的改變率，自變量為分子數目.在建立常微分方程時，假設化學反應系統隨著時間的演化是連續且確定的，分子數目是以離散的整數數目而變化，而且時間演化不是確定性的過程.顯然這樣的物理假設是不現實的，預測結果也不理想.Gillespie于1976年最早提出將齊性空間中的化學反應系統的演化看做是具有離散狀態的連續時間的馬爾科夫過程.該過程的動力學特征可由化學主方程(CME)表示，化學主方程所表達的系統狀態的統計性質可由精確的隨機模擬算法(Stochastic Simulation Algorithm,SSA)[1-2]得到.

精確的SSA算法每次只能模擬一個反應，模擬速度較慢.為了提高模擬的速度，Gillespie等人提出加速的隨機模擬算法(τ-leaping)[3]，該算法通過計算生成一個跳躍時間τ，并假設在[t,t+τ)內系統狀態量的改變非常小；然后利用泊松隨機數確定每一個反應通道在時間區域內的反應次數k，計算這一時間區域內系統中反應時間的狀態改變量.除了τ-leaping算法，還有許多的算法從跳躍的步長[4-6]方面考慮，也有從跳躍期間反應發生次數[7-8]方面去考慮加速精確隨機模擬算法.

許多改進的隨機模擬算法都是以加速算法本身為目的的，而未對生化反應系統做更多的研究.RRA[9]算法的提出開拓出新的研究領域，它主要思想是找出反應系統中能代表系統中所有反應的一個反應.文獻[9]中最終討論出具有代表性的反應是2A→B.我們可以理解為整個反應系統中所有的反應由2A→B代表，并求解出2A→B反應的跳躍時間，利用泊松隨機數確定代表反應在時間區域內的反應次數kj,計算這一時間區域內系統中反應時間的狀態改變量kjv，最后更新系統反應狀態.

上述的諸多算法在模擬一個生化反應系統時，都未對算法輸出數據進行處理.我們知道要得到高的精度，往往要對一個模型進行多次重復的模擬.如果模擬的模型輸出的數據量非常大時，就會影響算法的模擬效率.FAPS[10]算法可以有效地解決這一難題，該算法是在其它的隨機模擬算法運行到最后一步才去執行，并改變其它算法數據的收集方式，減少算法的運行次數.為了提高RRA算法的運行效率，本文提出改進的RRA隨機模算法-(FAPSRRA);并通過兩個具體的生化反應系統證明我們提出的算法的有效性.文章的具體安排如下：第一節介紹了生化反應系統中精確的隨機模擬算法(SSA)算法，RRA算法，以及最后所有可能步進法(FAPS)；第二節給出了本文提出的改進的RRA算法(FAPSRRA)；第三節通過兩個具體的生化反應系統證明我們提出的算法的有效性；最后對文章做總結.

1背景知識

設一個生化反應系統中存在M個反應通道R1,R2，…,RM及N種化學物質S1,S2，…，SN,其t時刻的動力學狀態可由Xn(t),aj(t),vj,n=1…N,j=1…M表示，其中Xn(t)表示物質Sn在t時刻的分子數目，aj(t),vj分別表示在t時刻反應通道Rj的傾向函數與狀態改變向量.對給定的反應通道在下一個無窮小的時間區間(t+τ,t+τ+dt)內發生一次反應的概率為ajτ.

1.1SSA算法

SSA算法在模擬過程中假設所有的反應發生都是瞬時發生的，不考慮反應過程時間的消耗，即反應一旦發生就立即結束.我們在模擬生化反應系統中主要要研究以下兩個問題：1)下一個反應何時發生?2)究竟哪個反應會發生?研究出這兩個問題就可以得到系統狀態隨時間變化的軌道.由當前時刻系統的動力學狀態確定下一個反應Rμ發生在時間區間(t+τ,t+τ+dτ)內的概率是p(τ,μ|x,t)dτ，可求解出

p(τ,μ|x,t)=aμ(x)exp(-a0(x)τ)，

(1)

令

(2)

(3)

算法1SSA算法

(1)初始化：x=x(t0),t0=0,tfinal=T;

(3)產生兩個[0,1)上均勻分布的隨機數r1和r2，根據表達式(2)-(3)，計算出τ和μ；

(4)更新系統的狀態x=x+vj及反應時間t=t+τ；

(5)運行算法直至t>T終止.

1.2FAPS算法

在微生物學的領域中，生化反應系統都是一些大型復雜的系統.而在實際應用模型中，若采用SSA算法，模型的模擬速度非常慢，影響模擬的效率.針對這一問題，Lipshtat提出了所有可能步進方法(APS)[11].它是基于精確的隨機模擬算法的基礎上對其輸出的數據進行分析.FAPS[10]算法不同于APS算法，它是在隨機模擬算法運行到最后時刻時考慮所有可能反應狀態而不去實際地運行它們.在到達最后時刻的前一步，我們假設在SSA中系統的狀態是X(t)=x,有M個可能的狀態分別會以概率aμ(x)/a0(x),μ=1，…，M，發生.SSA算法僅根據點概率aμ(x)/a0(x)隨機地選取下一個反應Rj，然后更新x+vj的概率為

q(x+vj)→q(x+vj)+1

在FAPS算法中，我們考慮虛擬的所有可能的反應，更新概率如下：

q(x+vj)→q(x+vj)+(aj(x)/a0(x))τj,j=1…M.

所有需要的運行次數結束后，對概率進行正則化.

FAPS算法是在隨機模擬算法運行到最后時刻考慮其所有可能的反應狀態，因此它不影響隨機模擬算法的本身，僅是改變了數據的收集方式.我們已經知道FAPS算法是在其它算法到最后時刻才運行的，因此它和其它隨機模擬算法能起互相補充的作用，這樣就可以使算法在不失精度的前提下提高模擬的效率.

1.3RRA算法

RRA[9]算法包括兩個算法即RRA-τ算法和RRA-N算法.RRA算法在模擬生化反應系統時都是選取2A→B為最合適的代表反應.

RRA算法的傾向函數a0，

(4)

由(4)式可以求解出

(5)

(5)式中x0表示生化反應系統中物質的總數目，故x0取正值舍負值.

結合上述參數，對RRA算法歸納如下：

算法2RRA算法

(1)初始化：x=x(t0),t0=0,tfinal=T;

(4)計算生成反應系統中跳躍步長τ：

(5)生成在[t,t+τ)內反應的次數kj～B(ajτ,iseed),iseed為一個隨機數；

(6)更新系統的狀態x=x+kjv及反應時間t=t+τ；

(7)運行算法直至t>T終止.

2改進的RRA算法

2.1FAPS RRA算法

RRA-τ算法與RRA-N算法的運行步驟基本相同，則在下文討論的RRA算法均以RRA-τ算法為例.

FAPS算法特有一個性質是可以和其它隨機模擬算法相結合.考慮到這一點，我們將RRA算法與FAPS算法有效地結合，提出FAPS RRA算法.當RRA算法運行到需要到達時刻的最后時刻時，FAPS算法去運行且考慮反應系統下一個狀態的所有可能反應，更新概率為：

(6)

所有需要的運行次數結束后，對概率進行正則化.該算法總結如下：

算法3FAPS RRA算法

(1)初始化：x=x(t0),t0=0,tfinal=T,p(x)=0,q(x)=0;

(2)重復運行步驟(3)-(5)，直至t+τ≥T；

(5)生成在[t,t+τ)內反應的次數kj～B(ajτ,iseed),iseed為一個隨機數；

(6)更新系統的狀態x=x+kjv及反應時間t=t+τ；

(7)根據方程(6)，更新每個反應Rj,j=1…M發生的概率；

(8)完成L次運行后，q(x)歸一化為p(x)=q(x)/∑xq(x).

3數值模擬

為驗證所提出的FAPS RRA算法的模擬性能，本節通過模擬兩個具體的數值實例，分別從模擬的精度與速度上與RRA算法進行比較.這里的模擬工作均是用MATLAB實現的，其運行環境為Windows 7系統，2Gb內存，2.53G Hz處理器.

3.1降解聚合系統

該模型已在早期的文章(見[3，7])中使用過來檢驗他們的算法.降解聚合系統包括以下4個反應通道：

圖1　降解聚合系統反應中物質S2隨時間演化的濃度圖Fig.1　Concentration change of the molecular S2 for the decaying-dimerizing model

R1φ

R2

R3

R4

在模擬中，取初始分子數目分別為x1(0)=x2(0)=1000,x3(0)=0,各反應的反應速率系數分別為c1=0.1,c2=0.002，c3=0.5,c4=0.04，反應參數ε=0.3.我們分別用SSA算法，RRA算法和FAPSRRA算法在t∈[0,15]內對該模型進行模擬，結果如下：

圖1為分別使用SSA算法和RRA算法得到x2(t)在[0，15]內的分子數目隨時間變化的函數圖.從圖1中可以看出，兩條曲線的趨勢基本相同.對反應模型進行1000次模擬，SSA算法需要平均CPU消耗時間比RRA算法需要的平均CPU消耗時間多大約1倍.然后分別使用RRA算法和FASP RRA算法對該模型進行模擬得到圖2，3.

圖2　降解聚合系統反應中x2(15)的直方圖距離　　　　　　圖3　降解聚合系統反應中x3(15)的直方圖距離Fig.2　Histogram distance of x2(15) for the　　　　　　　　　　　Fig.3　Histogram distance of x3(15) for thedecaying-dimerizing model　　　　　　　　　　　　　　　　　　　　　 decaying-dimerizing model

圖2，3表示使用RRA算法和FAPS RRA算法分別運行10000次，得到S2(t)和S3(t)在t=15時刻的直方圖距離對應運行次數的函數關系圖.從圖2，3上可以看出直方圖距離隨著運行次數的增加而降低，直至某一誤差值.可以看出，對一個特定精度，FAPSRRA算法比RRA算法對物質S2大約需要近4倍的模擬次數，對物質S3大約需要2倍的運行次數.例如，對模型模擬10000次后，對物質S2(S3)使用FAPSRRA算法降至0.12(0.18)，而使用RRA算法降至0.41(0.31).

為了更直觀地描述圖2，3的數據關系，給出表1統計出使用RRA算法和FAPS RRA算法得到相同直方圖距離的CPU時間比較.分析表中統計的數據得出，當直方圖距離相同時，使用FAPS RRA算法在不失精度前提下能有效地提高RRA算法的模擬速率，提高模擬效率大約2倍.

表1　不同算法得到相同的直方圖距離的CPU時間比較

3.2線性聚合反應系統

本實驗模擬的化學反應系統為一個簡單的線性聚合反應系統(見文獻[6])：

R1.

(1)

用這個模型我們可以看到由SSA算法，RRA算法和FASP RRA算法得到物質S2的分子數目x2隨時間變化的函數圖和x2(1),x3(1)的直方圖距離，然后比較相同的精度時他們的運行需要的CPU時間.

我們模擬這個系統在時間間隔[0，1]內演化，其中初始條件為x(0)=(5000,5000,0),反應率常數c=(0.1,0.01,1)，反應參數ε=0.03.分別使用SSA算法，RRA算法和FAPSRRA算法得到圖4及圖5，6.

圖4表示分別使用SSA算法和RRA算法得到線性聚合反應中物質S2的分子數目隨時間變化的函數圖.從圖中可以看出，兩條曲線的趨勢基本相同，反應中物質S2的分子數目在[0，1]內隨著時間變化而減少，在t=1時刻，x2(1)=2252,這與在線性聚合反應中，物質S2在反應通道R1中是生成物，在反應通道R2和R3中是反應物相吻合.

圖5　線性聚合反應中x2(15)的直方圖距離　　　　　　　　　　　　　圖6　線性聚合反應中x3(15)的直方圖距離Fig.5　Histogram distance of x2(15) for the linear　　　　　　　　　　　　　Fig.6　Histogram distance of x3(15) for the linearprototype kinetic system　　　　　　　　　　　　　　　　　　　　　　 prototype kinetic system

圖5，6表示分別使用FAPS RRA算法和RRA算法得到線性聚合反應中物質S2和S3的分子數目x2(1)和x3(1)的直方圖距離.直方圖距離隨著時間增加而降低，直至降低到某一誤差值.由于對于任意給定的精度，RRA算法需要比FAPS RRA算法多將近4倍的運行次數.例如，為得到誤差0.4，使用RRA算法對物質S2(S3)平均需要511(808)次，而使用FAPS RRA算法對物質S2(S3)平均需要182(488)次.分析模擬運行時間得出，為得到相同的直方圖距離FAPS RRA算法需要的CPU時間遠小于RRA算法需要的CPU時間.

為進一步表明FAPS RRA算法的優越性，列出表2.表2為分析圖5，6中的數據，統計出分別使用FAPS RRA算法和RRA算法得到相同直方圖距離時對應的運行時間.從上表中可以看出FAPS RRA算法相比RRA算法在幾乎不失精度的情況下，平均消耗的CPU時間比RRA算法平均消耗的CPU時間少4倍多.

表2　不同算法得到相同的直方圖距離的CPU時間比較

4總結

參考文獻：

[1]GILLESPIE D T. A general method for numerically simulating the stochastic time evolution of coupled chemical reaction[J]. J Phys Chem, 1996,22(4):403-434.

[2]GILLESPIE D T. Exact stochastics simulation of coupled chemical reactions[J]. J Phys Chem, 1997,81(25_:2340-2361.

[3]GILLESPIE D， PETZOLD L R. Improved leap-size selection for accelerated stochastic simulation[J]. J Phys Chem, 2003,119(16):8229-8234.

[4]CHATTERJEE A, VLACHOS D. An overview of spatial microscopic and accelerated kinetic monte carlo methods[J]. J Comput-Aid Mater Des, 2007,14(2):253-308.

[5]CAO Yang, GILLESPIE D T, PETZOLD L R. Efficient step size selection for the tau-leaping simulation method[J]. J Phys Chem, 2006,124(4):044109.

[6]劉煥，彭新俊，周文，王翼飛.改進的τ-leoping算法在生化反應系統隨機模擬中的應用[J].應用科學學報，2009，27(3)：266-271.

[7]AUGER A, CHATELAIN P, KOUMOUSAKOS P. R-leaping: accelerating the stochastic simulation algorithm by reaction leaps[J]. J Phys Chem, 2006,125(8):084103.

[8]CAO Yang, XU Zhaouyi. K-leap method for accelerating stochastic simulation of coupled chemical reactions[J]. J Phys Chem, 2007,126(7):1.

[9]SHANTANU Kadam, KUMER Vanka. A new approximate method for the stochastic simulation of chemical systems: the representative reaction approach[J]. J Comput Chem, 2012,(33):276-285.

[10]ZHOU Wen, PENG Xinjun, LIU Xiang, et al. “Final all possible steps” approach for accelerating stochastic simulation of coupled chemical reactions[J]. Appl Math, 2008,29(3):379-387.

[11]LIPSHTAT A. An “all possible steps” approach to the accelerated use of Gillespie's algorithm[J]. J Chem Phys, 2007,26(18):184103.

[12]DANIEL T, GILLESPIE. Approximate accelerated stochastic simulation of chemically reacting systems[J]. Chem Phys, 2001,115(4):16-18.

A Stochastic Simulation Algorithm for Accelerating RRA

ZHOU Wen1,QU Jia1,CHAI Tian2

(1.College of Mathematics and Computer Science， Anhui Normal University, Wuhu 241000, China；2.Gangxia Experimental Primary School, Wuxi 214000, China)

Abstract:A “Final All Possible Step RRA” algorithm is proposed in this paper. The algorithm changes the data collection of RRA algorithm. Our proposed algorithm aims to minimize the running time while preserving the accuracy. So it can improve the efficiency of the RRA algorithm. The result shows that our proposed algorithm is more effective than RRA algorithm.

Key words:stochastic simulation algorithm; final all possible step method; RRA algorithm

中圖分類號：O644

文獻標志碼：A

文章編號：1001-2443(2016)02-0109-06

作者簡介：周文(1980-)，安徽桐城人，副教授，碩導，研究方向：生物信息學.

基金項目：國家自然科學基金項目(11302002)；國家自然科學基金資助項目(61202156)；安徽高校省級優秀青年人才基金重點項目(2010SQRL0256ZD).

收稿日期：2015-06-20

DOI：10.14182/J.cnki.1001-2443.2016.02.002

引用格式：周文，瞿佳，柴甜.加速RRA的隨機模擬算法[J].安徽師范大學學報：自然科學版，2016，39(2)：109-114.