随机反应网络中瞬态概率的乘积形式近似【都灵大学】

可在www.sciencedirect.com在线获取理论计算机科学电子笔记277（2011）3-14www.elsevier.com/locate/entcs随机反应网络中瞬态概率的乘积形式近似AlessioAngius和Andr'asHorva'th都灵大学计算机科学系，意大利都灵angius@di.unito.it，horvath@di.unito.it摘要大多数描述随机反应网络的马尔可夫链都有一个巨大的状态空间。 这使得精确的分析是不可行的，因此除了模拟之外，唯一可行的方法就是近似。 在本文中，我们得到了这种马尔可夫链的瞬时概率的乘积形式近似。近似可以解释为一组相互作用的时间非齐次马尔可夫链与一个链的系统的每个反应物。 因此，计算复杂性仅随反应物的数量线性增长，并且可以对具有巨大状态空间的马尔可夫链进行近似。几个数值例子来说明这种方法。关键词：随机反应网络;马尔可夫链;瞬态分析;乘积型逼近。1介绍在[9，10]中，Gillespie提供了一般化学反应系统演化的随机描述。这种描述对应于一个简单的随机过程，即连续时间马尔可夫链（CTMC）。因此，原则上，我们有可能通过构造CTMC的无穷小生成矩阵并计算其指数来分析这样的系统。 计算 矩阵的指数一般不是直接的[17]，但是对于表示马尔可夫链的矩阵，已经开发了有效的和数值稳定的技术[12，19]。 然而，如果马尔可夫链的状态数非常大或无限大，即使这些技术也可能失败，这通常是考虑反应系统时的情况。已经提出了几种近似技术来克服巨大的状态空间的问题。基于CTMC遵循的轨迹关系和基于微分方程的系统描述[13，14]的平均场方法提供了系统行为的确定性近似平均场方法提供的近似值可视为模型的近似平均行为这个想法可以推广到更高的阶1571-0661 © 2011 Elsevier B. V.在CC BY-NC-ND许可下开放访问。doi：10.1016/j.entcs.2011.09.0314A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 277（2011）3和人口水平的联合时刻，导致更精确的近似[8]。其他技术通过直接在模型的状态空间上操作来获得近似由于状态空间可以是无限的，所以自然要约束所考虑的状态集[7]。此外，随着系统的发展，可能需要应用动态边界，以便在任何瞬态时间考虑状态空间的最可能区域。由于即使在缩减的状态空间上计算也会很慢，因此最近提出了更快的近似均匀化方法[16，23]。 处理巨大状态空间的另一种可能性在于聚合。一个自然的选择是聚合附近的状态[22，5]，而通过应用Escherow等价可以获得更复杂的聚合[4，6]。另一种分析表示反应网络的大型（或无限）马尔可夫链的方法是模拟。由于巨大的状态空间和大量的反应可以在很短的时间间隔内发生的事实，即使模拟也不是简单的。从[9]开始，一些论文提出了提高反应系统模拟效率这些方法中使用最多的是显式[11]和隐式[18] tau-leaping，它使用近似在单步中“跳过”许多反应，以及慢尺度随机模拟算法[ 3 ]，其目的是通过区分快速和缓慢反应来在本文中，我们提出了一种新的近似技术。我们保持模型的状态空间不变（即，我们不执行归约和聚合），但我们通过对模型的瞬时概率进行假设来简化分析特别地，我们假设瞬时概率是乘积形式的。这一假设使我们能够有一个紧凑的描述的瞬态概率，从而分析系统的状态空间，否则将是禁止的大小。与我们的技术最相关的结果是在[2]中提出的，其中作者提供了一个马尔可夫排队网络具有瞬时乘积形式的充分必要条件。条件是网络由无限服务器队列组成。显然，并非所有反应网络都对应于无限服务器队列网络，因此乘积形式假设一般不成立。然而，我们将通过数值实验表明，反应网络越接近无限服务器的网络，近似值就越好。本文的组织结构如下。在第2节中，简要描述了所考虑的反应系统在第3节中，我们介绍我们的近似方法。近似的性质在第4节中讨论。第5节提供了数值示例。第六节得出结论A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 277（2011）35ΣΣ.ΣM.Σdt=−p（X，t）λn=DTλncnm（四）！2模型配方我们考虑M个反应物（也称为物种），R1，R2，...，RM，通过N个反应相互作用：ManmRm−λ→nm=1MbnmRm，1 ≤n ≤ N。（一）m=1第n个反应消耗了a nm单位的反应物Rm，并产生了b nm单位的反应物Rm。a nm和bnm都是非负整数值，并将被组织成向量a n= |a n1，.，的nM|b n= |b n1，.，b nM|.我们用c nm= b nm−a nm表示反应物m上反应n的总效应，相应的矢量用c n = b nm − anm表示。|c n1，.，c nM|.第n个反应的速度由λn给出。有两种经典的方法将时间行为与（1）中的反应联系起来。随机方法将连续时间马尔可夫链（CTMC）与系统联系起来[9]。CTMC是离散状态，即，给定反应物在任意时刻t的量由整数给出，链的状态由一个整数向量。在状态X=|X1，...，X M|当Xm≥nm，1 ≤m≤M时，反应n是可能的。 我们将把关系式≥应用于向量，这意味着X≥an当且仅当Xm≥anm，1≤m≤M.如果反应n在状态X中是可能的，那么它的强度由下式给出：Mλnm=1Xmanm（二）并且相应的转换使CTMC从状态X到状态X+Cn。 我们将用p（X，t）表示CTMC在时间t处于状态X的概率这个量满足以下著名的查普曼-柯尔莫哥洛夫普通微分方程（ODE）dp（X，t）- 是的Xmn：X≥anm=1一nmn：X−<$cn≥anp（X−cn，t）λnm=1Xm−cnm的nm.（三）第二种确定性方法不是将反应描述为离散的过渡，而是假设反应以有限的量进行修改所涉及的物种。因此，反应物的量由连续值给出，并且常微分方程描述了系统的演化。通过应用质量作用动力学[20]，ODE是N MaRm（t）<$Ri（t）nin=1i=1一niM+6A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 277（2011）30、1个jλ3 0，j（j−1）λ2m=1DT11212DT21232λ32λ33λ30， 2λ2（j+ 1）λ3λ31， 0λ1λ3二、零2λ1（i−1）λ1λ3i，0iλ1Fig. 1. Lotka-Volterra模型的马尔可夫链的边界和一般状态（每个状态由猎物和捕食者的数量其中Rm（t）∈R≥0是反应物Rm在时间t的量。在设定了初始值Rm（0），1≤m≤M之后，（4）中的常微分方程组可以通过数值积分求解，并得到确定性的时间行为。上述确定性方法恰好是随机方法的CTMC此外，库尔茨还表明，[14]随着初始种群水平的增加和反应强度的相应调整（产生一系列所谓的水平依赖马尔可夫链），CTMC遵循的轨迹趋于确定性方法的ODE所在整个文件中，我们将使用著名的Lotka-Volterra模型来说明的方法。该模型由Lotka [15]和Volterra [21]独立提出，使用三个反应来描述两个种群在竞争中的进化这三种反应是猎物的生长：R1−λ→1 2R1，捕食者的增长：R1+R2−λ→22R2，捕食者死亡：R2−λ→3∅该模型的马尔可夫链如图1所示。模型的平均场近似由常微分方程提供R1（t）= λ R（t）−λ R（t）R（t），R2（t）= λ R（t）R（t）−λR（t）这导致沿着闭合曲线的振荡，除非系统在平衡状态下启动。3乘积形式近似为了得到所提出的近似，我们假设模型的瞬态概率是乘积形式的。这意味着，用p（x，m，t）表示在时间t有x个反应物m单位的概率，瞬态概率状态的性质可以写为p（X，t）=Mp（X m，m，t）。量p（x，m，t）0，0i−1，j i −1，j+1i，j−1（i−1）λ1jλ3i、jijλ2（j+ 1）λ3i，j+1（i+1）（j−1）λ2iλ1i+1，j−1i+1，jA. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 277（2011）37MΣJm=1DTDT=DT⎣−i=1λnp（Xi−cni，i，t）- 是的（六）λnp（x−cnm，m，t）Σn：x−cnm≥anm一nmi=1，i/=m满足ODEdp（x，m，t）=d<$X：Xm=xp（X，t）=<$dp（X，t）DT（五）通过应用（3）和产品形式假设，dp（x，m，t）Σ ⎡ΣλnX：Xm=xn：X≥an.Xi一ni- 是的Xi−cnin：X−cn≥ani=1（6）中求和的顺序可以颠倒，这导致：dp（x，m，t）dt=−λnn：x≥anmxanmp（x，m，t）Mi=1，i=/f（i，ani，t）+M- 是 的x−cnm其中，量f（i，j，t）与时间t时第i个反应物的量的第j阶阶乘矩密切相关，定义为：∞f（i，j，t）=K=j.k∈p（k，i，t）.（八）4近似的性质（7）中描述我们的产物形式近似的常微分方程可以解释为M个相互作用的时间非齐次马尔可夫链，其中对应于给定反应物的链通过考虑（8）中给出的其他链的量来模拟反应物参与的反应。图2描述了Lotka-Volterra模型的解释。由此可见，时间非齐次马尔可夫链的数值解技术，如[1]中提出的技术，可以应用于计算瞬态概率。让我们用rm，1≤m≤M表示概率p（x，m，t）不可忽略的x值的然后，描述近似值为mr m。这个数量随着反应的数量线性增长因此，该方法具有很好的扩展性。如第2节所述，如果我们考虑一系列具有递增初始状态的水平依赖马尔可夫链，则瞬态行为倾向于模型的平均场近似[13，14]。这同样适用于我们提出的近似。增M.MX：Xm=xp（Xi，i，t）+一nif（i，ani，t）（7）8A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 277（2011）3加初始人口水平会引起一系列的水平依赖，A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 277（2011）390121009998−→∅λ12λ1（i−1）λ1iλ1我f（2，1，t）λ22f（2，1，t）λ23f（2，1，t）λ2if（2，1，t）λ2（i+1）f（2，1，t）λ2f（1，1，t）λ22f（1，1，t）λ2（j−1）f（1，1，t）λ2jf（1，1，t）λ2λ32λ33λ3jλ3（j+ 1）λ3图二. Lotka-Volterra模型的近似时间非齐次马氏链;上半部分：食饵，下半部分：捕食者。100f（2， 2，t）99f（2， 2，t）98f（2， 2，t）2f（2， 2，t）f（2， 2，t）100⎝2⎠100⎝2⎠98⎝2⎠96⎝2⎠图三. R1+ 2 R21模型的时间非齐次马氏链的逼近;上半部分：R1，下半部分：R2.相互作用，时间非齐次马尔可夫链和瞬态行为，这一系列往往平均场近似。这意味着，在所提出的乘积形式近似和平均场近似之间的关系与原始马尔可夫链模型和平均场近似之间的关系相同。5数值算例作为第一个例子，我们考虑一个非常简单的模型，以显示近似无法提供准确的结果。 该模型由单反应R1+2R2−→1π，起始态为（100，100）。很明显模型满足不变量2（100−R1）= 100−R2，因此状态空间由51个状态组成，它们是（100，100），（99，98），（98，96），...，（50，0）。图3中描述了两个交互的时间非齐次马尔可夫链，它们代表了我们的乘积形式近似。在原始模型和它的近似模型之间有两个重要的区别近似不保持原始模型的不变量。 此外，在原模型中，反应物R1的量不能减少到50以下，而在近似模型中，它可以发生。因此，该近似适用于152个微分方程，并且需要比原始模型更多的计算然而，一般情况并非如此在模型包含更多的反应物，并没有非常严格的不变量的情况下，近似需要比原来的模型少得多的计算在图4中，我们描述了所涉及的量作为时间的函数。该近似精确地描述了两种反应物的平均值，并提供了很好的估计R2的方差，而它失败的R1的方差。 图5显示了数量的分布反应物的几个不同的瞬态时间。对于发行版，012J10f（1，1，t）2010A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 277（2011）3−→∅−→∅1009080706050403020100精 确 平 均 值 ，R1 近 似 值 平 均值 ， R1 精 确 平均 值 ， R2 近 似平 均 值 ， 平 均值，R20.0005 0.001 0.0015 0.0020.0025时间353025201510500.0005 0.001 0.0015 0.002 0.0025时间图四、R 1 + 2 R 2 − 1 → λ模型的M平均（左）和方差（右t）。0.40.350.30.250.20.150.10.0500 2040时间60 801000.40.350.30.250.20.150.10.0500 2040时间60 80 100图五. R1+ 2 R21模型：不同瞬态时间的R1（左）R2（右）分布（对于R2每个奇数值的概率为零;为了得到一个清晰的图形，我们只绘制非零概率）。0.00140.00120.0010.00080.00060.00040.00020产品形式近似，R1平均场近似值，R100.00050.0010.00150.0020.0025时间0.020.0180.0160.0140.0120.010.0080.0060.0040.0020产品形式近似，R2平均场近似，R200.00050.0010.00150.0020.0025时间见图6。R1+ 2 R21模型：建议的乘积形式近似和平均场近似提供的平均值的相对误差;左R1，右：R2。近似对于R2是好的，而对于R1是坏的。结果并不令人惊讶，因为该模型远非产品形式，并且对于R1，它引入了在原始模型中是不可能的。 在图6中，我们比较了通过乘积形式近似获得的平均值的精度和平均场方法的精度。可以看出，我们提出的方法，即使模型对它不利，也比平均场近似法提供了更精确的平均值作为第二个例子，我们考虑Lotka和Volterra的食饵-捕食者模型。 首先，我们在状态（2000， 2000）中启动模型，并使用反应速率 λ1 = 10，λ2 = 0。01，λ3 =10。通过模拟和乘积形式近似获得的捕食者数量的均值和方差如图所示精确方差，R1近似值方差，R1精确方差，R2近似方差，方差，R20.0001，准确0.0001，约0.0005，准确0.0005，约0.002，准确0.002，约0.0001，准确0.0001，约0.0005，准确0.0005，约0.002，准确0.002，约A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 277（2011）311300025002000150010005000在大约我是说，R2模拟。平均值，R2电话：+86-0511 - 8888888传真：+86-0511-8888888时间800000700000600000500000400000300000200000100000000.511.52时间2.533.54图第七章 Lotka-Volterra模型，λ1= 10，λ2= 0。01，λ3 = 10，初始状态（2000， 2000）：通过模拟和乘积形式近似的捕食者数量的平均300250200150100500在大约我是说，R2模拟。平均值，R2电话：+86-0511 - 8888888传真：+86-0511-8888888时间80000700006000050000400003000020000100000在大约方差，R2sim。方差，R20 0.5 1 1.522.5 3 3.5 4时间见图8。Lotka-Volterra模型，λ1= 10，λ2= 0。1，λ3 = 10和初始状态（200， 200）：通过模拟和乘积形式近似的捕食者数量的平均值（左）和方差（右）。图7.即使人口水平很高，平均值很快就会偏离稳定的振荡模式。乘积形式近似预测均值的稳定振荡，并提供与均值场方法非常相似的均值（图7中未示出，因为无法将其与乘积形式方法提供的均值区分开）。乘积形式近似提供的方差模式给出了原始模型中发生的事情的更精确的描述。在图8中，我们描述了与图7中相同的量，但从状态（200，200）开始，速率λ1= 10，λ2 = 0。1，λ3 = 10。具有这些参数的模型的平均场近似值与之前使用的参数相同。但随着食饵和捕食者数量的减少，模型偏离稳定振荡的速度加快。乘积形式近似无法捕捉这种行为，并且提供了均值和方差的不精确估计。图9显示了通过模拟和提出的近似值获得的捕食者灭绝概率随时间的变化。近似不能给出精确估计的原因是系统的行为强烈依赖于人口水平的相关性，而这并不是由乘积形式概率所捕获的我们的第三个例子是Lotka-Volterra模型的扩展版本。我们考虑两种类型的猎物和两种类型的捕食者。反应集是在大约我是说，R2模拟。平均值，R212A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 277（2011）3150.450.40.350.30.250.20.150.10.0500 0.5 1 1.5 2时间2.5 3 3.5 4见图9。Lotka-Volterra模型，λ1= 10，λ2= 0。1，λ3 = 10和初始状态（200， 200）：通过模拟和乘积形式近似的捕食者灭绝。6005004003002001000在 大 约 我 是说 ， R1 模 拟。平 均值 ， R1 近似 。 我是 说 ，R2 模 拟。 平均值，R20 0.5 1 1.52时间350003000025000200001500010000500000 0.5 1时间1.5 2见图10。 扩展的Lotka-Volterra模型：两种类型数量的均值（左）和方差（右）通过模拟和乘积形式近似的方法以下内容：猎物的生长：R1−1→02R1，R2−1→02R2，捕食者的生长：R1+R30−。0→152R3，R1+R4−0。→032R4，R2+R3−0。→022R3，R2+R40−。0→252R4，捕食者的死亡：R3−1→0，R4−→。平均场方法将稳定振荡与模型联系起来.由于该模型比原始的Lotka-Volterra模型有更多的物种和更多的反应，因此在更长的在图10和11中，我们显示了系统中所有物种的均值和方差。未显示平均场方法提供的平均值，因为它与乘积形式近似提供的平均值无法区分对于这个模型的产品形式的方法给出了很好的估计的平均值，并说明了物种的方差的行为。请注意，由于涉及四个物种，即使不考虑概率可忽略的状态，原始马尔可夫链也是巨大的。产品形式的方法需要解决一个有大约6000个方程的常微分方程系统，在普通的笔记本电脑上使用odeToJava包1大约需要2分钟。作为最后一个例子，我们考虑的模型是接近网络的信息，1可在http://www.netlib.org/ode/上查阅，由M.Patterson和R.J. 斯皮特里在大约灭绝，R2模拟。R2消光在大 约方差 ，R1sim。方差，R1近似 值 方 差 ， R2sim。方差，R2A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 277（2011）313400350300250200150100500在 大 约 我 是说 ， R3 模 拟。平 均值 ， R3 近似 。我 是说 R4模 拟 器 平 均值，R40 0.5 1 1.52时间350003000025000200001500010000500000 0.5 1时间1.5 2见图11。 扩展的Lotka-Volterra模型：两种类型数量的均值（左）和方差（右）捕食者的模拟和产品形式近似。14 7121086在大约平均值，R2，λ =0.14 sim。平均值，R2，λ =0.1在大约平均值，R2，λ =0.52sim。平均值，R2，λ =0.5在大约平均值，R2，λ =10sim。平均值，R2，λ =10 2 4 6 810时间65432100 2 4 6 8 10时间见图12。 无限队列的扰动网络：R2（左）和R4（右）的平均值。16141210864200 2 4 6 810时间765432100 2 4 6 8 10时间图十三. 无限队列的扰动网络：R2（左）和R4（右）的方差。Nite服务器队列以及因此它们的瞬时概率以乘积形式很好地近似。我们认为反应200. 5 2 1 1 1λ−→R1，R1−→R2，只有最后的反应才能使系统不再是无限服务器的网络。与该反应相关的速率越大，产物形式近似越差初始状态是|0，0，0，0|.在图12和13中，我们显示了R2和R4的均值和方差。两种反应物的平均值近似良好，而R4的方差捕获良好，R2的方差被低估。计算时间不到两秒。在大 约方差 ，R3sim。方差，R3近似 值 方 差 ， R4sim。方差，R4在大约平均值，R4，λ =0.1在大约平均值，R4，λ =0.5SIM.平均值，R4，λ=0.1SIM.平均值，R4，λ=0.5在大约平均值，R4，λ =1SIM.平均值，R4，λ =1在大约var.，R2，λ=0.1在大约var.，R2，λ=0.5SIM. var.，R2，λ=0.1SIM. var.，R2，λ=0.5在大约var.，R2，λ =1SIM. var.，R2，λ =1在大约var.，R4，λ=0.1在大约var.，R4，λ=0.5SIM. var.，R4，λ=0.1SIM. var.，R4，λ =0.5在大约var.，R4，λ =1SIM. var.，R4，λ =114A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 277（2011）36结论本文导出了表示反应网络的马氏链的瞬时概率的乘积形式近似。该方法的计算量仅与反应物的数量呈线性增长，因此它可以应用于相应的马尔可夫链的精确分析是不可行的反应网络。我们在几个例子上测试了该方法，发现如果平均场方法能很好地捕捉系统的平均行为，或者瞬态系统行为接近产物形式，那么所提出的近似提供了反应物数量的方差和分布的良好图像。引用[1] Arns，M.，P. Buchholz和A. Panchenko，关于非齐次连续时间马尔可夫链的数值分析，Informs Journalon Computing22（2010），pp. 416-432[2] 布舍里河J.和P. Taylor，暂态产品形式分布在网络中，离散事件动态系统：理论和应用3（1993），pp.375-396.[3] 曹玉，D. T. Gillespie和L. R. Petzold，慢尺度随机模拟算法，J Chem Phys122（2005）。[4] Chandy，K. M.，联合Herzog和L. S.吴，网络的参数分析，IBM期刊，R。D. 19（1975），pp.36比42[5] Ciocchetta，F.，A. Degasstrom，J. Hillston和M.关于从bio-pepa导出的ctmc和ode模型的一些研究。电子笔记理论。Comput. Sci. 229（2009），pp. 145-163[6] Cordero，F.，A. Hor v'ath，D. 马尼尼湖Napione，M. D. Pierro，S. Pav an，A. Picco，A. Veglio，M. 塞里诺F. Bussolino和G. Balbo，Simplification of a complex signal transduction model by the application ofinvariants and conventional equivalent server，Submitted to Theoretical Computer Science（2011）.[7] Dayar，T.，L. Mikeev和V. Wolf，关于随机Lotka-Volterra模型的数值分析，在：Proc. 数 值 算 法 的计算机方面的 研 讨 会 （CANA10），2010年，pp。第289-296页。[8] Engblom，S.，计算主方程的高维解的矩，应用数学。Comput180（2006），pp.498-515[9] Gillespie，D. T.，耦合化学反应的精确随机模拟，J. Phys. Chem.81（1977），pp. 2340-2361。[10] Gillespie，D. T.，化学主方程的严格推导，Physica A188（1992），pp. 404[11] Gillespie，D. T.，化学反应系统的近似加速随机模拟，J Chem Phys115（2001），pp. 1716-1733年。[12] Jensen，A.，Marko链作为Marko过程研究的辅助工具，Skandinavisk Aktuarietidskrift36（1953），pp. 87比91[13] Kurtz，T. G.，Solutions of ordinary diffusion equations as limits of pure jump Markov processes，Journal of Applied Probability1（1970），pp.49比58[14] Kurtz，T. G.，化学反应，J Chem Phys57（1972），pp.2976-2978年。[15] Lotka，A. J.，“Elements of Mathematical Biology,” Williams and Wilkins Company,[16] Mateescu，M.，V. Wolf，F. Didier和T. A. Henzinger，化学主方程的快速自适应均匀化，IET系统生物学4（2010），pp. 441-452.A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 277（2011）315[17] 莫勒角和C. V. Loan，Nineteen doubtful ways to calculate the exponential of a matrix，25 years later，SIAM Review45（2003），pp.pp. 3-49[18] Rathinam，M.，L. R. Petzold，Y. Cao和D. T. Gillespie，Sti-100-ness in stochastic chemicalreacting systems：The implicit tau-leaping method，J Chem Phys119（2003），pp.12784-12794。[19] 斯图尔特，W。J.， 普林斯顿大学出版社，1995年。[20] Voit，E. O.，“Computational Analysis of Biochemical Systems,” Cambridge University Press,[21] Volterra，V.，波动在丰富的一个物种认为数学，自然（1924年）。[22] 张杰，L. T. Watson和Y.曹，化学主方程的自适应聚集方法，国际。J. 计算生物学和药物设计2（2009），pp。134-148。[23] 张杰，L. T. Watson和Y. Cao，一种改进的化学主方程，计算机数学与应用59（2010），pp.573