ギレスピーアルゴリズム

確率論において、ギレスピーアルゴリズム（またはドゥーブ–ギレスピーアルゴリズム、確率的シミュレーションアルゴリズム、SSA ）は、反応速度がわかっている確率方程式システムの統計的に正しい軌道（可能な解）を生成します。これはJoseph L. Doobら（1945年頃）によって考案され、 Dan Gillespieが1976年に発表し、1977年に限られた計算能力を使用して化学または生化学の反応システムを効率的かつ正確にシミュレートするために使用した論文で普及しました（確率的シミュレーションを参照）。^[1]コンピュータの速度が速くなるにつれて、このアルゴリズムはますます複雑なシステムをシミュレートするために使用されるようになりました。このアルゴリズムは、試薬の数が少なく、すべての反応を追跡することが計算的に実行可能な、細胞内の反応をシミュレートするのに特に便利です。数学的には、動的モンテカルロ法の変形であり、運動論的モンテカルロ法に似ています。計算システム生物学で頻繁に使用されます。^[要出典]

このアルゴリズムに至るプロセスには、いくつかの重要なステップがあります。1931 年、アンドレイ・コルモゴロフは、ジャンプによって進行する確率過程の時間発展に対応する微分方程式を導入しました。これは今日ではコルモゴロフ方程式 (マルコフジャンプ過程)として知られています (自然科学では、簡略化されたバージョンがマスター方程式として知られています) 。1940 年に、コルモゴロフ方程式が (適切な) 確率を解として受け入れる条件を発見したのは、ウィリアム・フェラーでした。彼は、定理 I (1940 年の作業) で、次のジャンプまでの時間は指数分布し、次のイベントの確率はその率に比例することを立証しました。このようにして、彼はコルモゴロフ方程式と確率過程の関係を確立しました。その後、ドゥーブ (1942、1945) は、フェラーの解を純粋なジャンプ過程の場合を超えて拡張しました。この方法は、 David George Kendall (1950)がManchester Mark 1コンピュータを使用してコンピュータに実装し、後にMaurice S. Bartlett (1953) が伝染病の発生に関する研究で使用しました。Gillespie (1977) は、物理的な議論を利用して別の方法でアルゴリズムを取得しました。

反応室内には、有限数の分子が存在します。極小の時間ごとに、1 つの反応が起こる可能性があります。反応速度は、各化学種の分子数によって決まります。

単純に言えば、時間を離散化して反応室の軌道をシミュレートし、各タイムステップをシミュレートすることができます。ただし、反応が発生しない時間が長く続く可能性があります。Gillespie アルゴリズムは、何らかの反応が発生するまでランダムに待機時間をサンプリングし、その後別のランダム サンプルを取得してどの反応が発生したか を判断します。

重要な前提は

• それぞれの反応は時間的にマルコフ的である
• 反応間に相関関係はない

2 つの仮定を前提とすると、ある反応のランダムな待機時間は指数分布し、指数率は個々の反応の速度の合計になります。

従来の連続的かつ決定論的な生化学反応速度式は、数百万の分子の相互作用を必要とするバルク反応に依存しているため、細胞反応を正確に予測できません。通常、これらは結合した常微分方程式のセットとしてモデル化されます。対照的に、Gillespie アルゴリズムでは、すべての反応が明示的にシミュレートされるため、反応物の少ないシステムの離散的かつ確率的なシミュレーションが可能です。単一の Gillespie シミュレーションに対応する軌跡は、マスター方程式の解である確率質量関数からの正確なサンプルを表します。

アルゴリズムの物理的基礎は、反応容器内での分子の衝突です。衝突は頻繁に起こると想定されますが、適切な方向とエネルギーでの衝突はまれです。反応環境は十分に混合されていると想定されます。

レビュー（Gillespie、2007）では、直接法、第一反応法、第一ファミリー法という3つの異なるが同等の定式化が概説されており、前者2つは後者の特別なケースである。直接法と第一反応法の定式化は、いわゆる「確率的化学反応速度論の基本前提」に基づいて通常のモンテカルロ逆変換ステップを実行することに重点を置いており、数学的には関数である。

${\displaystyle p(\tau,j\mid{\boldsymbol{x}},t)=a_{j}({\boldsymbol{x}})\exp\left(-\tau\sum_{j}a_{j}({\boldsymbol{x}})\right),}$

