Algoritmo de Gillespie - Gillespie algorithm

Na teoria da probabilidade , o algoritmo de Gillespie (ou ocasionalmente o algoritmo Doob-Gillespie ) gera uma trajetória estatisticamente correta (solução possível) de um sistema de equações estocásticas para as quais as taxas de reação são conhecidas. Foi criado por Joseph L. Doob e outros (por volta de 1945), apresentado por Dan Gillespie em 1976, e popularizado em 1977 em um artigo onde ele o usa para simular sistemas químicos ou bioquímicos de reações com eficiência e precisão usando poder computacional limitado (ver simulação estocástica ). À medida que os computadores se tornaram mais rápidos, o algoritmo foi usado para simular sistemas cada vez mais complexos. O algoritmo é particularmente útil para simular reações dentro das células, onde o número de reagentes é baixo e manter o controle da posição e do comportamento de moléculas individuais é computacionalmente viável. Matematicamente, é uma variante de um método de Monte Carlo dinâmico e semelhante aos métodos de Monte Carlo cinéticos . É muito usado em biologia de sistemas computacionais .

História

O processo que levou ao algoritmo reconhece várias etapas importantes. Em 1931, Andrei Kolmogorov introduziu as equações diferenciais correspondentes à evolução no tempo dos processos estocásticos que procedem por saltos, hoje conhecidas como equações de Kolmogorov (processo de salto de Markov) (uma versão simplificada é conhecida como equação mestre nas ciências naturais). Foi William Feller , em 1940, quem encontrou as condições sob as quais as equações de Kolmogorov admitiam as probabilidades (adequadas) como soluções. Em seu Teorema I (trabalho de 1940), ele estabelece que o tempo para o próximo salto foi distribuído exponencialmente e a probabilidade do próximo evento é proporcional à taxa. Como tal, ele estabeleceu a relação das equações de Kolmogorov com os processos estocásticos . Mais tarde, Doob (1942,1945) estendeu as soluções de Feller para além do caso dos processos de salto puro. O método foi implementado em computadores por David George Kendall (1950) usando o computador Manchester Mark 1 e mais tarde usado por Maurice S. Bartlett (1953) em seus estudos de surtos de epidemias. Gillespie (1977) obtém o algoritmo de uma maneira diferente, fazendo uso de um argumento físico.

Ideia por trás do algoritmo

As equações de taxas bioquímicas contínuas e determinísticas tradicionais não prevêem com precisão as reações celulares, uma vez que dependem de reações em massa que requerem as interações de milhões de moléculas. Eles são normalmente modelados como um conjunto de equações diferenciais ordinárias acopladas. Em contraste, o algoritmo de Gillespie permite uma simulação discreta e estocástica de um sistema com poucos reagentes porque cada reação é explicitamente simulada. Uma trajetória correspondente a uma única simulação de Gillespie representa uma amostra exata da função de massa de probabilidade que é a solução da equação mestre .

A base física do algoritmo é a colisão de moléculas dentro de um vaso de reação. Presume-se que as colisões são frequentes, mas as colisões com a orientação e energia adequadas são raras. Portanto, todas as reações dentro da estrutura de Gillespie devem envolver no máximo duas moléculas. As reações envolvendo três moléculas são consideradas extremamente raras e são modeladas como uma sequência de reações binárias. Presume-se também que o ambiente de reação está bem misturado.

Algoritmo

Uma revisão recente (Gillespie, 2007) descreve três formulações diferentes, mas equivalentes; os métodos direto, de primeira reação e de primeira família, em que os dois primeiros são casos especiais do último. A formulação dos métodos de reação direta e de primeira reação é centrada na execução das etapas usuais de inversão de Monte-Carlo na chamada "premissa fundamental da cinética química estocástica", que matematicamente é a função

${\ displaystyle p (\ tau, j | {\ boldsymbol {x}}, t) = a_ {j} ({\ boldsymbol {x}}) \ exp (\ tau \ sum _ {j} a_ {j} ( {\ boldsymbol {x}}))}$,

