Algoritmo de Gillespie

Método para sistemas de ecuaciones estocásticas.

En teoría de la probabilidad , el algoritmo de Gillespie (o el algoritmo de Doob-Gillespie o algoritmo de simulación estocástica , el SSA ) genera una trayectoria estadísticamente correcta (posible solución) 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 Dan Gillespie en 1976 y popularizado en 1977 en un artículo donde lo utiliza para simular sistemas de reacciones químicos o bioquímicos de manera eficiente y precisa utilizando un poder computacional limitado (ver 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 células, donde la cantidad de reactivos es baja y realizar un seguimiento de cada reacción es computacionalmente factible. Matemáticamente, es una variante del método dinámico de Monte Carlo y similar a los métodos cinéticos de Monte Carlo . Se utiliza mucho en biología de sistemas computacionales . ^{[ cita necesaria ]}

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 se desarrollan mediante saltos, hoy conocidas 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 encontró las condiciones bajo las cuales las ecuaciones de Kolmogorov admitían probabilidades (propias) como soluciones. En su Teorema I (obra de 1940) establece que el tiempo hasta el siguiente salto se distribuyó exponencialmente y la probabilidad del siguiente evento es proporcional a la velocidad. Como tal, estableció la relación de las ecuaciones de Kolmogorov con los procesos estocásticos . Posteriormente, Doob (1942, 1945) amplió las soluciones de Feller más allá del caso de los procesos de salto puro. El método fue implementado en computadoras por David George Kendall (1950) usando la computadora Manchester Mark 1 y luego utilizado por Maurice S. Bartlett (1953) en sus estudios sobre brotes epidémicos. Gillespie (1977) obtiene el algoritmo de otra manera haciendo uso de un argumento físico.

Idea detrás del algoritmo

Las ecuaciones de velocidad bioquímicas continuas y deterministas tradicionales no predicen con precisión las reacciones celulares, ya que se basan en reacciones masivas que requieren las interacciones de millones de moléculas. Por lo general, se modelan como un conjunto de ecuaciones diferenciales ordinarias acopladas. Por el contrario, el algoritmo de Gillespie permite una simulación discreta y estocástica de un sistema con pocos reactivos porque 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 masa de probabilidad que es la solución de la ecuación maestra .

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

Algoritmo

Una revisión (Gillespie, 2007) describe tres formulaciones diferentes, pero equivalentes; los métodos directo, de primera reacción y de primera familia, siendo los dos primeros casos especiales del segundo. La formulación de los métodos directo y de primera reacción se centra en realizar los habituales pasos 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

${\displaystyle 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 términos son funciones de propensión de una reacción elemental, cuyo argumento es , el vector de especies cuenta. El parámetro es el tiempo hasta la siguiente reacción (o tiempo de estancia) y es el tiempo actual. Parafraseando a Gillespie, esta expresión se lee como "la probabilidad, dada , de que la siguiente reacción del sistema ocurra en el intervalo de tiempo infinitesimal , y será de estequiometría correspondiente a la reacción enésima". Esta formulación proporciona una ventana a los métodos directo y de primera reacción al implicar que es una variable aleatoria distribuida exponencialmente y es "una variable aleatoria entera estadísticamente independiente con probabilidades puntuales ". ${\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}})}$

Por lo tanto, el método de generación de Monte Carlo consiste simplemente en dibujar dos números pseudoaleatorios, y en , y calcular ${\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),}$

y

${\displaystyle j={}}$el entero más pequeño que satisface ${\displaystyle \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 estancia y la siguiente reacción, Gillespie establece el algoritmo del método directo como

1. Inicialice el tiempo y el estado del sistema
2. Con el sistema en el estado en el momento , evalúe todos y su suma
3. Calcule el valor anterior de y
4. Efectúe la siguiente reacción reemplazando y
5. Registre como desee. Regrese al paso 2 o finalice la simulación. ${\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)}$

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