ここで、各項は素反応の傾向関数であり、その引数は、種の数のベクトルです。パラメータは次の反応までの時間（または滞在時間）で、 は現在の時間です。Gillespie を言い換えると、この式は「 が与えられた場合に、システムの次の反応が無限小の時間間隔 で発生し、番目の反応に対応する化学量論になる確率」と解釈されます。この定式化により、 は指数分布のランダム変数であり、 は「点確率を持つ統計的に独立した整数ランダム変数」 であることが示唆され、直接法と最初の反応法への窓口が提供されます。${\displaystyle a}$${\displaystyle {\boldsymbol {x}}}$${\displaystyle \tau}$${\displaystyle t}$${\displaystyle {\boldsymbol {X}}(t)={\boldsymbol {x}}}$${\displaystyle [t+\tau ,t+\tau +d\tau ]}$${\displaystyle j}$${\displaystyle \tau}$${\displaystyle j}$${\displaystyle a_{j}({\boldsymbol {x}})/\sum _{j}a_{j}({\boldsymbol {x}})}$

したがって、モンテカルロ法は、2つの擬似乱数とを上に描画し、計算する だけです。${\displaystyle r_{1}}$${\displaystyle r_{2}}$${\displaystyle [0,1]}$

${\displaystyle \tau ={\frac {1}{\sum _{j}a_{j}({\boldsymbol {x}})}}\log \left({\frac {1}{r_{1}}}\right),}$

そして

${\displaystyle j={}}$を満たす最小の整数${\displaystyle \sum _{j'=1}^{j}a_{j'}({\boldsymbol {x}})>r_{2}\sum _{j}a_{j}({\boldsymbol {x}}).}$

この生成法を滞在時間と次の反応に利用して、直接法アルゴリズムはギレスピーによって次のように述べられている。

1. 時間とシステムの状態を初期化します。 
2. システムが時刻 の状態にある状態で、すべてのとそれらの合計を評価します。 3.および
の上記の値を計算します。
4. および を置き換えて次の反応を実行します。
5.必要に応じて記録します。ステップ 2 に戻るか、シミュレーションを終了します。
${\displaystyle t=t_{0}}$${\displaystyle {\boldsymbol {x}}={\boldsymbol {x}}_{0}}$${\displaystyle {\boldsymbol {x}}}$${\displaystyle t}$${\displaystyle a_{j}({\boldsymbol {x}})}$${\displaystyle \sum _{j}a_{j}({\boldsymbol {x}})}$${\displaystyle \tau}$${\displaystyle j}$${\displaystyle t\leftarrow t+\tau }$${\displaystyle {\boldsymbol {x}}\leftarrow {\boldsymbol {x}}+\nu _{j}}$${\displaystyle ({\boldsymbol {x}},t)}$

ここで、 は、与えられた状態変化ベクトル の成分を追加することを表します。このアルゴリズム ファミリは計算コストが高いため、次の反応法 (Gibson & Bruck)、タウ リーピング、および豊富な反応物が決定論的な動作でモデル化されるハイブリッド手法など、多くの修正と適応が存在します。適応された手法は、マスター方程式に接続するため、アルゴリズムの背後にある理論の正確さを一般に損ないますが、大幅に改善された時間スケールで妥当な実現を提供します。アルゴリズムの正確なバージョンの計算コストは​​、反応ネットワークの結合クラスによって決まります。弱く結合されたネットワークでは、他の反応の影響を受ける反応の数は、小さな定数によって制限されます。強く結合されたネットワークでは、1 つの反応の発火が原理的に他のすべての反応に影響を与える可能性があります。弱く結合されたネットワーク用の定数時間スケーリングを備えたアルゴリズムの正確なバージョンが開発され、非常に多数の反応チャネルを持つシステムの効率的なシミュレーションが可能になりました (Slepoy Thompson Plimpton 2008)。遅延を伴うランダムな生化学イベントの非マルコフ特性を考慮した一般化 Gillespie アルゴリズムは、Bratsun ら (2005) およびそれとは独立して Barrio ら (2006)、さらに (Cai 2007) によって開発されました。詳細については、以下の引用文献を参照してください。 ${\displaystyle \nu _{j}}$${\displaystyle j^{\text{th}}}$${\displaystyle \nu }$