onde cada um dos termos são funções de propensão de uma reação elementar, cujo argumento é , o vetor de contagem de espécies. O parâmetro é o tempo para a próxima reação (ou tempo de permanência) e é, obviamente, o tempo atual. Parafraseando Gillespie, essa expressão é lida como "a probabilidade, dada , de que a próxima reação do sistema ocorra no intervalo de tempo infinitesimal , e será de estequiometria correspondente à ª reação". Esta formulação fornece uma janela para os métodos direto e de primeira reação, implicando que é uma variável aleatória distribuída exponencialmente e é "uma variável aleatória inteira estatisticamente independente com probabilidades pontuais ". ${\ 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}})}$

Assim, o método de geração de Monte-Carlo é simplesmente desenhar dois números pseudo-aleatórios e , em seguida, 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} }}\direito)}$,

e

${\ displaystyle j =}$o menor inteiro satisfaz .${\ displaystyle \ sum _ {j '= 1} ^ {j} a_ {j'} ({\ boldsymbol {x}})> r_ {2} \ sum _ {j} a_ {j} ({\ boldsymbol { x}})}$

Utilizando este método de geração para o tempo de permanência e próxima reação, o algoritmo do método direto é declarado por Gillespie como

1. Initialize the time ${\displaystyle t=t_{0}}$ and the system's state ${\displaystyle {\boldsymbol {x}}={\boldsymbol {x}}_{0}}$
2. With the system in state ${\displaystyle {\boldsymbol {x}}}$ at time ${\displaystyle t}$, evaluate all the ${\displaystyle a_{j}({\boldsymbol {x}})}$ and their sum ${\displaystyle \sum _{j}a_{j}({\boldsymbol {x}})}$
3. Effect the next reaction by replacing ${\displaystyle t\leftarrow t+\tau }$ and ${\displaystyle {\boldsymbol {x}}\leftarrow {\boldsymbol {x}}+\nu _{j}}$
4. Record ${\displaystyle ({\boldsymbol {x}},t)}$ as desired. Return to step 1, or else end the simulation.

Esta família de algoritmos é computacionalmente cara e, portanto, existem muitas modificações e adaptações, incluindo o próximo método de reação (Gibson & Bruck), salto tau , bem como técnicas híbridas onde reagentes abundantes são modelados com comportamento determinístico. As técnicas adaptadas geralmente comprometem a exatidão da teoria por trás do algoritmo, uma vez que se conecta à equação mestre, mas oferecem realizações razoáveis ​​para escalas de tempo muito melhoradas. O custo computacional de versões exatas do algoritmo é determinado pela classe de acoplamento da rede de reação. Em redes fracamente acopladas, o número de reações que é influenciado por qualquer outra reação é limitado por uma pequena constante. Em redes fortemente acopladas, o disparo de uma única reação pode, em princípio, afetar todas as outras reações. Uma versão exata do algoritmo com escala de tempo constante para redes fracamente acopladas foi desenvolvida, permitindo a simulação eficiente de sistemas com um grande número de canais de reação (Slepoy Thompson Plimpton 2008). O algoritmo de Gillespie generalizado que leva em consideração as propriedades não markovianas de eventos bioquímicos aleatórios com atraso foi desenvolvido por Bratsun et al. 2005 e independentemente Barrio et al. 2006, bem como (Cai 2007). Veja os artigos citados abaixo para mais detalhes.

