加速RRA的隨機(jī)模擬算法

周　文，瞿　佳，柴　甜

(1.安徽師范大學(xué) 數(shù)學(xué)計(jì)算機(jī)科學(xué)學(xué)院，安徽 蕪湖　241000；2.無(wú)錫市港下實(shí)驗(yàn)小學(xué)，江蘇 無(wú)錫，214000)

?

加速RRA的隨機(jī)模擬算法

周文1，瞿佳1，柴甜2

(1.安徽師范大學(xué) 數(shù)學(xué)計(jì)算機(jī)科學(xué)學(xué)院，安徽 蕪湖241000；2.無(wú)錫市港下實(shí)驗(yàn)小學(xué)，江蘇 無(wú)錫，214000)

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

關(guān)鍵詞：隨機(jī)模擬算法；最后所有可能步進(jìn)法；RRA算法

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

精確的SSA算法每次只能模擬一個(gè)反應(yīng)，模擬速度較慢.為了提高模擬的速度，Gillespie等人提出加速的隨機(jī)模擬算法(τ-leaping)[3]，該算法通過(guò)計(jì)算生成一個(gè)跳躍時(shí)間τ，并假設(shè)在[t,t+τ)內(nèi)系統(tǒng)狀態(tài)量的改變非常?。蝗缓罄貌此呻S機(jī)數(shù)確定每一個(gè)反應(yīng)通道在時(shí)間區(qū)域內(nèi)的反應(yīng)次數(shù)k，計(jì)算這一時(shí)間區(qū)域內(nèi)系統(tǒng)中反應(yīng)時(shí)間的狀態(tài)改變量.除了τ-leaping算法，還有許多的算法從跳躍的步長(zhǎng)[4-6]方面考慮，也有從跳躍期間反應(yīng)發(fā)生次數(shù)[7-8]方面去考慮加速精確隨機(jī)模擬算法.

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

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

1背景知識(shí)

設(shè)一個(gè)生化反應(yīng)系統(tǒng)中存在M個(gè)反應(yīng)通道R1,R2，…,RM及N種化學(xué)物質(zhì)S1,S2，…，SN,其t時(shí)刻的動(dòng)力學(xué)狀態(tài)可由Xn(t),aj(t),vj,n=1…N,j=1…M表示，其中Xn(t)表示物質(zhì)Sn在t時(shí)刻的分子數(shù)目，aj(t),vj分別表示在t時(shí)刻反應(yīng)通道Rj的傾向函數(shù)與狀態(tài)改變向量.對(duì)給定的反應(yīng)通道在下一個(gè)無(wú)窮小的時(shí)間區(qū)間(t+τ,t+τ+dt)內(nèi)發(fā)生一次反應(yīng)的概率為ajτ.

1.1SSA算法

SSA算法在模擬過(guò)程中假設(shè)所有的反應(yīng)發(fā)生都是瞬時(shí)發(fā)生的，不考慮反應(yīng)過(guò)程時(shí)間的消耗，即反應(yīng)一旦發(fā)生就立即結(jié)束.我們?cè)谀M生化反應(yīng)系統(tǒng)中主要要研究以下兩個(gè)問(wèn)題：1)下一個(gè)反應(yīng)何時(shí)發(fā)生?2)究竟哪個(gè)反應(yīng)會(huì)發(fā)生?研究出這兩個(gè)問(wèn)題就可以得到系統(tǒng)狀態(tài)隨時(shí)間變化的軌道.由當(dāng)前時(shí)刻系統(tǒng)的動(dòng)力學(xué)狀態(tài)確定下一個(gè)反應(yīng)Rμ發(fā)生在時(shí)間區(qū)間(t+τ,t+τ+dτ)內(nèi)的概率是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)產(chǎn)生兩個(gè)[0,1)上均勻分布的隨機(jī)數(shù)r1和r2，根據(jù)表達(dá)式(2)-(3)，計(jì)算出τ和μ；

(4)更新系統(tǒng)的狀態(tài)x=x+vj及反應(yīng)時(shí)間t=t+τ；

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

1.2FAPS算法

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

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

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

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

所有需要的運(yùn)行次數(shù)結(jié)束后，對(duì)概率進(jìn)行正則化.

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

1.3RRA算法

RRA[9]算法包括兩個(gè)算法即RRA-τ算法和RRA-N算法.RRA算法在模擬生化反應(yīng)系統(tǒng)時(shí)都是選取2A→B為最合適的代表反應(yīng).

RRA算法的傾向函數(shù)a0，

(4)

由(4)式可以求解出

(5)

(5)式中x0表示生化反應(yīng)系統(tǒng)中物質(zhì)的總數(shù)目，故x0取正值舍負(fù)值.

結(jié)合上述參數(shù)，對(duì)RRA算法歸納如下：

算法2RRA算法

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

(4)計(jì)算生成反應(yīng)系統(tǒng)中跳躍步長(zhǎng)τ：

(5)生成在[t,t+τ)內(nèi)反應(yīng)的次數(shù)kj～B(ajτ,iseed),iseed為一個(gè)隨機(jī)數(shù)；

(6)更新系統(tǒng)的狀態(tài)x=x+kjv及反應(yīng)時(shí)間t=t+τ；

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

2改進(jìn)的RRA算法

2.1FAPS RRA算法

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

