加速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.