Formulações de propensão parcial, desenvolvidas independentemente por Ramaswamy et al. (2009, 2010) e Indurkhya e Beal (2010), estão disponíveis para construir uma família de versões exatas do algoritmo cujo custo computacional é proporcional ao número de espécies químicas na rede, ao invés do (maior) número de reações. Essas formulações podem reduzir o custo computacional para escalonamento de tempo constante para redes fracamente acopladas e escalar no máximo linearmente com o número de espécies para redes fortemente acopladas. Uma variante de propensão parcial do algoritmo de Gillespie generalizado para reações com atrasos também foi proposta (Ramaswamy Sbalzarini 2011). O uso de métodos de propensão parcial é limitado a reações químicas elementares, ou seja, reações com no máximo dois reagentes diferentes. Cada reação química não elementar pode ser decomposta de forma equivalente em um conjunto de reações químicas elementares, às custas de um aumento linear (na ordem da reação) no tamanho da rede.

Considere a seguinte implementação orientada a objetos dos métodos de reação direta e de primeira reação, que estão contidos em uma classe Python 3:

from math import log
from random import random

class SSA:
    """Container for SSAs"""

    def __init__(self, model):
        """Initialize container with model and pseudorandom number generator"""
        self.model = model
        self.random = random

    def direct(self):
        """Indefinite generator of direct-method trajectories"""
        while True:
            while not self.model.exit():

                # evaluate weights and partition
                weights = [
                    (rxn, sto, pro(self.model))
                    for (rxn, sto, pro) in self.model.reactions
                ]
                partition = sum(w[-1] for w in weights)

                # evaluate sojourn time (MC step 1)
                sojourn = log(1.0 / self.random()) / partition
                self.model["time"].append(self.model["time"][-1] + sojourn)

                # evaluate the reaction (MC step 2)
                partition = partition * self.random()
                while partition >= 0.0:
                    rxn, sto, pro = weights.pop(0)
                    partition -= pro
                for species, delta in sto.items():
                    self.model[species].append(self.model[species][-1] + delta)

                self.model.curate()
            yield self.model
            self.model.reset()

    def first_reaction(self):
        """Indefinite generator of 1st-reaction trajectories"""
        while True:
            while not self.model.exit():

                # evaluate next reaction times
                times = [
                    (
                        log(
                            1.0 / self.random()
                        ) / pro(self.model),
                        sto
                    )
                    for (rxn, sto, pro) in self.model.reactions
                ]
                times.sort()

                # evaluate reaction time
                self.model["time"].append(
                    self.model["time"][-1] + times[0][0]
                )

                # evaluate reaction
                for species, delta in times[0][1].items():
                    self.model[species].append(
                        self.model[species][-1] + delta
                    )

                self.model.curate()
            yield self.model
            self.model.reset()

Os membros directe first_reactionsão geradores indefinidos, o que significa que eles produzirão continuamente trajetórias, cada trajetória sendo uma simulação completa do sistema, em um loop até que um sinal quebre esse loop. Uma advertência importante para qualquer implementação desse algoritmo é que ele não fornece nenhuma lógica para distinguir entre eventos elementares possíveis e impossíveis. Portanto, para realmente utilizar qualquer uma das versões do algoritmo em um loop e produzir algum número de trajetórias para análise, essa classe requer que um modelo seja passado a ela na instanciação. Em adição às espécies de habitação populações e funções de propensão, um modelo tem três métodos que o algoritmo Gillespie faz chamadas para: curate, exit, e reset. A finalidade desses métodos é especificada a seguir e, essencialmente, governa quando o algoritmo é encerrado. O raciocínio por trás da introdução da classe de modelo é principalmente para evitar misturar as propriedades cinéticas do processo específico que está sendo simulado com a lógica interna do algoritmo de Gillespie.

O modelo a seguir é uma subclasse dicionário com membros públicos curate, exite resetque, respectivamente,

• determinar quais reações são e não são válidas no final de cada iteração de uma dada trajetória;
• retornar Truese não houver mais reações possíveis;
• retornar os hashings do modelo aos seus valores originais (ou seja, suas condições iniciais) no final de uma determinada trajetória.

