轨迹抽样模拟加速CTMC模拟方法

可在www.sciencedirect.com在线获取理论计算机科学电子笔记296（2013）183-197www.elsevier.com/locate/entcs较少随机样本的马尔可夫链模拟Dimitrios Milios1 斯蒂芬·吉尔摩2英国爱丁堡大学信息学院摘要我们提出了一种加速CTMC模拟方法，这是确切的意义上说，它产生的所有涉及的过渡。我们称我们的方法为轨迹抽样模拟，因为它从状态序列的分布和给定某个特定序列的时间分布中抽样。从轨迹空间而不是过渡空间采样意味着我们需要生成更少的随机数，该操作通常在计算上是昂贵的。从时间分布中采样涉及到用几何分布来近似支配逗留时间的指数分布。适当地选取近似参数可以保证模拟的随机过程与原始马尔可夫链的模拟几乎一致。我们的方法不依赖于系统的属性，它可以作为一种替代更有效的方法时，这些是不适用因保留字：马尔可夫链，仿真，轨迹，随机变量1介绍连续时间马尔可夫链（CTMCs）已经被用于描述表现出随机行为的系统多年。随机模拟是探索大规模CTMC瞬态和稳态特性的传统方法然而，存在诸如生物化学反应网络的模型，其状态空间即使对于这种方法也太大。标准CTMC模拟方法被称为直接法（DM）[9]。在非常大的模型的情况下，它具有很高的计算成本，因为它模拟了每一个可能发生的转变。已经提出了几种加速方法，它们要么是精确的，要么是近似的。精确方法通常涉及对标准算法的优化，例如下一个反应方法[8]，优化的DM[5]，对数1电子邮件：D. sms.ed.ac.uk2电子邮件：Stephen. ed.ac.uk1571-0661 © 2013 Elsevier B. V.在CC BY-NC-ND许可下开放访问。http://dx.doi.org/10.1016/j.entcs.2013.07.012184D. Milios，S.Gilmore/Electronic Notes in Theoretical Computer Science 296（2013）183[17][19][ 19][19 ][ 19]这些方法中的大多数涉及使用适当的数据结构，以便有效地生成仿真事件。例如，优化的DM利用依赖图来避免重新计算保持不变的速率。近似模拟方法试图跳过一些模拟事件，与精确方法相比，这导致过程明显更快。例如，τ-跳跃[10]将时间提前预先选择的τ，在此期间可能发生许多类似地，在R-跳跃[1]中，随机模拟通过预先定义的数量的转换触发来加速。K-leap[3]是一种通过特定数量的事件来推进时间的方法。所有这些方法都假设一次转换只会对系统的状态造成很小的改变。其他方法如[20]和[4]利用时间尺度分离的概念假设该模型可以分为两个子系统：慢和快。快速子系统的行为是近似的，而只有缓慢的子系统进行模拟。这种假设可能并不总是适用于任意模型，这意味着要么引入了显著的误差，要么近似方法未能加速模拟过程。我们提出的轨迹采样模拟（TSS）算法是DM的一种改进，可以被描述为几乎精确的，在这个意义上说，它可以是任意精确的。在DM的情况下，每个步骤需要从两个分布中采样：状态分布和时间分布，两者都以当前状态为条件。以类似的方式，TSS涉及从状态序列的分布中采样。这减少了生成的随机样本的数量，这意味着更快的模拟算法。该算法仍然是精确的，因为没有跳过过渡。同样的方法被扩展到给定某个特定序列的时间分布中进行采样。这是通过用离散随机变量近似指数分布的逗留时间来实现的。时间离散化允许我们考虑时间分布作为一个离散状态马尔可夫链，因此采用TSS技术。这种修改基本上使我们的方法近似，因为CTMC的随机行为的一部分被抑制。然而，我们已经证明，我们的方法在极限收敛到原来的过程的解决方案，一个事实，解释了术语我们表明，一个适当的选择的近似参数可以导致一个行为非常接近原来的CTMC，并在一个合理的加速在同一时间。我们的实现基于优化的DM（ODM），因此ODM被用作效率比较的基线我们的方法与[2]中的K-skip方法I有关，或简称为K-skip。虽然他们从状态序列分布中采样的策略是相似的，但他们从时间分布中采样的方法是不同的。为了减少决定逗留时间的随机样本，假设后续状态的退出率相似，它们近似于具有Gamma分布的k个指数随机变量的和这种假设对于许多生化系统是合理的，但是它可能会为某些模型引入错误，正如我们在实验部分所证明的那样D. Milios，S.Gilmore/Electronic Notes in Theoretical Computer Science 296（2013）183185Σ⎪⎪⎩ΣΣ−我们已经实现了K-skip在原来的文件中的描述，为了产生一些比较结果。我们用于K-skip的错误参数是0。01，这是原始作品中使用的最小值。我们还强调了[2]中没有考虑的一些计算问题，这些问题是由一个随机数用于产生整个轨迹的事实引起的。在第2节中，我们介绍了本文中使用的一些概念。第3节和第4节包含我们工作的理论细节第5节讨论了一些实施问题。实验结果见第6节。最后，我们总结了第7节的结论。2预赛CTMC可以表示为三元组（S，Q，π0），其中S是有限状态集，Q∈R|S| ×|S|是生成矩阵，π0是S上的初始概率分布。每个s∈S都与一个指数分布的随机变量Ls<$Exp（λs）相关联，其中λs=s′×=Qsss′是离开状态s的速率。CTMC的跳链是具有概率矩阵P的离散时间马尔可夫，其中：Qss′/Qs，s=/P′=100，ssJ和Qsi= 0sJ和Qs= 0其中Q=Q’（1）ss0，J，s ssS=s⎪Qs= 0s′×=s1，s =sJ且Qs= 0CTMC中的转换与取决于当前状态s的两个随机变量相关联：确定下一个状态的Xs和确定在s中花费的时间量的Ls。DM涉及从Xs和Ls采样以生成下一个事件。Xs的分布是以s为条件的，其概率质量函数由跳矩阵P的第s行给出。我们假设状态的排序为s sJ，如果s对应于转移矩阵中索引小于sJ的一行。如果sk−1是系统在k−1次跃迁后的状态，则从Xsk−1采样涉及使用均匀样本U<$U（0， 1）并以概率选择下一个状态skPr（X sk−1 = s k）= Pr（a sk Pr（s0：k−1）sk′1）来确定转变的任何Un个样本对于k= 1的特殊情况，序列概率将等于第一步的转移概率然后，我们可以计算下一步骤所需的均匀样本U k+1，并递归地更新as0：k+1和bs0：k+1，以使用等式（8）获得新的累积序列概率。不过，这种策略可能不是最佳的，因为它需要跟踪两个累积概率。一个更简洁、更有效的解决方案是将U k写成之前的均匀样本Uk−1。定理3.3如果Uk−1<$U（0， 1）用于绘制状态样本sk−1，则sk由下式确定：U=Uk−1−ask−1（十）kbsk−1 — ask−1证据 给定一个确定序列的均匀样本U，样本U k和U k−1可以写成（9）中指定的。如果我们解决了w.r.t. 在这两种情况下，我们得到以下等式：UkPr（s0：k−1） +as0：k−1=Uk−1Pr（s0：k−2） +as0：k−2其产生：U=Uk−1Pr（s0：k−2）−as0：k−1−as0：k−2（十一）Pr（s0：k−1）Pr（s0：k−1）使用（8），上面第二个分数的分子可以写为：as0：k−1−as0：k−2=P r（s0：k−2）ask−1+as0：k−2−Pr（s0：k−3）ask−2−as0：k−3我们还知道，s0：k−2 =Pr（s0：k−3）ask−2+as0：k−3 由于（8），所以我们可以将（11）重 写 为：U=Uk−1P r（s0：k−2）−ask−1Pr（s0：k−2）Pr（s0：k−1）Pr（s0：k−1）s0：k−1s0：k−1D. Milios，S.Gilmore/Electronic Notes in Theoretical Computer Science 296（2013）183189Σ根据（5）中序列概率的定义，我们有Pr（s0：k−1）=Pr（s0：k−2）Psk−2sk−1，这意味着：U=Uk−1−ask−1kPsk−2sk−1最后，我们可以将Psk−2sk−1写为累积概率的差，以获得等式（10）。Q从一些初始过渡开始，我们可以递归地生成均匀分布在0和1之间的随机样本的 整 个 序 列 如 果 前一 步 使 用 了 一 个 样 本 Uk−1<$U （ 0 ， 1 ） ， 则 我们知 道ak−11，模拟也仍然是精确的。然而，使用几何近似意味着我们有一个稍微改变的过程来近似原始的过程。为了评估这种近似的质量，我们构建了各种奖励的直方图（即，物种种群）在所使用的模型，因为它会被im-D. Milios，S.Gilmore/Electronic Notes in Theoretical Computer Science 296（2013）183195√实际上，比较这种大小的模型的整个状态空间分布。然后我们计算直方图距离[6]，这是真实分布和近似分布的直方图重要的是要注意，即使模拟是精确的，直方图距离也总是大于零，因为模拟产生的经验分布总是不同的。为了确定计算出的距离是否有意义，必须将其与相应的自我距离进行比较。直方图自距离取决于绘制的样本数小于自距离的直方图距离值意味着对于给定数量的样本，两个分布实际上是相同的。根据[6]，给定N个样本的平均直方图自距离的上限与分布无关，并且可以使用（4K）/（πN），其中K是数字 直方图中的间隔。 对于下面的例子，我们考虑了 K= 50。表2总结了所考虑模型中几个种属种群和时间点的直方图距离对于p= 1的TSS，一些距离略大于自距离（用斜体字表示的值）。这意味着我们有一个相当好的近似值，但对于所考虑的样本数量，可以观察到使用固定时间然而，当使用p = 0的TSS时，近似质量更好。1，正如所料。几乎在所有情况下，与真实分布的直方图距离与估计的自距离处于相同水平或小于自距离。这意味着观察到的误差在模拟过程固有引入的误差的限度内。 这些发现支持参数p = 0的TSS的主张。1用于几何近似是几乎精确的加速模拟方法。虽然K-skip被证明比我们的LacY和Goldbeter模型更有效，但表2表明在某些情况下它并不准确。LacY模型的大多数直方图距离都大于自距离或为两个版本的TSS计算的相应距离。似乎后续状态的速率相似的假设可能会引入一些错误，这使得K-skip不太适合某些模型。我们的方法一般化的系统，这个假设是无效的。此外，我们使用几何近似指定了每个事件发生的持续时间，这对某些系统来说可能很重要7结论轨迹抽样仿真需要较少的随机样本来生成马尔可夫链轨迹。这是通过使用单个随机数来确定整个转换序列来实现的我们已经证明，选择下一个跃迁所需的随机数可以用选择前一个跃迁的随机数来表示这将导致递归更新单个随机数196D. Milios，S.Gilmore/Electronic Notes in Theoretical Computer Science 296（2013）183表2106次模拟运行的直方图距离（自距离：0。0080）(a) LacY模型时间乳糖K-skipPLac产品乳糖TSS（p=PLac第一章产品TS乳糖S（p=PLac0。第一章产品2500.00700.00900.00640.00620.00050.00710.00450.00110.00545000.00400.00740.00870.00420.00110.00830.00300.00040.00717500.00410.00770.00740.00340.00040.00860.00500.00010.007610000.00440.00860.00810.00400.00040.00870.00320.00020.0079(b) 戈德贝特时间活性MK-skipActive XCTSS（p=1）活性MActive XCTSS（p = 0. 第一章active M active X C2.50.00650.00680.00740.00710.00390.00390.00480.00400.00365.00.00670.00590.00550.00670.00660.00540.00660.00530.00527.50.00700.00880.00370.00680.00820.00540.00800.00790.006010.00.00710.00630.00670.00550.00480.00810.00410.00610.0062它决定了整个状态序列在CTMC的情况下，使用第二随机数来确定该序列的长度。同样的概念已经被用于用具有控制这种近似的质量的参数p的几何分布来近似指数分布的时间。我们模拟了两个不同性质的生化模型，以评估我们的方法的效率和准确性。实验结果表明，我们的方法比ODM快约15-20%，而观察到的误差可以忽略不计。找到了一种类似的K-skip方法I更有效，但在某些情况下不太准确。因此，TSS可以被认为是K-跳跃的替代方案，在这种情况下，这可能是不合适的。对于要采样的轨迹的长度k，也有一些实际的考虑。k值太大可能导致丢失可能的状态序列，而值太小将使轨迹采样仿真退化到ODM。我们在实验中使用k= 10，但在较大模型的情况下，我们必须将k设置为较小的值。我们认为k= 5甚至适用于非常大的模型。例如，给定一个有500个反应的模型，我们有：500kΩ 3。125× 1013 253.确认作者得到了综合系统生物学中心（CISB）的支持，该中心由BBSRC和EPSRC资助，参考BB/D 019621/1。引用[1] A. 俄歇， P. 夏特兰，和 P. Koumoutsakos。R-leaping：通过反应跳跃加速随机模拟算法。化学物理学报，125（8）：084103，2006.[2] X. Cai和J. Wen.耦合化学反应随机模拟的高效精确和K-skip方法。化学物理学报，131（6）：064108，2009。[3] X. Cai和Z.徐加速耦合化学反应随机模拟的K-leap方法化学物理学报，126（7）：074102，Feb.2007年D. Milios，S.Gilmore/Electronic Notes in Theoretical Computer Science 296（2013）183197[4] Y. Cao，中国粘蝇D.T. Gillespie和L.R. 佩佐德慢尺度随机模拟算法。化学物理学报，122（1）：14116，2005.[5] Y. Cao，H. Li和L.佩佐德化学反应系统随机模拟算法的有效公式。化学物理学报，121（9）：4059[6] Y. Cao和L.佩佐德 化学反应系统随机模拟的精度限制和误差测量计算物理杂志，212（1）：62006年。[7] F.乔凯塔和J.希尔斯顿。Bio-PEPA：生物系统建模和分析框架。Theoretical Computer Science，410（33-34）：3065[8] M.吉布森和J.布鲁克。 多物种化学系统的高效精确随机模拟很多渠道。物理化学杂志，104（9）：1876[9] D. T.吉莱斯皮耦合化学反应的精确随机模拟。The Journal of Physical Chemistry，81（25）：2340[10] D. T.吉莱斯皮化学反应系统的近似加速随机模拟。化学物理学报，115（4）：1716[11] A.戈德比特有丝分裂振子中细胞周期蛋白和cdc2激酶的最小级联模型。美国国家科学院院刊，88（20）：9107[12] J. Himmelspach和A. M.乌马赫JAMES II建模与仿真框架。2009年高性能计算系统生物学国际研讨会，第101[13] A.詹森马尔可夫链作为马尔可夫过程研究的辅助工具。斯坎德Aktuarietidskrift，36：87 - 91 ，1953.[14] A. M.基尔泽克STOCKS：Gillespie算法对生化系统的随机动力学模拟。Bioinformatics，18（3）：470[15] P. 勒厄耶TestU01：一个用于随机数生成器经验测试的C库。ACM Transactions on MathematicalSoftware，33（4）：1[16] L 'Ecuyer和E.布斯特在Java中使用SSJ进行模拟。2005年冬季模拟会议论文集，第611-620页[17] H. Li和L.佩佐德化学反应系统离散随机模拟的对数直接法。加州大学圣巴巴拉分校计算机科学系技术报告，2006年。[18] M. Matsumoto和T.西村梅森扭转器：一个623维均匀伪随机数发生器。ACM Transactions on Modelingand Computer Simulation，8（1）：31998年[19] E. Mjolsness，D. Orendor R.J.，P. Chatelain，and P. Koumoutsakos.一种精确的加速随机模拟算法。化学物理杂志，130（14）：1441102009年[20] C. V.Rao和A. P·阿金随机化学动力学和准稳态假设：应用于吉莱斯皮算法。化学物理学报，118（11）：4999