Algoritmo de Gillespie

En teoría de la probabilidad , el algoritmo de Gillespie (o algoritmo de Doob-Gillespie o algoritmo de simulación estocástica , SSA ) genera una trayectoria estadísticamente correcta (solución posible) de un sistema de ecuaciones estocásticas para el cual se conocen las velocidades de reacción . Fue creado por Joseph L. Doob y otros (alrededor de 1945), presentado por Daniel Gillespie en 1976 y popularizado en 1977 en un artículo donde lo utiliza para simular sistemas de reacciones químicas o bioquímicas de manera eficiente y precisa utilizando una potencia computacional limitada (véase simulación estocástica ). ^{[ 1 ]} A medida que las computadoras se han vuelto más rápidas, el algoritmo se ha utilizado para simular sistemas cada vez más complejos. El algoritmo es particularmente útil para simular reacciones dentro de las células, donde el número de reactivos es bajo y el seguimiento de cada reacción individual es computacionalmente factible. Matemáticamente, es una variante de un método de Monte Carlo dinámico y similar a los métodos de Monte Carlo cinético . Se utiliza ampliamente en biología de sistemas computacional .

Historia

El proceso que condujo al algoritmo reconoce varios pasos importantes. En 1931, Andrei Kolmogorov introdujo las ecuaciones diferenciales correspondientes a la evolución temporal de procesos estocásticos que proceden mediante saltos, conocidas hoy como ecuaciones de Kolmogorov (proceso de salto de Markov) (una versión simplificada se conoce como ecuación maestra en las ciencias naturales). Fue William Feller , en 1940, quien halló las condiciones bajo las cuales las ecuaciones de Kolmogorov admitían probabilidades (propias) como soluciones. En su Teorema I (trabajo de 1940) establece que el tiempo hasta el siguiente salto se distribuye exponencialmente y que la probabilidad del siguiente evento es proporcional a la tasa. De este modo, estableció la relación de las ecuaciones de Kolmogorov con los procesos estocásticos . Posteriormente, Doob (1942, 1945) extendió las soluciones de Feller más allá del caso de procesos de salto puro. El método fue implementado en computadoras por David George Kendall (1950) utilizando la computadora Manchester Mark 1 y posteriormente empleado por Maurice S. Bartlett (1953) en sus estudios sobre brotes epidémicos. Gillespie (1977) obtiene el algoritmo de una manera diferente mediante un argumento físico.

Idea

Matemáticas

En una cámara de reacción, existe un número finito de moléculas. En cada instante infinitesimal de tiempo, puede tener lugar una única reacción. La velocidad de reacción está determinada por el número de moléculas de cada especie química.

Ingenuamente, podemos simular la trayectoria de la cámara de reacción discretizando el tiempo y simulando cada paso temporal. Sin embargo, pueden existir largos periodos de tiempo en los que no se produzca ninguna reacción. El algoritmo de Gillespie toma una muestra aleatoria de un tiempo de espera hasta que se produzca alguna reacción, y luego toma otra muestra aleatoria para determinar qué reacción se ha producido.

Los supuestos clave son que

Cada reacción es markoviana en el tiempo.
no hay correlación entre las reacciones

Dadas estas dos suposiciones, el tiempo de espera aleatorio para una reacción determinada sigue una distribución exponencial, donde la tasa exponencial es la suma de las tasas de las reacciones individuales.

Validez en simulaciones bioquímicas

Las ecuaciones de velocidad bioquímicas tradicionales, continuas y deterministas, no predicen con precisión las reacciones celulares, ya que se basan en reacciones a gran escala que requieren la interacción de millones de moléculas. Suelen modelarse como un conjunto de ecuaciones diferenciales ordinarias acopladas. En cambio, el algoritmo de Gillespie permite una simulación discreta y estocástica de un sistema con pocos reactivos, puesto que cada reacción se simula explícitamente. Una trayectoria correspondiente a una única simulación de Gillespie representa una muestra exacta de la función de probabilidad que es la solución de la ecuación maestra .