class SSAModel(dict):
    """Container for SSA model"""

    def __init__(
        self, initial_conditions, propensities, stoichiometry
    ):
        """Initialize model"""
        super().__init__(**initial_conditions)
        self.reactions = list()
        self.excluded_reactions = list()
        for reaction,propensity in propensities.items():
            if propensity(self) == 0.0:
                self.excluded_reactions.append(
                    (
                        reaction,
                        stoichiometry[reaction],
                        propensity
                    )
                )
            else:
                self.reactions.append(
                    (
                        reaction,
                        stoichiometry[reaction],
                        propensity
                    )
                )

    def exit(self):
        """Return True to break out of trajectory"""

        # return True if no more reactions
        if len(self.reactions) == 0: return True

        # return False if there are more reactions
        else: return False

    def curate(self):
        """Validate and invalidate model reactions"""
        
        # evaluate possible reactions
        reactions = []
        while len(self.reactions) > 0:
            reaction = self.reactions.pop()
            if reaction[2](self) == 0:
                self.excluded_reactions.append(reaction)
            else:
                reactions.append(reaction)
        self.reactions = reactions

        # evaluate impossible reactions
        excluded_reactions = []
        while len(self.excluded_reactions) > 0:
            reaction = self.excluded_reactions.pop()
            if reaction[2](self) > 0:
                self.reactions.append(reaction)
            else:
                excluded_reactions.append(reaction)
        self.excluded_reactions = excluded_reactions

    def reset(self):
        """Clear the trajectory"""

        # reset species to initial conditions
        for key in self: del self[key][1:]

        # reset reactions per initial conditions
        self.curate()

Exemplos

Ligação reversível de A e B para formar dímeros AB

Um exemplo simples pode ajudar a explicar como o algoritmo de Gillespie funciona. Considere-se um sistema de moléculas de dois tipos, A e B . Neste sistema, uma ea B de forma reversível ligam-se conjuntamente para formar AB dímeros de tal modo que duas reacções são possíveis: ou A e B reagir reversivelmente para formar uma AB dímero, ou uma AB dissocia-se em dímeros A e B . A constante da taxa de reação para uma determinada molécula A reagindo com uma determinada molécula B é , e a taxa de reação para a quebra de um dímero AB é . ${\ displaystyle k _ {\ mathrm {D}}}$${\ displaystyle k _ {\ mathrm {B}}}$

Se no tempo t houver uma molécula de cada tipo, então a taxa de formação do dímero é , enquanto se houver moléculas do tipo A e moléculas do tipo B , a taxa de formação do dímero é . Se houver dímeros, a taxa de dissociação do dímero é . ${\ 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}}}$

A taxa de reação total,, no tempo t é então 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}}}$

Portanto, agora descrevemos um modelo simples com duas reações. Esta definição é independente do algoritmo de Gillespie. Vamos agora descrever como aplicar o algoritmo de Gillespie a este sistema.

No algoritmo, avançamos no tempo em duas etapas: calculando o tempo para a próxima reação e determinando qual das possíveis reações é a próxima reação. As reações são assumidas como completamente aleatórias, então se a taxa de reação em um tempo t for , então o tempo, δ t , até que a próxima reação ocorra é um número aleatório retirado da função de distribuição exponencial com média . Assim, avançamos no tempo de t para t + δ t . ${\ displaystyle R _ {\ mathrm {TOT}}}$${\ displaystyle 1 / R _ {\ mathrm {TOT}}}$

Gráfico do número de moléculas A (curva preta) e dímeros AB em função do tempo. Como começamos com 10 moléculas A e B no tempo t = 0, o número de moléculas B é sempre igual ao número de moléculas A e, portanto, não é mostrado.

A probabilidade de que essa reação seja uma molécula A se ligando a uma molécula B é simplesmente a fração da taxa total devido a esse tipo de reação, ou seja,

a probabilidade de que a reação seja ${\ displaystyle P ({\ ce {{A} + B -> AB}}) = k_ {D} n_ {A} n_ {B} / R _ {{\ ce {TOT}}}}$