Ramaswamy ら (2009、2010) と Indurkhya と Beal (2010) の両者が独立して開発した部分的傾向定式化は、計算コストが (より大きな) 反応数ではなく、ネットワーク内の化学種の数に比例するアルゴリズムの正確なバージョンのファミリーを構築するために使用できます。これらの定式化により、計算コストを弱く結合したネットワークでは定数時間スケーリングに、強く結合したネットワークでは種の数に最大限線形にスケーリングできます。遅延のある反応に対する一般化 Gillespie アルゴリズムの部分的傾向バリアントも提案されています (Ramaswamy Sbalzarini 2011)。部分的傾向法の使用は、基本化学反応、つまり最大で 2 つの異なる反応物による反応に限定されます。すべての非基本化学反応は、ネットワーク サイズが線形 (反応の順序で) 増加することを犠牲にして、基本化学反応のセットに等価的に分解できます。

簡単な例を挙げると、Gillespie アルゴリズムの仕組みが説明しやすくなります。2 種類の分子AとBのシステムを考えてみましょう。このシステムでは、AとB が可逆的に結合してABダイマーを形成します。そのため、2 つの反応が可能です。A と B が可逆的に反応してABダイマーを形成するか、ABダイマーがAとBに解離します。特定の単一の A 分子が特定の単一のB分子と反応する場合の反応速度定数は で、 ABダイマーが分解する場合の反応速度は です。 ${\displaystyle k_{\mathrm {D} }}$${\displaystyle k_{\mathrm {B} }}$

時刻tに各タイプの分子が 1 つずつ存在する場合、二量体形成速度は ですが、タイプAの分子とタイプBの分子が存在する場合、二量体形成速度は です。二量体が存在する場合、二量体の解離速度は です。 ${\displaystyle k_{\mathrm {D} }}$${\displaystyle n_{\mathrm {A} }}$${\displaystyle n_{\mathrm {B} }}$${\displaystyle k_{\mathrm {D} }n_{\mathrm {A} }n_{\mathrm {B} }}$${\displaystyle n_{\mathrm {AB} }}$${\displaystyle k_{\mathrm {B} }n_{\mathrm {AB} }}$

時刻tにおける全反応速度は次のように表される。 ${\displaystyle R_{\mathrm {TOT} }}$

${\displaystyle R_{\mathrm {TOT} }=k_{\mathrm {D} }n_{\mathrm {A} }n_{\mathrm {B} }+k_{\mathrm {B} }n_{\mathrm {AB} }}$

これで、2 つの反応を持つ単純なモデルについて説明しました。この定義は、Gillespie アルゴリズムとは独立しています。次に、このシステムに Gillespie アルゴリズムを適用する方法について説明します。

このアルゴリズムでは、次の反応までの時間を計算し、次の反応が可能な反応のどれであるかを決定するという 2 つのステップで時間を進めます。反応は完全にランダムであると想定されるため、時間tにおける反応率が である場合、次の反応が発生するまでの時間 δ tは、平均 の指数分布関数から抽出された乱数です。したがって、時間をtからt + δ tに進めます。 ${\displaystyle R_{\mathrm {TOT} }}$${\displaystyle 1/R_{\mathrm {TOT} }}$

この反応がA分子とB分子の結合である確率は、単にこのタイプの反応による全速度の割合である。すなわち、

反応が起こる確率は${\displaystyle P({\ce {{A}+ B -> AB}})=k_{D}n_{A}n_{B}/R_{{\ce {TOT}}}}$

次の反応でAB二量体が解離する確率は、その値から 1 を引いた値です。したがって、この 2 つの確率では、とを 1 減らして 1 増やすことによって二量体を形成するか、または二量体を解離して と を1増やして1 減らすことになります。 ${\displaystyle n_{\mathrm {A} }}$${\displaystyle n_{\mathrm {B} }}$${\displaystyle n_{\mathrm {AB} }}$${\displaystyle n_{\mathrm {A} }}$${\displaystyle n_{\mathrm {B} }}$${\displaystyle n_{\mathrm {AB} }}$

これで、時間をt + δ tまで進め、単一の反応を実行できました。Gillespie アルゴリズムは、システムを必要な時間だけシミュレートするために必要な回数だけ、この 2 つのステップを繰り返します (つまり、反応の数だけ)。t =0でおよびから始まり、および である Gillespie シミュレーションの結果が右側に示されています。これらのパラメーター値では、平均して 8 つの二量体とAおよびBが 2 つありますが、分子の数が少ないため、これらの値の変動は大きくなります。Gillespie アルゴリズムは、これらの変動が重要なシステムを研究するためによく使用されます。 ${\displaystyle n_{\mathrm {A} }=n_{\mathrm {B} }=10}$${\displaystyle n_{\mathrm {AB} }=0}$${\displaystyle k_{\mathrm {D} }=2}$${\displaystyle k_{\mathrm {B} }=1}$${\displaystyle n_{\mathrm {AB} }}$