La base física del algoritmo reside en la colisión de moléculas dentro de un recipiente de reacción. Se asume que las colisiones son frecuentes, pero que aquellas con la orientación y energía adecuadas son poco frecuentes. Se asume que el entorno de reacción está bien mezclado.

Algoritmo

Una revisión ^{[ 1 ]} describe tres formulaciones diferentes, pero equivalentes; los métodos directo, de primera reacción y de primera familia, donde los dos primeros son casos especiales del último. La formulación de los métodos directo y de primera reacción se centra en realizar los pasos habituales de inversión de Monte Carlo sobre la llamada "premisa fundamental de la cinética química estocástica ", que matemáticamente es la función

p(\tau ,j\mid {\boldsymbol {x}},t)=a_{j}({\boldsymbol {x}})\exp \left(-\tau \sum _{j}a_{j}({\boldsymbol {x}})\right),

donde cada uno de los $a$ Los términos son funciones de propensión de una reacción elemental, cuyo argumento es ${\boldsymbol {x}}$ , el vector de recuentos de especies. El $\tau$ El parámetro es el tiempo hasta la siguiente reacción (o tiempo de permanencia), y $t$ es la hora actual. Parafraseando a Gillespie, esta expresión se lee como "la probabilidad, dado ${\boldsymbol {X}}(t)={\boldsymbol {x}}$ que la siguiente reacción del sistema ocurrirá en el intervalo de tiempo infinitesimal. $[t+\tau ,t+\tau +d\tau ]$ y tendrá una estequiometría correspondiente a la $j$ la reacción". Esta formulación proporciona una ventana a los métodos directos y de primera reacción al implicar $\tau$ es una variable aleatoria con distribución exponencial , y $j$ es "una variable aleatoria entera estadísticamente independiente con probabilidades puntuales $a_{j}({\boldsymbol {x}})/\sum _{j}a_{j}({\boldsymbol {x}})$ ".

Por lo tanto, el método de generación de Monte Carlo consiste simplemente en extraer dos números pseudoaleatorios, $r_{1}$ y $r_{2}$ en $[0,1]$ y calcular

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

y

j={}

el entero más pequeño que satisface

\sum _{j'=1}^{j}a_{j'}({\boldsymbol {x}})>r_{2}\sum _{j}a_{j}({\boldsymbol {x}}).

Utilizando este método de generación para el tiempo de permanencia y la siguiente reacción, Gillespie enuncia el algoritmo del método directo como:

1. Inicializar la hora $t=t_{0}$ y el estado del sistema ${\boldsymbol {x}}={\boldsymbol {x}}_{0}$  2. Con el sistema en estado ${\boldsymbol {x}}$ en ese momento $t$ , evaluar todos los $a_{j}({\boldsymbol {x}})$ y su suma $\sum _{j}a_{j}({\boldsymbol {x}})$  3. Calcula el valor anterior de $\tau$ y $j$  4. Realizar la siguiente reacción reemplazando $t\leftarrow t+\tau$ y ${\boldsymbol {x}}\leftarrow {\boldsymbol {x}}+\nu _{j}$  5. Registro $({\boldsymbol {x}},t)$ Como se desee. Vuelva al paso 2 o finalice la simulación.

dónde $\nu _{j}$ representa agregar el $j^{\text{th}}$ componente del vector de cambio de estado dado $\nu$ Esta familia de algoritmos es computacionalmente costosa, por lo que existen numerosas modificaciones y adaptaciones, incluyendo el método de reacción siguiente (Gibson y Bruck), el salto tau , así como técnicas híbridas donde los reactivos abundantes se modelan con un comportamiento determinista. Las técnicas adaptadas generalmente comprometen la exactitud de la teoría subyacente al algoritmo en su conexión con la ecuación maestra, pero ofrecen implementaciones razonables para escalas de tiempo mucho mejores. El costo computacional de las versiones exactas del algoritmo está determinado por la clase de acoplamiento de la red de reacciones. En redes débilmente acopladas, el número de reacciones influenciadas por cualquier otra reacción está limitado por una pequeña constante. En redes fuertemente acopladas, la activación de una sola reacción puede, en principio, afectar a todas las demás reacciones. Se ha desarrollado una versión exacta del algoritmo con escala de tiempo constante para redes débilmente acopladas, lo que permite la simulación eficiente de sistemas con un número muy grande de canales de reacción (Slepoy Thompson Plimpton 2008). El algoritmo generalizado de Gillespie, que considera las propiedades no markovianas de los eventos bioquímicos aleatorios con retardo, fue desarrollado por Bratsun et al. (2005) e independientemente por Barrio et al. (2006) y Cai (2007). Consulte los artículos citados a continuación para obtener más detalles.