A probabilidade de que a próxima reação seja a dissociação de um dímero AB é de apenas 1 menos isso. Assim, com essas duas probabilidades, ou formamos um dímero reduzindo e em um e aumentamos em um, ou dissociamos um dímero e aumentamos e em um e diminuímos em um. ${\ displaystyle n _ {\ mathrm {A}}}$${\ displaystyle n _ {\ mathrm {B}}}$${\ displaystyle n _ {\ mathrm {AB}}}$${\ displaystyle n _ {\ mathrm {A}}}$${\ displaystyle n _ {\ mathrm {B}}}$${\ displaystyle n _ {\ mathrm {AB}}}$

Agora avançamos o tempo para t + δ t e realizamos uma única reação. O algoritmo de Gillespie apenas repete essas duas etapas quantas vezes forem necessárias para simular o sistema pelo tempo que quisermos (ou seja, para quantas reações). O resultado de uma simulação de Gillespie que começa com e em t = 0, e onde e , é mostrado à direita. Para esses valores de parâmetro, em média existem 8 dímeros e 2 de A e B, mas devido ao pequeno número de moléculas, as flutuações em torno desses valores são grandes. O algoritmo de Gillespie é freqüentemente usado para estudar sistemas onde essas flutuações são 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}}}$

Esse foi apenas um exemplo simples, com duas reações. Sistemas mais complexos com mais reações são tratados da mesma maneira. Todas as taxas de reação devem ser calculadas em cada etapa de tempo e escolhida com probabilidade igual à sua contribuição fracionária para a taxa. O tempo é então avançado como neste exemplo.

A epidemia SIR sem dinâmica vital

O modelo SIR é uma descrição biológica clássica de como certas doenças permeiam uma população de tamanho fixo. Em sua forma mais simples, existem membros da população, em que cada membro pode estar em um dos três estados - suscetível, infectado ou recuperado - a qualquer momento, e cada um desses membros transita irreversivelmente através desses estados de acordo com o gráfico direcionado abaixo . Podemos denotar o número de membros suscetíveis como , o número de membros infectados como e o número de membros recuperados como . Portanto, para qualquer momento. ${\ displaystyle N}$${\ displaystyle n_ {S}}$${\ displaystyle n_ {I}}$${\ displaystyle n_ {R}}$${\ displaystyle N = n_ {S} + n_ {I} + n_ {R}}$

Além disso, um determinado membro suscetível fará a transição para o estado infectado ao entrar em contato com qualquer um dos membros infectados, de modo que a infecção ocorre com frequência . Um determinado membro do estado infectado se recupera sem dependência de nenhum dos três estados, o que é especificado pela taxa β . Dado este esquema básico, é possível construir o seguinte sistema não linear. ${\ displaystyle n_ {I}}$${\ displaystyle \ alpha n_ {I}}$

Trajetórias produzidas pelo método direto da epidemia de SIR (linhas laranja), com uma solução numérica (linhas pretas) sobreposta.

${\ displaystyle {\ frac {dn_ {S}} {dt}} = - {\ frac {\ alpha n_ {S}} {N}} n_ {I}}$,

${\ displaystyle {\ frac {dn_ {I}} {dt}} = \ left ({\ frac {\ alpha n_ {S}} {N}} - \ beta \ right) n_ {I}}$,

${\ displaystyle {\ frac {dn_ {R}} {dt}} = \ beta n_ {I}}$.

Este sistema não possui solução analítica. Uma abordagem poderia ser usar o algoritmo de Gillespie, simular o sistema muitas vezes e usar uma técnica de regressão, como mínimos quadrados, para ajustar um polinômio em todas as trajetórias. Conforme o número de trajetórias aumenta, tal regressão polinomial se comportará assintoticamente como a solução numérica (linhas pretas). Além de estimar a solução de um problema intratável como a epidemia de SIR, a natureza estocástica de cada trajetória permite calcular outras estatísticas além de . Ao executar o script abaixo, às vezes são obtidas realizações da epidemia de SIR que diferem drasticamente da solução numérica. Por exemplo, quando todas as pessoas são curadas (por acaso) muito cedo ou muito tarde. ${\ displaystyle \ mathrm {E} [n | t]}$