Las formulaciones de propensión parcial, desarrolladas de forma independiente por Ramaswamy et al. (2009, 2010) e Indurkhya y Beal (2010), están disponibles para construir una familia de versiones exactas del algoritmo cuyo costo 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 costo computacional a un escalamiento de tiempo constante para redes débilmente acopladas y a una escala como máximo lineal con el número de especies para redes fuertemente acopladas. También se ha propuesto una variante de propensión parcial del algoritmo generalizado de Gillespie para reacciones con retrasos (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. Cada reacción química no elemental se puede descomponer de manera equivalente en un conjunto de reacciones elementales, a expensas de un aumento lineal (en el 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. Considere 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: A y B reaccionan reversiblemente para formar un dímero AB , o un dímero AB se disocia en A y B. La constante de velocidad de reacción para una determinada molécula A que reacciona con una determinada molécula B es , y la velocidad de reacción para la ruptura de un dímero AB es . ${\displaystyle k_{\mathrm {D} }}$ ${\displaystyle k_{\mathrm {B} }}$

Si en el momento t hay una molécula de cada tipo entonces la velocidad de formación de dímeros es , mientras que si hay moléculas de tipo A y moléculas de tipo B , la velocidad de formación de dímeros es . Si hay dímeros, entonces la tasa de disociación de los dímeros es . ${\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} }}$

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

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

Entonces, ahora hemos descrito un modelo simple con dos reacciones. Esta definición es independiente del algoritmo de Gillespie. Ahora 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 velocidad de reacción en un momento t es , entonces el tiempo, δ t , hasta que ocurra la siguiente reacción es un número aleatorio extraído de la función de distribución exponencial con media . Por tanto, adelantamos el tiempo de t a t + δ t . ${\displaystyle R_{\mathrm {TOT} }}$ ${\displaystyle 1/R_{\mathrm {TOT} }}$

Gráfico del número de moléculas A (curva negra) y dímeros AB en función del tiempo. Como comenzamos con 10 moléculas A y B en el instante t = 0, el número de moléculas B siempre es igual al número de moléculas A , por lo que no se muestra.

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 ${\displaystyle 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 solo 1 menos eso. Entonces, con estas dos probabilidades formamos un dímero reduciendo y en uno y aumentando en uno, o disociamos un dímero y aumentando y en uno y disminuyendo en uno. ${\displaystyle n_{\mathrm {A} }}$ ${\displaystyle n_{\mathrm {B} }}$ ${\displaystyle n_{\mathrm {AB} }}$ ${\displaystyle n_{\mathrm {A} }}$ ${\displaystyle n_{\mathrm {B} }}$ ${\displaystyle n_{\mathrm {AB} }}$

Ahora hemos avanzado el tiempo hasta t + δ t y hemos realizado una única 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, durante tantas reacciones). El resultado de una simulación de Gillespie que comienza con y en t =0, y donde y , se muestra a la derecha. Para estos valores de parámetros, en promedio hay 8 dímeros y 2 de A y B , pero debido al pequeño número de moléculas, las fluctuaciones alrededor de estos valores son grandes. El algoritmo de Gillespie se utiliza a menudo para estudiar sistemas en los que estas fluctuaciones son importantes. ${\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} }}$

Ese fue sólo un ejemplo simple, 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 elegirse una con probabilidad igual a su contribución fraccionaria a la velocidad. Luego se avanza el tiempo como en este ejemplo.

Autoensamblaje estocástico

El modelo de Gard describe el autoensamblaje de lípidos en agregados. Utilizando simulaciones estocásticas se muestra el surgimiento de múltiples tipos de agregados y su evolución.

Referencias

1. Gillespie, Daniel T. (1 de mayo de 2007). "Simulación estocástica de la cinética química". Revista Anual de Química Física . 58 (1): 35–55. Código Bib : 2007ARPC...58...35G. doi : 10.1146/annurev.physchem.58.032806.104637. ISSN  0066-426X. PMID  17037977.

Otras lecturas