FAPS算法特有一個(gè)性質(zhì)是可以和其它隨機(jī)模擬算法相結(jié)合.考慮到這一點(diǎn)，我們將RRA算法與FAPS算法有效地結(jié)合，提出FAPS RRA算法.當(dāng)RRA算法運(yùn)行到需要到達(dá)時(shí)刻的最后時(shí)刻時(shí)，F(xiàn)APS算法去運(yùn)行且考慮反應(yīng)系統(tǒng)下一個(gè)狀態(tài)的所有可能反應(yīng)，更新概率為：

(6)

所有需要的運(yùn)行次數(shù)結(jié)束后，對(duì)概率進(jìn)行正則化.該算法總結(jié)如下：

算法3FAPS RRA算法

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

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

(5)生成在[t,t+τ)內(nèi)反應(yīng)的次數(shù)kj～B(ajτ,iseed),iseed為一個(gè)隨機(jī)數(shù)；

(6)更新系統(tǒng)的狀態(tài)x=x+kjv及反應(yīng)時(shí)間t=t+τ；

(7)根據(jù)方程(6)，更新每個(gè)反應(yīng)Rj,j=1…M發(fā)生的概率；

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

3數(shù)值模擬

為驗(yàn)證所提出的FAPS RRA算法的模擬性能，本節(jié)通過(guò)模擬兩個(gè)具體的數(shù)值實(shí)例，分別從模擬的精度與速度上與RRA算法進(jìn)行比較.這里的模擬工作均是用MATLAB實(shí)現(xiàn)的，其運(yùn)行環(huán)境為Windows 7系統(tǒng)，2Gb內(nèi)存，2.53G Hz處理器.

3.1降解聚合系統(tǒng)

該模型已在早期的文章(見[3，7])中使用過(guò)來(lái)檢驗(yàn)他們的算法.降解聚合系統(tǒng)包括以下4個(gè)反應(yīng)通道：

圖1　降解聚合系統(tǒng)反應(yīng)中物質(zhì)S2隨時(shí)間演化的濃度圖Fig.1　Concentration change of the molecular S2 for the decaying-dimerizing model

R1φ

R2

R3

R4

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

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

圖2　降解聚合系統(tǒng)反應(yīng)中x2(15)的直方圖距離　　　　　　圖3　降解聚合系統(tǒng)反應(yīng)中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算法分別運(yùn)行10000次，得到S2(t)和S3(t)在t=15時(shí)刻的直方圖距離對(duì)應(yīng)運(yùn)行次數(shù)的函數(shù)關(guān)系圖.從圖2，3上可以看出直方圖距離隨著運(yùn)行次數(shù)的增加而降低，直至某一誤差值.可以看出，對(duì)一個(gè)特定精度，F(xiàn)APSRRA算法比RRA算法對(duì)物質(zhì)S2大約需要近4倍的模擬次數(shù)，對(duì)物質(zhì)S3大約需要2倍的運(yùn)行次數(shù).例如，對(duì)模型模擬10000次后，對(duì)物質(zhì)S2(S3)使用FAPSRRA算法降至0.12(0.18)，而使用RRA算法降至0.41(0.31).

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

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

3.2線性聚合反應(yīng)系統(tǒng)

本實(shí)驗(yàn)?zāi)M的化學(xué)反應(yīng)系統(tǒng)為一個(gè)簡(jiǎn)單的線性聚合反應(yīng)系統(tǒng)(見文獻(xiàn)[6])：

R1.

(1)

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

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

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

圖5　線性聚合反應(yīng)中x2(15)的直方圖距離　　　　　　　　　　　　　圖6　線性聚合反應(yīng)中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算法得到線性聚合反應(yīng)中物質(zhì)S2和S3的分子數(shù)目x2(1)和x3(1)的直方圖距離.直方圖距離隨著時(shí)間增加而降低，直至降低到某一誤差值.由于對(duì)于任意給定的精度，RRA算法需要比FAPS RRA算法多將近4倍的運(yùn)行次數(shù).例如，為得到誤差0.4，使用RRA算法對(duì)物質(zhì)S2(S3)平均需要511(808)次，而使用FAPS RRA算法對(duì)物質(zhì)S2(S3)平均需要182(488)次.分析模擬運(yùn)行時(shí)間得出，為得到相同的直方圖距離FAPS RRA算法需要的CPU時(shí)間遠(yuǎn)小于RRA算法需要的CPU時(shí)間.

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

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

4總結(jié)

參考文獻(xiàn)：

[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]劉煥，彭新俊，周文，王翼飛.改進(jìn)的τ-leoping算法在生化反應(yīng)系統(tǒng)隨機(jī)模擬中的應(yīng)用[J].應(yīng)用科學(xué)學(xué)報(bào)，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

中圖分類號(hào)：O644

文獻(xiàn)標(biāo)志碼：A

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

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

基金項(xiàng)目：國(guó)家自然科學(xué)基金項(xiàng)目(11302002)；國(guó)家自然科學(xué)基金資助項(xiàng)目(61202156)；安徽高校省級(jí)優(yōu)秀青年人才基金重點(diǎn)項(xiàng)目(2010SQRL0256ZD).

收稿日期：2015-06-20

DOI：10.14182/J.cnki.1001-2443.2016.02.002

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