As trajetórias apresentadas na figura acima foram obtidas com a seguinte implementação Python do método direto, junto com uma classe de modelo cujos membros interagem com a maquinaria do método direto para realizar simulações gerais com as teorias subjacentes ao algoritmo de Gillespie. Estes são apresentados acima. Além disso, um solucionador ODE de SciPy é invocado para obter a solução numérica do sistema de equações diferenciais, ou seja, uma representação de . ${\ displaystyle \ mathrm {E} [n | t]}$

from matplotlib import pyplot
from numpy import linspace
from scipy.integrate import odeint

# initial species counts and sojourn times
initital_conditions = {
    "s": [480],
    "i": [20],
    "r": [0],
    "time": [0.0],
}

# propensity functions
propensities = {
    0: lambda d: 2.0 * d["s"][-1] * d["i"][-1] / 500,
    1: lambda d: 1.0 * d["i"][-1],
}

# change in species for each propensity
stoichiometry = {
    0: {"s": -1, "i": 1, "r": 0},
    1: {"s": 0, "i": -1, "r": 1},
}

# instantiate the epidemic SSA model container
epidemic = SSAModel(
    initital_conditions,
    propensities,
    stoichiometry
)

# instantiate the SSA container with model
epidemic_generator = SSA(epidemic)

# make a nice, big figure
pyplot.figure(figsize=(10,10), dpi=500)

# make a subplot for the susceptible, infected and recovered individuals
axes_s = pyplot.subplot(311)
axes_s.set_ylabel("susceptible individuals")

axes_i = pyplot.subplot(312)
axes_i.set_ylabel("infected individuals")

axes_r = pyplot.subplot(313)
axes_r.set_ylabel("recovered individuals")
axes_r.set_xlabel("time (arbitrary units)")

# simulate and plot 30 trajectories
trajectories = 0
for trajectory in epidemic_generator.direct():
    axes_s.plot(trajectory["time"], trajectory["s"], color="orange")
    axes_i.plot(trajectory["time"], trajectory["i"], color="orange")
    axes_r.plot(trajectory["time"], trajectory["r"], color="orange")
    trajectories += 1
    if trajectories == 30:
        break

# numerical solution using an ordinary differential equation solversir
t = linspace(0, 14, num=200)
y0 = (480, 20, 0)
alpha = 2.0
beta = 1.0

def differential_SIR(n_SIR, t, alpha, beta):
    dS_dt = -alpha * n_SIR[0] * n_SIR[1] / 500
    dI_dt = ((alpha * n_SIR[0] / 500) - beta) * n_SIR[1]
    dR_dt = beta * n_SIR[1]
    return dS_dt, dI_dt, dR_dt

solution = odeint(differential_SIR, y0, t, args=(alpha, beta))
solution = [[row[i] for row in solution] for i in range(3)]

# plot numerical solution
axes_s.plot(t, solution[0], color="black")
axes_i.plot(t, solution[1], color="black")
axes_r.plot(t, solution[2], color="black")

pyplot.show()

Leitura adicional