• Gillespie, Daniel T. (1977). "Simulación estocástica exacta de las reacciones químicas acopladas". El diario de la química física . 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". Revista de Física Computacional . 22 (4): 403–434. Código bibliográfico : 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) . Revista de Química Física A. 104 (9): 1876–1889. Código Bib : 2000JPCA..104.1876G. doi :10.1021/jp993732q.
• Doob, Jacob L. (1942). "Temas de la teoría de las cadenas de Markoff". Transacciones de la Sociedad Matemática Estadounidense . 52 (1): 37–64. doi : 10.1090/S0002-9947-1942-0006633-7 . JSTOR  1990152.
• Doob, Jacob L. (1945). "Cadenas de Markoff - Caso numerable". Transacciones de la Sociedad Matemática Estadounidense . 58 (3): 455–473. doi :10.2307/1990339. JSTOR  1990339.
• Prensa, William H.; Teukolsky, Saúl A.; Vetterling, William T.; Flannery, Brian P. (2007). "Sección 17.7. Simulación estocástica de redes de reacciones químicas". Recetas numéricas: el arte de la informática científica (3ª ed.). Nueva York, Nueva York: Cambridge University Press. ISBN 978-0-521-88068-8. Archivado desde el 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 calificación puramente discontinuos". Transacciones de la Sociedad Matemática Estadounidense . 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 ". Revista de la Royal Statistical Society, Serie B. 12 (1): 116-119. JSTOR  2983837.
• Bartlett, Maurice S. (1953). "Procesos estocásticos o estadísticas de cambio". Revista de la 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 que reaccionan químicamente: el método implícito del salto tau". Revista de Física Química . 119 (24): 12784–12794. Código Bib : 2003JChPh.11912784R. doi :10.1063/1.1627296.
• Sinitsyn, Nikolai A.; Hengartner, Nicolás; Nemenman, Ilya (2009). "Grano grueso adiabático y simulaciones de redes bioquímicas estocásticas". Actas de la Academia Nacional de Ciencias de los Estados Unidos de América . 106 (20): 10546–10551. Código bibliográfico : 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". Revista de Física Química . 122 (5): 054103. Código bibliográfico : 2005JChPh.122e4103S. doi :10.1063/1.1835951. PMID  15740306.
• (Slepoy Thompson Plimpton 2008): Slepoy, Alexander; Thompson, Aidan P.; Plimpton, Steven J. (2008). "Un algoritmo de Monte Carlo cinético de tiempo constante para la simulación de grandes redes de reacciones bioquímicas". Revista de Física Química . 128 (20): 205101. Código bibliográfico : 2008JChPh.128t5101S. doi : 10.1063/1.2919546. PMID  18513044.
• (Bratsun et al. 2005): Bratsun, Dmitri; Volfson, Dmitri; Apresurado, 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. Código Bib : 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 Biología Computacional . 2 (9): 1017. Código Bib : 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 retrasos". Revista de Física Química . 126 (12): 124108. Código bibliográfico : 2007JChPh.126l4108C. doi :10.1063/1.2710253. PMID  17411109.
• (Barnes Chu 2010): Barnes, David J.; Chu, Dominique (2010). Introducción a la Modelización para Biociencias . Springer Verlag. Código Bib : 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 exacta altamente eficientes para redes de reacciones químicas". Revista de Física Química . 130 (24): 244104. arXiv : 0906.1992 . Código Bib :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) . Revista de Física Química . 132 (4): 044102. Código bibliográfico : 2010JChPh.132d4102R. doi : 10.1063/1.3297948. PMID  20113014.
• (Indurkhya Beal 2010): Indurkhya, Sagar; Beal, Jacob S. (2005). Isalan, Mark (ed.). "El factoraje de reacción y los gráficos de actualización bipartita aceleran el algoritmo de Gillespie para sistemas bioquímicos a gran escala". MÁS UNO . 5 (1): e8125. Código Bib : 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) . Revista de Física Química . 134 (1): 014106. Código bibliográfico : 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". Revista Anual de Química Física . 58 (9): 094103. Código bibliográfico : 2013JChPh.138i4103Y. doi : 10.1063/1.4792207. PMID  23485273.
• Gillespie, Daniel T. (2007). "Simulación estocástica de la cinética química". Revista Anual de Química Física . 58 : 35–55. Código Bib : 2007ARPC...58...35G. doi : 10.1146/annurev.physchem.58.032806.104637. PMID  17037977.