可在www.sciencedirect.com在线获取理论计算机科学电子笔记275（2011）23-36www.elsevier.com/locate/entcs生物模型的Monte Carlo EM方法AlessioAngius1Andr'asHorva'th2都灵大学计算机科学系，意大利都灵摘要在模拟生物现象时，通常的情况是，所涉及的相互作用的结构和效果是已知的，但相互作用的速率既不知道，也不容易通过实验确定。本文在一般和抽象的背景下讨论了反应网络速率参数的估计问题。特别是，我们考虑的情况下，所研究的现象是随机的，连续时间马尔可夫链（CTMC）是适合其建模。此外，我们假设所研究的系统的演化不能被连续地观察到，而只能在离散的采样点之间发生大量的反应随机反应网络的参数估计通常采用极大似然原理。 在本文中，我们描述了如何期望最大化（EM）方法，是一种在不完全数据情况下的极大似然估计技术，可用于估计反应网络的动力学速率。 特别地，由于底层CTMC的巨大状态空间，方便的是使用EM方法的这种变体，即MonteCarlo EM（MCEM）方法，其利用模拟来分析模型。 我们表明，在质量作用动力学的情况下，MCEM方法的应用导致高效且令人惊讶的简单估计过程。我们提供的例子来说明该方法的特点，并表明它是适用的情况下，涉及多个物种的反应系统关键词：参数估计，质量作用动力学，最大似然，期望最大化方法1背景正如Gillespie在[5]中所描述的，生化系统的时间行为可以用随机过程来描述，特别是用连续时间马尔可夫链（CTMC）来描述。为了对描述所研究现象的CTMC模型进行完整描述并能够进行分析，动力学速率的估计是在这种情况下，参数估计本质上是一个优化问题，其目的是找到一组参数，使得模型能够以高概率再现实验观察结果。这个问题并非微不足道，原因有几的1电子邮件地址：angius@di.unito.it2电子邮件地址：horvath@di.unito.it1571-0661 © 2011 Elsevier B. V.在CC BY-NC-ND许可下开放访问。doi：10.1016/j.entcs.2011.09.00324A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）23所研究的现象可能非常复杂，其中几种试剂通过许多反应相互作用此外，认为过程是完美和连续可观察的往往是不现实的。特别地，测量技术通常不能将系统行为观察为连续过程，并且仅在有限的一组时刻提供对系统状态的观察。此外，连续的时刻可以彼此相距很远，以至于在它们之间发生相当多的反应。这意味着我们必须面对一个不完整数据的优化问题。在不完全数据情况下，A. Dempster在[3]中。EM方法的基本思想是重建期望中的缺失数据，并应用优化来寻找使重建的完整数据的概率最大化的参数。通常情况下，丢失数据的精确重建是一项艰巨的任务。在这些情况下，如Wei和Tanner在[15]中提出的，可以使用模拟来完成数据，这种方法称为Monte Carlo EM（MCEM）方法。本文将MCEM方法应用于随机反应网络中动力学速率的估计。特别是，我们认为随机反应网络根据质量作用动力学的发展，并表明，MCEM方法导致一个简单而有效的估计过程。有几项工作是通过应用优化方法来估计反应网络的动力学速率[10，4，12]。然而，大多数这些作品不考虑随机性，但适用于一个确定性的观点，正在研究的现象的演变。理论上，可以将确定性模型获得的速率转换为可用于随机设置的速率，但正如[11]中所这一观察导致试图根据[5]引入的随机观点给出动力学速率的估计。在[7，1]中使用贝叶斯推断方法，在[11，2，13]中应用最大似然方法。相对于上面提到的方法，我们的方法的优势在于它只适用于有限的信息（即，有可能将其应用于非常不频繁的观察，在这些观察之间发生数千个反应），并且不涉及繁重的优化任务。本文的组织结构如下。在下一节中，我们提供参考随机模型。然后简要介绍了EM方法的一般情况。随后，我们描述了应用的MCEM方法的随机反应网络的参数估计。最后一节是专门的数值说明所提出的方法。在最后一节中，我们得出了结论。2考虑模型本文考虑一个连续时间马尔可夫链（CTMC），它描述了M试剂通过R反应的相互作用.系统的状态由提供试剂量的M个整数的向量给出。反应i，1≤i≤R的反应由M个整数的向量ei描述，如果反应i，A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）2325.ΣK1i出现在状态x中，则下一个状态是xJ=x+ei。CTMC的跃迁代表反应，其强度取决于系统的状态。我们用fi（x）表示提供反应i在状态x中的强度的函数。在整个论文中，我们假设强度的形式为fi（x）=kiMj=1xjai，j（一）其中常数ai ，j提供了第i个反应的化学计量，其中ai ，j∈N并且ki是动力学速率常数（即，它描述交互的速度）。因此，该模型对应于质量作用动力学，CTMC正是由Gillespie [5]的经典算法模拟的作为模型的一个例子，考虑以下反应系统[6]，其也将用于数值示例中DNAof−→ DNA上，DNAon−k→2去掉DNA，DNAon−k→3 mRNA+ DNA开启。上述一组反应描述了DNA通过聚合酶结合/解结合和聚合酶结合（即，启动）DNA转录成mRNA。开关由第一个反应将DNA关闭转换为DNA打开来描述，开关由第二个反应将DNA关闭转换为DNA打开来描述，转录是由于第三个反应产生mRNA而使DNA的实际量保持不变。 当涉及三种试剂时，系统的状态是三重x = |x1，x2，x3|分别描述了DNA关闭、DNA开启和mRNA的量。描述三个反应的效应的向量为：|-1，1，0|，e2= |1，−1，0 |e 3 =| 0，0，1|.假设质量作用动力学，与反应相关的强度为：f1（x）=k1x1，f2（x）=k2x2和f3（x）=k3x2。基于以上描述，系统的随机模拟是一个简单的任务。3 EM方法期望最大化（EM）方法是在不完全数据的情况下用于最大似然参数估计的算法。EM算法的输入由一组样本和一个随机模型组成，该随机模型由一组用λ表示的参数表征。EM方法是迭代的，即，它从参数集λ0的初始猜测开始，然后以这样一种方式逐步改进它，即模型提供的行为越来越类似于样本描述的行为。i步后的参数集用λi表示。每次迭代由两个步骤组成，分别称为期望步骤（E-步骤）和最大化步骤（M-步骤）。26A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）23E步骤的作用是计算期望中的缺失信息。形式上，用Y表示不完整数据集，用Z表示完整数据集，E步骤计算条件期望E [Z|Y，λ i]。在我们的问题的上下文中，Y包含离散时间点的样本，E步骤的目的是计算，给定λi中的当前参数集，最典型的完整轨迹Z，该轨迹通过Y中给定的观察样本。然后应用M步来找到一组新的参数λi+1，使得在E步期间生成的轨迹的可能性最大。一旦找到新的参数集λi+1，它就被用作下一次迭代的起点。在许多情况下，包括本文所考虑的情况，的EM方法在于这样一个事实，即找到这样的参数，最大化不完整数据的可能性比找到使完整数据的可能性最大化的参数要困难得多。换句话说，M步骤中所需的优化比原始优化问题的负担更小。特别是，在本文中考虑的问题，M-步骤是非常简单的情况下的质量作用动力学。反过来，在许多情况下，包括我们的情况，从理论和计算的角度来看，E步骤可能很难对于其中E步骤所需的计算特别复杂的那些情况，可以应用EM方法的变体。在这个变体中，E步骤中条件期望的精确计算被模拟代替。这种方法被称为蒙特卡罗EM（MCEM）方法，并且在以精确方式执行E步骤太耗时或甚至不可行的情况下特别有用。对于本文所考虑的问题，由于所涉及的CTMC的巨大的状态空间，唯一可行的方法是由MCEM方法提供MCEM方法的收敛特性比EM方法差，并且模拟会引入参数的波动。 然而，如果迭代次数很高，仍然可以证明该方法的收敛性得到保持[15]。4用于生化系统的问题公式化我们假设我们得到一个反应网络和在离散时刻对所涉及物种数量的实验 时间瞬间的观测值表示为t0 = 0，t1，t2，.，t N和由y0，y1，.，其中yi是提供系统在ti处的状态的整数向量。(We这里仅考虑单个观测序列，但是扩展到多个观测序列是直接的。我们假设所有或一些动力学速率是未知的，即，在提供反应强度的函数fi（x）中存在未知常数。这些未知常数的集合将用λ表示，我们将写fi（x，λ）以明确强度对未知数的依赖性。我们的目标是通过MCEM方法给出未知动力学速率的极大似然估计在A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）2327在下面的两个小节中，我们描述E步骤和M步骤。e步骤给定样本集和参数的当前估计，E步骤旨在构建经过所观察到的状态的最可能的完整轨迹。这需要找到每两个连续样本点之间的最可能轨迹。因此，E步需要找到CTMC上最可能的随机游走。这在理论上是可能的，但不幸的是，由于在我们的上下文中，所考虑的CTMC几乎总是具有巨大的状态空间，因此无法执行以精确的方式。如所预期的，在这种情况下，E步可以通过模拟来求解，这提供了最可能轨迹的良好近似。通过仿真生成CTMC的轨迹是简单的。由于我们被给予N+ 1个样本，我们需要N个子迹来连接观察时刻。由于CTMC是巨大的，并且当前的参数集可以是远离真实的参数集，单个模拟运行不太可能精确地（或甚至接近）到达所观察的状态。因此，E步骤由以下两个阶段组成。(i) 生成痕迹。 对于每个区间[t i，t i+1]，0 ≤ i ≤ N − 1，我们从yi开始生成K个长度为t i+1−t i的随机游动，并选择到达“距离”最接近y i +1的随机游动。样本点和随机游走最后一个状态之间的“距离”的概念(ii) 改善跟踪。在这个阶段，我们通过修改子迹来改进它们。我们从子迹中随机选取一个反应，并检查它是否可以被另一个反应取代，以使子迹更接近观察状态。只有当所有剩余的反应仍然是可能的时，这种替代才被接受.连续反应之间的时间保持不变。修改的程度由参数ρ ∈ [0，1]决定，该参数定义了我们试图替换的反应的比例在估算过程中，ρ所定义的比例起着微妙的作用。当当前估计远离真实参数时，需要更高的ρ，以便提出合理的子迹并使估计过程具有更快的收敛性。相反，当估计已经很好时，必须使用较低的ρ（甚至ρ= 0），以便不太改变由实际估计引起的随机行为。在复杂性方面，随机游走的生成成本与反应发生的次数成线性关系，并且最佳随机游走的存储非常便宜，因为每个子迹都可以通过伪随机数生成器的种子来唯一识别。由于这些原因，最佳迹线的改进代表了该方法的最昂贵的阶段，因为每次试验都迫使最佳迹线的改进。28A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）23Σ⎛Σ展开并检查所有后续反应的启动3。因此，这一阶段的实施可能并非微不足道。第一个提示可以是使用预处理，在随机选择用于替换的“候选者”之后，将它们按出现顺序排序。以这种方式，整个跟踪可以被展开，并且对于每次尝试，检查的数量变得越来越小。第二个提示是一次验证多个反应，如果控制不满意，则丢弃（不增加成本）所有替代品。这种解决方案可能会丢弃一些有效的替代，但在计算时间方面更方便m步在M步中，我们必须找到下一组估计值λi+1，它使E步中生成的完整轨迹由E步生成的第i个子迹，在时间i_int val[ti-1，ti]中重建状态y i-1和状态y i之间的最可能轨迹，将被表示为ySi，并且它的形式为yi−1 =si，1ri，1，ui，1−→s我2ri，2，ui，2−→··· −→ui，Hisi，Hi−→s i，Hi其中Hi表示第i个子迹的长度，si，j，ri，j和ui，j分别是第i个子迹的状态、反应和停留时间我们有这个Hiui，j= ti− ti− 1。j=1最后一个箭头表示没有反应，这表示过程在至少ui，Hi时间单位内保持在状态si，Hi状态si，Hi越接近状态yi，第i个子迹越好地反映所观察到的行为。从CTMC的理论可以得出，第i个子迹的似然性（由Li表示）可以计算为：Hi−1Li=j=1fri，j（si，j，λ）e−f（si，j，λ）ui，j·e−f（si，Hi，λ）ui，Hi（2）哪里Rf（x，λ）=fk（x，λ）k=1是状态x下反应强度的总和。（2）中的乘积给出了第i个子迹的转移的可能性，而最后一个指数项是过程至少在ui，Hi次内不离开状态si，Hi的概率3请注意，迹线的展开几乎与其生成一样昂贵⎝A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）2329NN⎛ΣΣN.ΣΣ M.ΣNΣ MΣΣ.拉乌=i=1。i、j单位 所有子迹的可能性简单地由下式给出：NL=Lii=1我们必须找到使乘积最大化的λ为了找到L的最大值，取其对数是有用的，在对数中，乘积被转换为和，ln（L）=i=1ln（Li）=i=1Hi−1j=1ln（fri，j（si，j，λ））−Hij=1f（si，j，λ）ui，j≠0.（三）自然地，找到（3）的最大值的难度取决于函数fi（x，λ），1≤i≤R。如前所述，我们考虑反应强度对应于质量作用动力学（1）的情况此外，我们假设化学计量的反应（描述为ai，j，1≤i≤R， 1≤j≤M）是已知的，这是参数估计问题的典型情况。因此，要估计的参数是动力学速率常数，即，λ ={k1，...，k R}。在不失一般性的情况下，我们将注意力集中在找到使（3）最大化的k1上应用（1），（3）关于k1的导数为：白蛋白（L） =1号线Σ⎛HΣi−1Mk=11xka1，k你好 。 xi=1j=1I{ri，j是反应1}kMk=1xk−a1、kj=1k=1a1、kui，j其中，如果参数为真，则I为1，否则为0用fi，j表示数字反应j发生在第i个子迹中的次数，白蛋白（L）1996年，你好 。 xK1号线=i=1 1号线−j=1k=1a1、ku i，ji。（四）很容易检验，当k1的值为0时，（4）使L最大化。因此，估计数为nfi，1因此，在质量作用动力学的情况下，M步骤中所需的优化归结为（5）中给出的显式公式处理fi（x，λ）的其他形式超出了本文的范围。我们在这里只提到，对于一般的fi（x，λ）函数，M步所需的优化可能变得更复杂，但即使在这种情况下，也可以将K1Ni=1Hij=1Mk=1XKa1、kK30A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）23优化问题划分为更小的子问题。为了说明这一点，让我们用νi表示fi（x，λ）所依赖的变量可以合理地假设，A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）2331集合νi，1≤i≤R，相互不相交。在这种情况下，对数似然函数（在（3）中给出）根据属于νi的变量的导数不取决于属于其他集合的变量vj，j i。 这意味着涉及所有变量的原始优化问题可以通过R更小的优化问题，更小的维度。5说明性数值示例在本节中，我们将展示使用MCEM方法获得的数值结果。我们将该方法应用于两个模型。对于这两种情况，通过模拟“在计算机上”生成样本在这些测试中，我们把模型的生物意义放在一边，我们的目的是说明该方法，并表明它能够重建参数集。对于所有情况，我们提供表格来比较原始值及其估计值。此外，为了直观地比较原始值和估计值的行为，我们根据相应的常微分方程（ODE）和相应的随机设置来确定系统的演化常微分方程有助于快速了解在随机环境中获得的估计值的优劣。MCEM方法已在原型JAVA工具中实现。所有实验都是在具有4Gb RAM的Intel迅驰双核上进行的5.1基因转录模型我们的第一个例子是在第2节中已经介绍的模型，描述了DNA的结合和解结合及其转录成mRNA[6]。我们假设具有单个DNA单元，并且初始条件是[DNAoff]= 1和[DNAon] =[mRNA] = 0。为了在不同情况下评估该方法，我们使用具有不同“粒度”水平的模型我们使用不同水平的离散化，以获得离散模型从原来的连续浓度。离散化步骤将由h表示。初始状态的CTMC建模的三个反应，联系我们|1/h，0，0|. 为了让模型在同一时间尺度上发展，与h无关，动力学速率常数必须取决于h。具体地说，反应的强度是k1h，k2h和k3h，其中k1 = 0。027，k2 = 0. 166且k3 = 0。4. h的效应是双重的：h越小，状态空间越大，模型行为的变量越少。事实上，当h趋于0时，模型的行为趋于相应的常微分方程的解[8，14]。我们已经产生了“在硅片”样本，在这样一种方式，连续采样点之间有大约25000个反应。样本数量为25个。这意味着模型中发生的事情只有一小部分可用于估计过程。在K= 20的情况下进行E步骤，即，在E步骤的第一阶段生成了20条迹线。在第二阶段，参数ρ起着至关重要的作用。从一个高的ρ值开始，然后随着时间的推移逐渐降低，32A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）23估计变得更加可靠。我们选择从ρ = 0开始。4，并随着生成的迹线越来越接近观测值而减小。表1报告了原始参数和在大约一秒的CPU时间内计算的100次迭代后通过MCEM方法获得的参数。参数的初始猜测在[0：10]范围内是随机的。关于结果的一个可能的观察可能是，k1和k2的一些估计远离原始值。这是由于少数不常见的离散样本反映了比值 k1/（k1+k2）×k3（决定mRNA的增加），在某种程度上反映了比值k1/（k1+k2）（决定DNA的量），而不是参数的值。然而，具有更高采样频率的更多样本可以缓解这个问题。还可以观察到，随着更精细的离散化，结果更准确。在同一张图上说明所有情况的一种方法是将结果应用于代表模型的常微分方程中。这在图1中描述。可以看出，所有情况都很好地捕捉到mRNA的渐近增加率，并且h =0。估计值很好地再现了原始模型。常微分方程不能提供一个准确的表示CTMC轨迹，因为这个原因，我们也考虑了随机设置。图2描绘了100的平均值和方差。000条模拟轨迹。作为最后的观察，重要的是要指出，通过使用较小的h值，该方法变得更好。K1K2K3原始0.0270.1660.4估计值，h=10.65951.53860.1879估计值，h=0.10.12790.63850.3431估计值，h=0.010.10980.66830.4115估计值，h=0.0010.03970.25910.4045表1基因转录模型的参数估计结果5.2DFG降解途径为了测试具有更多变量的方法，我们使用描述N-（脱氧-D-果糖-1-基）-甘氨酸（DFG）降解途径控制的模型[9]。该模型可在www.sbml.org网站上的数据库中找到。它涉及14种试剂通过16个反应相互作用，具有质量作用动力学（反应报告见表2）。将原ODE模型转化为离散步长h = 0的CTMC。0001，并从初始浓度[DFG] =9和所有其他试剂的量0开始。注意，h的这种选择导致具有巨大状态空间的CTMC。样本包含20个观察结果，连续样本之间约有10000个反应。在K= 20和ρ从0.2开始并在整个计算过程中降低的情况下进行E步骤。结果在表3中给出。可以看出，该方法给出了很好的估计几乎所有涉及的参数。图5描绘了根据具有原始和估计参数的ODE的相应系统的模型的一些试剂的演变。图4中描述了随机设置，为了便于综合，我们仅报告方差为A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）2333原始h=1h=0.1h=0.01h=0.0011 0.350.95 0.30.9 0.250.85 0.20.80.750.70.650 2 4 6 8100.60.50.150.10.050原始h=1h=0.1h=0.01h =0.0010 2 4 6 8 100.40.30.20.100 2 4 6 8 10Fig. 1.基因转录模型：用原始参数和不同h值的估计值计算的ODE（左边是DNA，右边 是DNA，下面是mRNA大多数物种的原始行为与估计值提供的行为之间存在良好的对应关系在图3中，我们描述了估计值重现原始“计算机模拟”生成轨迹的可能性经过350次交互后，估计值的似然性非常接近原始参数的似然性反应K1DFG → E1DFGk2E2→DFGk3Gly+Cn→E1k4Gly +DG 3→K5DG3→ CnDG3k6wF A→E2k7Gly +DG 1→DG1k8Cn→K9DG1→ AAE1k10Gly+Man→E1k11Gly +Glu→Mank12 Glu→k13Glu→ DG3Gly+Ck14Meln →Ck15AA+FA+MGn →Ek16Gly+Fru2→表2DFG降解途径壳体3K1K2K3K4K5K6K7K8原始结果0.0050.00480.0150.01750.0150.01190.0790.0580.0900.0680.0270.0100.2120.2520.1810.706情况K9k10K11K12k13k14k15k16原始结果1.9081.8470.0700.06510.1130.1228.0E-47.8E-40.0020.0140.0030.0030.0150.01470.0130.0122表3DFG降解途径的参数估计结果原始h=1 h=0.1h=0.01h=0.00134A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）23原始h=1h=0.1h=0.01h=0.001原始h=1h=0.1h=0.01h=0.00110.250.950.90.850.80.750.70.650 2 4 6 8100.20.150.10.0500 2 4 6 8 100.35 0.250.30.250.20.150.10.0500 2 4 6 8100.20.150.10.0500 2 4 6 8 100.610.90.50.80.40.70.60.30.50.40.2原始0.3原始0.1h=1h=0.10.2h=1h=0.1h=0.01h=0.0010.1h=0.01h=0.001000 2 4 6 8 10 0 2 4 6 8 10图二、基因转录模型：DNA的数量的平均值（左）和变异（右），用原始参数和h的不同值的估计计算的mRna。-250000-300000-350000-400000-450000-500000-550000-600000-650000-700000-750000结果可能样本0 50 100150200250300 350图3.第三章。DFG降解途径：似然估计值作为迭代次数的函数原始h=1h=0.1h=0.01h=0.001原始h=1h=0.1h=0.01h=0.001A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）23352.521.510.50原始AA结果AA0 20 40 60 801001.61.41.210.80.60.40.200 204060 80 1002.521.510.50原始DFG结果DFG0 20 40 60 801000.50.450.40.350.30.250.20.150.10.050原始Glu结果Glu0 20 40 60 80 1002.521.510.5原始Gly1.61.41.210.80.60.40.2原始梅尔0结果Gly0 20 40 60 8010000 2040结果Mel60 80 1001.61.41.210.80.60.40.20原始MG结果MG0 20 40 60 80 100图四、参与DFG降解途径的一些试剂的变化结论本文将MCEM方法应用于随机反应网络中动力学速率的估计。我们已经表明，由此产生的技术是有效的，并导致令人惊讶的简单计算的情况下，质量作用动力学。所提出的方法的优点是，它可以适用于即使是一个有限的集原始Cn结果Cn36A. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）2398765432100 204060 801000.140.120.10.080.060.040.020原始Cn结果Cn0 20 40 60 80 1000.450.40.350.30.250.20.150.10.0500 204060 801000.50.450.40.350.30.250.20.150.10.050原始Glu结果Glu0 20 40 60 80 1000.35 0.250.30.250.20.150.10.050原始Gly结果Gly0 20 40 60 801000.20.150.10.050原始Mel结果Mel0 20 40 60 80 10021.81.61.41.210.80.60.40.20原始MG结果MG0 20 40 60 80 100图五、DFG降解途径中涉及的一些试剂的ODE对模拟现象的观察几个数值例子说明了该方法的计算特性原始AA结果AA原始DFG结果DFGA. Angius，A.Horváth/Electronic Notes in Theoretical Computer Science 275（2011）2337引用[1] 博 伊 斯 河 J. ， D. Wilkinson 和 T. B. L. Kirkwood ， Bayesian inference for a discrete observedstochastic kinetic model，Statistics and Computing18（2008），pp. 125比135[2] 巴罗斯河，G. Warnes和R. ChoughanHanumara，生化途径的统计建模，技术报告06/11，部门。生物统计学和计算生物学，罗切斯特大学（2006年）。[3] Dempster，A. P.，N. M. Laird和D. B. Rubin，通过em算法从不完整数据中获得最大似然，皇家统计学会杂志，B系列39（1977），pp.1-38[4] 加 德 卡 湾G.， R. Gunawan 和 F. J. Doyle 3rd， Iterative approach to model identification of biologicalnetworks，BMC Bioinformatics6（2005）。[5] Gillespie，D. T.，耦合化学反应的精确随机模拟，J. Phys. Chem.81（1977），pp. 2340-2361。[6] 戈 尔丁 岛， J. Paulsson ， S. M. Zawilski和E. C. Cox ， Real-time kinetics of gene activity in individualbacteria，Cell123（2005）.1025-1036。[7] Golightly ， A. 和 D. Wilkinson ， Bayesian inference for stochastic kinetic models using a difficultapproximation，Biometrics61（2005），pp.781-788[8] Kurtz ， T. G. ， Solutions of ordinary diffusion equations as limits of pure jump Markov processes ，Journal of Applied Probability1（1970），pp.49比58[9] 马 丁 斯 ， S. 一 、 A. T. Martinus 和 M. A. V. Boekel ， Amadori n- （ 1-deoxy-fructos-1- yl ） -glycinedegradation pathways的动力学建模。第ii部分-动力学分析，碳水化合物研究338（2003），pp. 1665- 1678年。[10] Moles 角 ， P. Mendes 和 J. Banga ， 生 化 途 径 中 的 参 数 估 计 ： 全 局 优 化 方 法 的 比 较 ， Genome Res.13（2003），pp. 2467-2474[11] Reinker ， S. ， R. Altman 和 J. Timmer ， 随 机 化 学 反 应 中 的 参 数 估 计 ， IEEE Proceedings SystemsBiology153（2006），pp. 168比178[12] Sugimoto，M.，S. Kikuchi和M. Tomita，Reverse engineering of biochemical equations from time-coursedata by means of genetic programming，Biosystems80（2005），pp. 155-164。[13] 田，T.， S. Xu， J. Gao和K. Burrage， 模拟最大似然法估计基因表达的动力学速率，生物信息学23（2007），pp. 84比91[14] Tribastone，M.，“随机过程代数模型的可扩展分析”，博士。论文，爱丁堡大学信息学院（2010年）。[15] Wei，G. C. G.和M. A. Tanner，em算法和穷人数据增强算法的蒙特卡洛实现699-704
我的内容管理
收起
我的收益 登录查看自己的收益
我的积分
登录查看自己的积分
我的C币
登录后查看C币余额
我的收藏 
我的下载 
下载帮助 