これは、2 つの反応がある単純な例です。より多くの反応があるより複雑なシステムも、同様に扱われます。すべての反応速度は、各時間ステップで計算され、速度への寄与率に等しい確率で 1 つが選択されます。その後、この例のように時間が進められます。

Gard モデルは、脂質が凝集体へと自己組織化する様子を説明しています。確率論的シミュレーションを使用して、複数の種類の凝集体の出現とその進化を示します。

• Gillespie, Daniel T. (1977). 「結合化学反応の正確な確率シミュレーション」. The Journal of Physical Chemistry . 81 (25): 2340–2361. CiteSeerX  10.1.1.704.7634 . doi :10.1021/j100540a008. S2CID  2606191.
• Gillespie, Daniel T. (1976). 「結合化学反応の確率的時間発展を数値的にシミュレートする一般的な方法」。Journal of Computational Physics . 22 (4): 403–434. Bibcode :1976JCoPh..22..403G. doi :10.1016/0021-9991(76)90041-3.
• Gibson, Michael A.; Bruck, Jehoshua (2000). 「多くの種と多くのチャネルを持つ化学システムの効率的で正確な確率的シミュレーション」(PDF) . Journal of Physical Chemistry A . 104 (9): 1876–1889. Bibcode :2000JPCA..104.1876G. doi :10.1021/jp993732q.
• Doob, Jacob L. (1942). 「マルコフ連鎖理論の話題」.アメリカ数学会誌. 52 (1): 37–64. doi : 10.1090/S0002-9947-1942-0006633-7 . JSTOR  1990152.
• Doob, Jacob L. (1945). 「マルコフ連鎖 - 可算なケース」.アメリカ数学会誌. 58 (3): 455–473. doi :10.2307/1990339. JSTOR  1990339.
• Press, William H.; Teukolsky, Saul A.; Vetterling, William T.; Flannery, Brian P. (2007)。「セクション 17.7. 化学反応ネットワークの確率的シミュレーション」。数値レシピ: 科学計算の技法(第 3 版)。ニューヨーク、NY: Cambridge University Press。ISBN 978-0-521-88068-8. 2011年8月11日時点のオリジナルよりアーカイブ。2011年8月17日閲覧。
• コルモゴロフ、アンドレイ N. (1931)。 「確率論における分析方法について」。数学アンナレン。104 : 415–458。土井：10.1007/BF01457949。S2CID  119439925。
• フェラー、ウィリー(1940)。 「純粋に不連続なマルコフ過程の積分微分方程式について」。アメリカ数学会誌。48 (3): 4885–15。doi : 10.2307/1990095。JSTOR  1970064 。
• ケンドール、デイビッド G. (1950)。「単純な「誕生と死」のプロセスの人工的実現」。王立統計学会誌、シリーズ B。12 ( 1): 116–119。JSTOR 2983837  。
• バートレット、モーリス・S. ( 1953) 。「確率過程または変化の統計」。王立統計学会誌、シリーズ C。2 ( 1 ): 44–64。doi :10.2307/2985327。JSTOR  2985327。
• Rathinam, Muruhan; Petzold, Linda R.; Cao, Yang; Gillespie, Daniel T. (2003). 「確率的化学反応システムの剛性: 暗黙のタウリーピング法」. Journal of Chemical Physics . 119 (24): 12784–12794. Bibcode :2003JChPh.11912784R. doi :10.1063/1.1627296.
• Sinitsyn , Nikolai A.; Hengartner, Nicolas; Nemenman, Ilya (2009). 「確率的生化学ネットワークの断熱粗視化とシミュレーション」。米国科学アカデミー紀要。106 (20): 10546–10551。Bibcode : 2009PNAS..10610546S。doi :  10.1073 / pnas.0809340106。PMC 2705573。PMID 19525397。 
• Salis, Howard; Kaznessis, Yiannis N. (2005). 「結合化学反応または生化学反応システムの正確なハイブリッド確率シミュレーション」。Journal of Chemical Physics . 122 (5): 054103. Bibcode :2005JChPh.122e4103S. doi :10.1063/1.1835951. PMID  15740306.
• (Slepoy Thompson Plimpton 2008): Slepoy, Alexander ; Thompson, Aidan P.; Plimpton, Steven J. (2008)。「大規模生化学反応ネットワークのシミュレーションのための定数時間運動モンテカルロアルゴリズム」。Journal of Chemical Physics。128 ( 20): 205101。Bibcode :2008JChPh.128t5101S。doi :10.1063/1.2919546。PMID 18513044  。
• (Bratsun et al. 2005): Bratsun, Dmitri; Volfson, Dmitri; Hasty, Jeff; Tsimring, Lev S. (2005). 「遺伝子調節における遅延誘導性確率振動」。米国科学アカデミー紀要。102 (41): 14593–8。Bibcode : 2005PNAS..10214593B。doi :  10.1073 / pnas.0503858102。PMC 1253555。PMID 16199522。 
• (Barrio et al. 2006): Barrio, Manuel; Burrage, Kevin; Leier, André; Tian, Tianhai (2006). 「hes1 の振動制御: 離散確率的遅延モデリングとシミュレーション」. PLOS Computational Biology . 2 (9): 1017. Bibcode :2006PLSCB...2..117B. doi : 10.1371 /journal.pcbi.0020117 . PMC 1560403. PMID  16965175. 
• (Cai 2007): Cai, Xiaodong (2007). 「遅延を伴う結合化学反応の正確な確率シミュレーション」. Journal of Chemical Physics . 126 (12): 124108. Bibcode :2007JChPh.126l4108C. doi :10.1063/1.2710253. PMID  17411109.
• (Barnes Chu 2010): Barnes, David J.; Chu, Dominique (2010). Introduction to Modeling for Biosciences . Springer Verlag. Bibcode :2010itmf.book.....B.
• (Ramaswamy González-Segredo Sbalzarini 2009): Ramaswamy, Rajesh; González-Segredo, Nélido; Sbalzarini, Ivo F. (2009). 「化学反応ネットワークのための新しいクラスの非常に効率的な正確な確率的シミュレーションアルゴリズム」。Journal of Chemical Physics . 130 (24): 244104. arXiv : 0906.1992 . Bibcode :2009JChPh.130x4104R. doi :10.1063/1.3154624. PMID  19566139. S2CID  4952205.
• (Ramaswamy Sbalzarini 2010): Ramaswamy, Rajesh; Sbalzarini, Ivo F. (2010). 「化学反応ネットワークの構成拒否確率シミュレーションアルゴリズムの部分的傾向バリアント」(PDF) . Journal of Chemical Physics . 132 (4): 044102. Bibcode :2010JChPh.132d4102R. doi :10.1063/1.3297948. PMID  20113014.
• (Indurkhya Beal 2010): Indurkhya, Sagar; Beal, Jacob S. (2005). Isalan, Mark (編). 「反応因数分解と二部更新グラフが大規模生化学システムのための Gillespie アルゴリズムを加速」. PLOS ONE . 5 (1): e8125. Bibcode :2010PLoSO...5.8125I. doi : 10.1371/journal.pone.0008125 . PMC 2798956 . PMID  20066048. 
• (Ramaswamy Sbalzarini 2011): Ramaswamy, Rajesh; Sbalzarini, Ivo F. (2011). 「遅延を伴う化学反応ネットワークの確率的シミュレーションアルゴリズムの部分的傾向定式化」(PDF) . Journal of Chemical Physics . 134 (1): 014106. Bibcode :2011JChPh.134a4106R. doi :10.1063/1.3521496. PMID  21218996. S2CID  4949530.
• (Yates Klingbeil 2013): Yates, Christian A.; Klingbeil, Guido (2013). 「確率的シミュレーションアルゴリズムにおける乱数のリサイクル」. Annual Review of Physical Chemistry . 58 (9): 094103. Bibcode :2013JChPh.138i4103Y. doi :10.1063/1.4792207. PMID  23485273.
• Gillespie, Daniel T. (2007). 「化学反応速度論の確率的シミュレーション」. Annual Review of Physical Chemistry . 58 : 35–55. Bibcode :2007ARPC...58...35G. doi :10.1146/annurev.physchem.58.032806.104637. PMID  17037977.