• Gillespie, Daniel T. (1977). "Simulação estocástica exata de reações químicas acopladas". The Journal of Physical Chemistry . 81 (25): 2340–2361. CiteSeerX  10.1.1.704.7634 . doi : 10.1021 / j100540a008 .
• Gillespie, Daniel T. (1976). "Um método geral para simular numericamente a evolução do tempo estocástico de reações 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). "Simulação Estocástica Exata Eficiente de Sistemas Químicos com Muitas Espécies e Muitos Canais" (PDF) . Journal of Physical Chemistry Uma . 104 (9): 1876–1889. Bibcode : 2000JPCA..104.1876G . doi : 10.1021 / jp993732q .
• Doob, Jacob L. (1942). "Tópicos da Teoria das Cadeias 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). "Cadeias de Markoff - Caso enumerável". 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). "Seção 17.7. Simulação Estocástica de Redes de Reações Químicas" . Receitas Numéricas: A Arte da Computação Científica (3ª ed.). New York, NY: Cambridge University Press. ISBN 978-0-521-88068-8.
• Kolmogorov, Andrey N. (1931). "Über die analytischen Methoden in der Wahrscheinlichkeitsrechnung" [Sobre métodos analíticos na teoria da probabilidade]. Mathematische Annalen . 104 : 415–458. doi : 10.1007 / BF01457949 . S2CID  119439925 .
• Feller, Willy (1940). "On the Integro-Differential Equations of Purely Discontinuous Markoff Processes" . Transactions of the American Mathematical Society . 48 (3): 4885–15. doi : 10.2307 / 1990095 . JSTOR  1970064 .
• Kendall, David G. (1950). "Uma realização artificial de um processo simples de" nascimento e morte ". Journal of Royal Statistical Society, Série B . 12 (1): 116-119. JSTOR  2983837 .
• Bartlett, Maurice S. (1953). "Processos estocásticos ou estatísticas de mudança". Journal of Royal Statistical Society, Série C . 2 (1): 44–64. JSTOR  2985327 .
• Rathinam, Muruhan; Petzold, Linda R .; Cao, Yang; Gillespie, Daniel T. (2003). "Rigidez em sistemas estocásticos de reação química: O 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). "Grosso-granulação adiabática e simulações de redes bioquímicas estocásticas" . Anais da Academia Nacional de Ciências dos Estados Unidos da América . 106 (20): 10546–10551. Bibcode : 2009PNAS..10610546S . doi : 10.1073 / pnas.0809340106 . PMC  2705573 . PMID  19525397 .
• Salis, Howard; Kaznessis, Yiannis N. (2005). "Simulação estocástica híbrida precisa de um sistema de reações químicas ou 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). "Um algoritmo de Monte Carlo cinético de tempo constante para simulação de grandes redes de reações 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; Apressado, Jeff; Tsimring, Lev S. (2005). Oscilações estocásticas induzidas por atraso na regulação do gene " . Anais da Academia Nacional de Ciências dos Estados Unidos da 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). "Regulação oscilatória de hes1 : Modelagem e Simulação de Atraso 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). "Simulação estocástica exata de reações químicas acopladas com atrasos". 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). Introdução à Modelagem para Biociências . Springer Verlag.
• (Ramaswamy González-Segredo Sbalzarini 2009): Ramaswamy, Rajesh; González-Segredo, Nélido; Sbalzarini, Ivo F. (2009). "Uma nova classe de algoritmos de simulação estocástica exata altamente eficiente para redes de reação química". 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). "Uma variante de propensão parcial do algoritmo de simulação estocástica de rejeição de composição para redes de reação química" (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.). "Fatores de reação e gráficos de atualização bipartida aceleram o algoritmo de Gillespie para sistemas bioquímicos em grande escala" . 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). "Uma formulação de propensão parcial do algoritmo de simulação estocástica para redes de reações químicas com atrasos" (PDF) . Journal of Chemical Physics . 134 (1): 014106. Bibcode : 2011JChPh.134a4106R . doi : 10.1063 / 1.3521496 . PMID  21218996 .
• (Yates Klingbeil 2013): Yates, Christian A .; Klingbeil, Guido (2013). "Reciclagem de números aleatórios no algoritmo de simulação estocástica" . Revisão Anual de Química Física . 58 (9): 094103. bibcode : 2013JChPh.138i4103Y . doi : 10.1063 / 1.4792207 . PMID  23485273 .
• Gillespie, Daniel T. (2007). "Simulação Estocástica de Cinética Química". The Journal of Physical Chemistry . 58 : 35–55. doi : 10.1146 / annurev.physchem.58.032806.104637 . PMID  17037977 .

links externos

• Um repositório contendo todo o código-fonte usado neste artigo