Las formulaciones de propensión parcial, desarrolladas independientemente por Ramaswamy et al. (2009, 2010) e Indurkhya y Beal (2010), permiten construir una familia de versiones exactas del algoritmo cuyo coste computacional es proporcional al número de especies químicas en la red, en lugar del (mayor) número de reacciones. Estas formulaciones pueden reducir el coste computacional a una escala de tiempo constante para redes débilmente acopladas y a una escala lineal, como máximo, con el número de especies para redes fuertemente acopladas. También se ha propuesto una variante de propensión parcial del algoritmo de Gillespie generalizado para reacciones con retardos (Ramaswamy Sbalzarini 2011). El uso de métodos de propensión parcial se limita a reacciones químicas elementales, es decir, reacciones con como máximo dos reactivos diferentes. Toda reacción química no elemental puede descomponerse de forma equivalente en un conjunto de reacciones elementales, a costa de un aumento lineal (del orden de la reacción) en el tamaño de la red.

Ejemplos

Unión reversible de A y B para formar dímeros AB

Un ejemplo sencillo puede ayudar a explicar cómo funciona el algoritmo de Gillespie. Consideremos un sistema de moléculas de dos tipos, $A$ y $B.$ En este sistema, $A$ y $B$ se unen reversiblemente para formar dímeros $AB,$ de modo que son posibles dos reacciones: o bien A y B reaccionan reversiblemente para formar un dímero $AB$ , o bien un dímero $AB$ se disocia en $A$ y $B.$ La constante de velocidad de reacción para una molécula A dada que reacciona con una molécula $B dada es$ $k_{\mathrm {D} }$ y la velocidad de reacción para la ruptura de un dímero $AB es$ ${\ Displaystyle k _ {\ mathrm {B}}}$ .

Si en el instante t hay una molécula de cada tipo, entonces la tasa de formación de dímeros es $k_{\mathrm {D} }$ , mientras que si hay $n_{\mathrm {A} }$ moléculas de tipo $A$ y ${\ Displaystyle n _ {\ mathrm {B}}}$ moléculas de tipo $B$ , la tasa de formación de dímeros es $k_{\mathrm {D} }n_{\mathrm {A} }n_{\mathrm {B} }$ . Si hay ${\ Displaystyle n _ {\ mathrm {AB}}}$ dímeros entonces la tasa de disociación del dímero es $k_{\mathrm {B} }n_{\mathrm {AB} }$ .

La velocidad de reacción total, $R_{\mathrm {TOT} }$ , en el instante t viene dado por

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

Así pues, hemos descrito un modelo sencillo con dos reacciones. Esta definición es independiente del algoritmo de Gillespie. A continuación, describiremos cómo aplicar el algoritmo de Gillespie a este sistema.

En el algoritmo, avanzamos en el tiempo en dos pasos: calculando el tiempo hasta la siguiente reacción y determinando cuál de las posibles reacciones es la siguiente. Se supone que las reacciones son completamente aleatorias, por lo que si la tasa de reacción en un tiempo t es $R_{\mathrm {TOT} }$ , entonces el tiempo, δ t , hasta que ocurra la siguiente reacción es un número aleatorio extraído de una función de distribución exponencial con media $1/R_{\mathrm {TOT} }$ . Por lo tanto, avanzamos el tiempo de t a t + δ t .

La probabilidad de que esta reacción sea una molécula $A$ que se une a una molécula $B$ es simplemente la fracción de la velocidad total debida a este tipo de reacción, es decir,

la probabilidad de que la reacción sea $P({\ce {{A}+ B -> AB}})=k_{D}n_{A}n_{B}/R_{{\ce {TOT}}}$

La probabilidad de que la siguiente reacción sea la disociación de un dímero $AB$ es simplemente 1 menos eso. Entonces, con estas dos probabilidades, formamos un dímero mediante reducción $n_{\mathrm {A} }$ y ${\ Displaystyle n _ {\ mathrm {B}}}$ en uno y aumentar ${\ Displaystyle n _ {\ mathrm {AB}}}$ por uno, o disociamos un dímero y aumentamos $n_{\mathrm {A} }$ y ${\ Displaystyle n _ {\ mathrm {B}}}$ por uno y disminuir ${\ Displaystyle n _ {\ mathrm {AB}}}$ por uno.

Ahora ambos hemos avanzado el tiempo a t + δ t y hemos realizado una sola reacción. El algoritmo de Gillespie simplemente repite estos dos pasos tantas veces como sea necesario para simular el sistema durante el tiempo que queramos (es decir, para tantas reacciones). El resultado de una simulación de Gillespie que comienza con $n_{\mathrm {A} }=n_{\mathrm {B} }=10$ y ${\ Displaystyle n _ {\ mathrm {AB}} = 0}$ en t =0, y donde $k_{\mathrm {D} }=2$ y $k_{\mathrm {B} }=1$ , se muestra a la derecha. Para estos valores de parámetros, en promedio hay 8 ${\ Displaystyle n _ {\ mathrm {AB}}}$ dímeros y 2 de $A$ y $B$ , pero debido a la pequeña cantidad de moléculas, las fluctuaciones alrededor de estos valores son grandes. El algoritmo de Gillespie se usa a menudo para estudiar sistemas donde estas fluctuaciones son importantes.

Este fue solo un ejemplo sencillo con dos reacciones. Los sistemas más complejos con más reacciones se manejan de la misma manera. Todas las velocidades de reacción deben calcularse en cada paso de tiempo, y se elige una con una probabilidad igual a su contribución fraccional a la velocidad. Luego, el tiempo avanza como en este ejemplo.

Referencias

1 2 Gillespie, Daniel T. (2007-05-01). "Simulación estocástica de la cinética química" . Annual Review of Physical Chemistry . 58 (1): 35– 55. Bibcode : 2007ARPC...58...35G . doi : 10.1146/annurev.physchem.58.032806.104637 . ISSN 0066-426X . PMID 17037977 .

Lecturas adicionales

Gillespie, Daniel T. (1977). "Simulación estocástica exacta de reacciones químicas acopladas". 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). "Un método general para simular numéricamente la evolución temporal estocástica de reacciones químicas acopladas". 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). "Simulación estocástica exacta y eficiente de sistemas químicos con muchas especies y muchos canales" (PDF) . Journal of Physical Chemistry A. 104 ( 9): 1876– 1889. Bibcode : 2000JPCA..104.1876G . doi : 10.1021/jp993732q .
Doob, Jacob L. (1942). "Temas en la teoría de las cadenas de Markoff" . Transactions of the American Mathematical Society . 52 (1): 37– 64. doi : 10.1090/S0002-9947-1942-0006633-7 . JSTOR 1990152 .
Doob, Jacob L. (1945). "Cadenas de Markoff: caso numerable". Transactions of the American Mathematical Society . 58 (3): 455– 473. doi : 10.2307/1990339 . JSTOR 1990339 .
Press, William H.; Teukolsky, Saul A.; Vetterling, William T.; Flannery, Brian P. (2007). «Sección 17.7. Simulación estocástica de redes de reacciones químicas» . Numerical Recipes: The Art of Scientific Computing (3.ª ed.). Nueva York, NY: Cambridge University Press. ISBN 978-0-521-88068-8Archivado del original el 11 de agosto de 2011. Consultado el 17 de agosto de 2011 .
Kolmogorov, Andrey N. (1931). "Über die analytischen Methoden in der Wahrscheinlichkeitsrechnung" [ Sobre los métodos analíticos en la teoría de la probabilidad ] . Annalen Matemáticas . 104 : 415– 458. doi : 10.1007/BF01457949 . S2CID 119439925 .
Feller, Willy (1940). "Sobre las ecuaciones integro-diferenciales de procesos de Markoff puramente discontinuos" . Transactions of the American Mathematical Society . 48 (3): 4885–15 . doi : 10.2307/1990095 . JSTOR 1970064 .
Kendall, David G. (1950). "Una realización artificial de un proceso simple de "nacimiento y muerte"". Journal of the Royal Statistical Society, Serie B. 12 ( 1): 116– 119. doi : 10.1111/j.2517-6161.1950.tb00048.x . JSTOR 2983837 .
Bartlett, Maurice S. (1953). "Procesos estocásticos o la estadística del cambio". Journal of the Royal Statistical Society, Serie C. 2 ( 1): 44– 64. doi : 10.2307/2985327 . JSTOR 2985327 .
Rathinam, Muruhan; Petzold, Linda R.; Cao, Yang; Gillespie, Daniel T. (2003). "Rigidez en sistemas estocásticos de reacción química: El método implícito de salto tau". Journal of Chemical Physics . 119 (24): 12784– 12794. Bibcode : 2003JChPh.11912784R . doi : 10.1063/1.1627296 .
Sinitsyn, Nikolai A.; Hengartner, Nicolas; Nemenman, Ilya (2009). "Adiabatic coarse-graining and simulations of stochastic biochemical networks" . Proceedings of the National Academy of Sciences of the United States of America . 106 (20): 10546– 10551. Bibcode : 2009PNAS..10610546S . doi : 10.1073/pnas.0809340106 . PMC 2705573. PMID 19525397 .
Salis, Howard; Kaznessis, Yiannis N. (2005). "Simulación estocástica híbrida precisa de un sistema de reacciones químicas o bioquímicas acopladas". 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). "Un algoritmo cinético de Monte Carlo de tiempo constante para la simulación de grandes redes de reacciones bioquímicas". 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). "Oscilaciones estocásticas inducidas por retraso en la regulación genética" . Actas de la Academia Nacional de Ciencias de los Estados Unidos de América . 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). "Regulación oscilatoria de hes1 : modelado y simulación de retardo estocástico discreto" . 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). "Simulación estocástica exacta de reacciones químicas acopladas con retardos". 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). Introducción al modelado para biociencias . Springer Verlag. Bibcode : 2010itmf.book.....B .
(Ramaswamy González-Segredo Sbalzarini 2009): Ramaswamy, Rajesh; González-Segredo, Nélido; Sbalzarini, Ivo F. (2009). "Una nueva clase de algoritmos de simulación estocástica exactos altamente eficientes para redes de reacciones químicas". 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). "Una variante de propensión parcial del algoritmo de simulación estocástica de rechazo de composición para redes de reacciones químicas" (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 (ed.). " Reaction Factoring and Bipartite Update Graphs Accelerate the Gillespie Algorithm for Large-Scale Biochemical Systems" . 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). "Una formulación de propensión parcial del algoritmo de simulación estocástica para redes de reacciones químicas con retrasos" ( 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). "Reciclaje de números aleatorios en el algoritmo de simulación estocástica" . Annual Review of Physical Chemistry . 58 (9): 094103. Bibcode : 2013JChPh.138i4103Y . doi : 10.1063/1.4792207 . PMID 23485273 .
Gillespie, Daniel T. (2007). "Simulación estocástica de la cinética química". Annual Review of Physical Chemistry . 58 : 35–55 . Bibcode : 2007ARPC...58...35G . doi : 10.1146/annurev.physchem.58.032806.104637 . PMID 17037977 .

[Gillepie2007-1] 1 2 Gillespie, Daniel T. (2007-05-01). "Simulación estocástica de la cinética química" . Annual Review of Physical Chemistry . 58 (1): 35– 55. Bibcode : 2007ARPC...58...35G . doi : 10.1146/annurev.physchem.58.032806.104637 . ISSN 0066-426X . PMID 17037977 .

[ 1 ]