Algorithme de Gillespie - Gillespie algorithm

En théorie des probabilités , l' algorithme de Gillespie (ou occasionnellement l' algorithme de Doob-Gillespie ) génère une trajectoire statistiquement correcte (solution possible) d'un système d'équations stochastiques pour lequel les taux de réaction sont connus. Il a été créé par Joseph L. Doob et d'autres (vers 1945), présenté par Dan Gillespie en 1976, et popularisé en 1977 dans un article où il l'utilise pour simuler des systèmes chimiques ou biochimiques de réactions de manière efficace et précise en utilisant une puissance de calcul limitée (voir simulation stochastique ). Comme les ordinateurs sont devenus plus rapides, l'algorithme a été utilisé pour simuler des systèmes de plus en plus complexes. L'algorithme est particulièrement utile pour simuler des réactions au sein des cellules, où le nombre de réactifs est faible et où le suivi de la position et du comportement des molécules individuelles est réalisable par ordinateur. Mathématiquement, c'est une variante d'une méthode de Monte Carlo dynamique et similaire aux méthodes de Monte Carlo cinétiques . Il est largement utilisé dans la biologie des systèmes informatiques .

Histoire

Le processus qui a conduit à l'algorithme reconnaît plusieurs étapes importantes. En 1931, Andrei Kolmogorov a introduit les équations différentielles correspondant à l'évolution temporelle des processus stochastiques qui procèdent par sauts, aujourd'hui appelées équations de Kolmogorov (processus de saut de Markov) (une version simplifiée est connue sous le nom d' équation maîtresse en sciences naturelles). C'est William Feller , en 1940, qui a trouvé les conditions dans lesquelles les équations de Kolmogorov admettaient des probabilités (propres) comme solutions. Dans son théorème I (1940), il établit que le temps jusqu'au prochain saut était distribué de manière exponentielle et que la probabilité de l'événement suivant est proportionnelle au taux. En tant que tel, il a établi la relation des équations de Kolmogorov avec les processus stochastiques . Plus tard, Doob (1942, 1945) a étendu les solutions de Feller au-delà du cas des processus de saut pur. La méthode a été mise en œuvre dans des ordinateurs par David George Kendall (1950) à l'aide de l' ordinateur Manchester Mark 1 et plus tard utilisée par Maurice S. Bartlett (1953) dans ses études sur les épidémies. Gillespie (1977) obtient l'algorithme d'une manière différente en utilisant un argument physique.

Idée derrière l'algorithme

Les équations de vitesse biochimiques continues et déterministes traditionnelles ne prédisent pas avec précision les réactions cellulaires car elles reposent sur des réactions en vrac qui nécessitent les interactions de millions de molécules. Ils sont généralement modélisés comme un ensemble d'équations différentielles ordinaires couplées. En revanche, l'algorithme de Gillespie permet une simulation discrète et stochastique d'un système avec peu de réactifs car chaque réaction est explicitement simulée. Une trajectoire correspondant à une seule simulation de Gillespie représente un échantillon exact de la fonction de masse de probabilité qui est la solution de l' équation maîtresse .

La base physique de l'algorithme est la collision de molécules dans un récipient de réaction. On suppose que les collisions sont fréquentes, mais les collisions avec la bonne orientation et l'énergie sont rares. Par conséquent, toutes les réactions dans le cadre de Gillespie doivent impliquer au plus deux molécules. Les réactions impliquant trois molécules sont supposées être extrêmement rares et sont modélisées comme une séquence de réactions binaires. On suppose également que l'environnement réactionnel est bien mélangé.

Algorithme

Une revue récente (Gillespie, 2007) décrit trois formulations différentes mais équivalentes ; les méthodes directes, de première réaction et de première famille, les deux premières étant des cas particuliers de la seconde. La formulation des méthodes directe et de première réaction est centrée sur la réalisation des étapes habituelles d'inversion de Monte-Carlo sur la soi-disant "prémisse fondamentale de la cinétique chimique stochastique", qui est mathématiquement la fonction

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

où chacun des termes sont des fonctions de propension d'une réaction élémentaire, dont l'argument est , le vecteur des espèces compte. Le paramètre est le temps jusqu'à la prochaine réaction (ou temps de séjour), et est bien sûr l'heure actuelle. Pour paraphraser Gillespie, cette expression se lit comme « la probabilité, étant donné , que la prochaine réaction du système se produira dans l'intervalle de temps infinitésimal , et sera de stoechiométrie correspondant à la e réaction ». Cette formulation ouvre une fenêtre sur les méthodes directes et de première réaction en impliquant qu'il s'agit d'une variable aléatoire à distribution exponentielle et qu'il s'agit « d'une variable aléatoire entière statistiquement indépendante avec des probabilités ponctuelles ». ${\style d'affichage a}$${\displaystyle {\boldsymbol {x}}}$${\style d'affichage \tau }$${\style d'affichage t}$${\displaystyle {\boldsymbol {X}}(t)={\boldsymbol {x}}}$${\displaystyle [t+\tau ,t+\tau +d\tau ]}$${\style d'affichage j}$${\style d'affichage \tau }$${\style d'affichage j}$${\displaystyle a_{j}({\boldsymbol {x}})/\sum _{j}a_{j}({\boldsymbol {x}})}$

Ainsi, la méthode de génération de Monte-Carlo consiste simplement à tirer deux nombres pseudo-aléatoires, et sur , et à calculer ${\style d'affichage r_{1}}$${\style d'affichage r_{2}}$${\style d'affichage [0,1]}$

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

et

${\style d'affichage j=}$le plus petit entier satisfaisant .${\displaystyle \sum _{j'=1}^{j}a_{j'}({\boldsymbol {x}})>r_{2}\sum _{j}a_{j}({\boldsymbol { X}})}$

En utilisant cette méthode de génération pour le temps de séjour et la réaction suivante, l'algorithme de la méthode directe est indiqué par Gillespie comme

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.

Cette famille d'algorithmes est coûteuse en calcul et donc de nombreuses modifications et adaptations existent, y compris la méthode de réaction suivante (Gibson & Bruck), tau-leaping , ainsi que des techniques hybrides où des réactifs abondants sont modélisés avec un comportement déterministe. Les techniques adaptées compromettent généralement l'exactitude de la théorie derrière l'algorithme lorsqu'il se connecte à l'équation maîtresse, mais offrent des réalisations raisonnables pour des échelles de temps considérablement améliorées. Le coût de calcul des versions exactes de l'algorithme est déterminé par la classe de couplage du réseau de réaction. Dans les réseaux faiblement couplés, le nombre de réactions influencées par toute autre réaction est limité par une petite constante. Dans les réseaux fortement couplés, un seul tir de réaction peut en principe affecter toutes les autres réactions. Une version exacte de l'algorithme avec une mise à l'échelle en temps constant pour les réseaux faiblement couplés a été développée, permettant une simulation efficace de systèmes avec un très grand nombre de canaux de réaction (Slepoy Thompson Plimpton 2008). L'algorithme de Gillespie généralisé qui tient compte des propriétés non-markoviennes d'événements biochimiques aléatoires avec retard a été développé par Bratsun et al. 2005 et indépendamment Barrio et al. 2006, ainsi que (Cai 2007). Voir les articles cités ci-dessous pour plus de détails.

Les formulations à propension partielle, telles que développées indépendamment par Ramaswamy et al. (2009, 2010) et Indurkhya et Beal (2010), sont disponibles pour construire une famille de versions exactes de l'algorithme dont le coût de calcul est proportionnel au nombre d'espèces chimiques dans le réseau, plutôt qu'au nombre (plus important) de réactions. Ces formulations peuvent réduire le coût de calcul à une mise à l'échelle en temps constant pour les réseaux faiblement couplés et à une échelle au plus linéaire avec le nombre d'espèces pour les réseaux fortement couplés. Une variante à propension partielle de l'algorithme généralisé de Gillespie pour les réactions avec retard a également été proposée (Ramaswamy Sbalzarini 2011). L'utilisation des méthodes de propension partielle est limitée aux réactions chimiques élémentaires, c'est-à-dire aux réactions avec au plus deux réactifs différents. Chaque réaction chimique non élémentaire peut être décomposée de manière équivalente en un ensemble de réactions élémentaires, au prix d'une augmentation linéaire (dans l'ordre de la réaction) de la taille du réseau.

Considérez l'implémentation orientée objet suivante des méthodes de réaction directe et de première réaction, qui sont contenues dans une 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()

Les membres directet first_reactionsont des générateurs indéfinis, ce qui signifie qu'ils produiront en continu des trajectoires, chaque trajectoire étant une simulation complète du système, dans une boucle jusqu'à ce qu'un signal interrompe cette boucle. Une mise en garde importante pour toute mise en œuvre de cet algorithme est qu'il ne fournit aucune logique pour faire la distinction entre les événements élémentaires possibles et impossibles. Ainsi, pour utiliser réellement l'une ou l'autre version de l'algorithme dans une boucle et produire un certain nombre de trajectoires pour l'analyse, cette classe nécessite qu'un modèle lui soit transmis lors de l'instanciation. En plus de loger les populations d'espèces et les fonctions de propension, un modèle a trois méthodes auxquelles l'algorithme de Gillespie fait appel : curate, exit, et reset. Le but de ces méthodes est spécifié ci-dessous, et régissent essentiellement lorsque l'algorithme se termine. Le raisonnement derrière l'introduction de la classe de modèle est principalement d'éviter de mélanger les propriétés cinétiques du processus spécifique simulé avec la logique interne de l'algorithme de Gillespie.

Le modèle ci-dessous est une sous-classe de dictionnaire avec les membres publics curate, exit, et resetque respectivement

• déterminer quelles réactions sont et ne sont pas valides à la fin de chaque itération d'une trajectoire donnée ;
• revenir Trues'il n'y a plus de réactions possibles ;
• remettre les hachages du modèle à leurs valeurs d'origine (c'est-à-dire leurs conditions initiales) à la fin d'une trajectoire donnée.

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()

Exemples

Liaison réversible de A et B pour former des dimères AB

Un exemple simple peut aider à expliquer le fonctionnement de l'algorithme de Gillespie. Considérons un système de molécules de deux types, A et B . Dans ce système, A et B de manière réversible se lient ensemble pour former AB dimères de telle sorte que les deux réactions sont possibles: ou bien A et B réagissent de manière réversible pour former un AB dimère, ou un AB dimères se dissocie en A et B . La constante de vitesse de réaction pour une seule molécule A donnée réagissant avec une seule molécule B donnée est , et la vitesse de réaction pour la rupture d' un dimère AB est . ${\displaystyle k_{\mathrm {D} }}$${\displaystyle k_{\mathrm {B} }}$

Si au temps t il y a une molécule de chaque type alors le taux de formation de dimère est , tandis que s'il y a des molécules de type A et des molécules de type B , le taux de formation de dimère est . S'il y a des dimères, le taux de dissociation des dimères est de . ${\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 vitesse de réaction totale, , au temps t est alors donnée par ${\displaystyle R_{\mathrm {TOT} }}$

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

Ainsi, nous avons maintenant décrit un modèle simple avec deux réactions. Cette définition est indépendante de l'algorithme de Gillespie. Nous allons maintenant décrire comment appliquer l'algorithme de Gillespie à ce système.

Dans l'algorithme, nous avançons dans le temps en deux étapes : calculer le temps jusqu'à la prochaine réaction et déterminer laquelle des réactions possibles est la prochaine réaction. Les réactions sont supposées être complètement aléatoires, donc si la vitesse de réaction à un instant t est , alors le temps, t , jusqu'à ce que la prochaine réaction se produise est un nombre aléatoire tiré de la fonction de distribution exponentielle avec la moyenne . Ainsi, nous avançons le temps de t à t + t . ${\displaystyle R_{\mathrm {TOT} }}$${\displaystyle 1/R_{\mathrm {TOT} }}$

Tracé du nombre de molécules A (courbe noire) et des dimères AB en fonction du temps. Comme nous avons commencé avec 10 molécules A et B au temps t = 0, le nombre de molécules B est toujours égal au nombre de molécules A et il n'est donc pas représenté.

La probabilité que cette réaction soit une molécule A se liant à une molécule B est simplement la fraction du taux total due à ce type de réaction, c'est-à-dire,

la probabilité que la réaction soit ${\displaystyle P({\ce {{A}+ B -> AB}})=k_{D}n_{A}n_{B}/R_{{\ce {TOT}}}}$

La probabilité que la prochaine réaction soit un dimère AB qui se dissocie n'est que de 1 moins cela. Donc avec ces deux probabilités soit on forme un dimère en réduisant et de un, et en augmentant de un, soit on dissocie un dimère et on augmente et de un et on diminue de un. ${\displaystyle n_{\mathrm {A} }}$${\displaystyle n_{\mathrm {B} }}$${\displaystyle n_{\mathrm {AB} }}$${\displaystyle n_{\mathrm {A} }}$${\displaystyle n_{\mathrm {B} }}$${\displaystyle n_{\mathrm {AB} }}$

Maintenant, nous avons à la fois avancé le temps jusqu'à t + t et effectué une seule réaction. L'algorithme de Gillespie répète simplement ces deux étapes autant de fois que nécessaire pour simuler le système aussi longtemps que l'on veut (c'est-à-dire pour autant de réactions). Le résultat d'une simulation Gillespie qui commence par et à t =0, et où et , est affiché à droite. Pour ces valeurs de paramètres, il y a en moyenne 8 dimères et 2 de A et B mais en raison du petit nombre de molécules, les fluctuations autour de ces valeurs sont importantes. L'algorithme de Gillespie est souvent utilisé pour étudier des systèmes où ces fluctuations sont 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} }}$

C'était juste un exemple simple, avec deux réactions. Les systèmes plus complexes avec plus de réactions sont traités de la même manière. Toutes les vitesses de réaction doivent être calculées à chaque pas de temps, et une doit être choisie avec une probabilité égale à sa contribution fractionnaire à la vitesse. Le temps est ensuite avancé comme dans cet exemple.

L'épidémie de SIR sans dynamique vitale

Le modèle SIR est une description biologique classique de la façon dont certaines maladies imprègnent une population de taille fixe. Dans sa forme la plus simple, il y a des membres de la population, chaque membre pouvant se trouver dans l'un des trois états - sensible, infecté ou récupéré - à tout moment, et chacun de ces membres passe de manière irréversible à travers ces états selon le graphique ci-dessous. . Nous pouvons noter le nombre de membres sensibles comme , le nombre de membres infectés comme , et le nombre de membres récupérés comme . Donc pour n'importe quel moment. ${\style d'affichage N}$${\style d'affichage n_{S}}$${\style d'affichage n_{I}}$${\style d'affichage n_{R}}$${\style d'affichage N=n_{S}+n_{I}+n_{R}}$

De plus, un membre sensible donné passera à l'état infecté en entrant en contact avec l'un des membres infectés, de sorte que l'infection se produit avec le taux . Un membre donné de l'état infecté se rétablit sans dépendre de l'un des trois états, ce qui est spécifié par le taux β . Compte tenu de ce schéma de base, il est possible de construire le système non linéaire suivant. ${\style d'affichage n_{I}}$${\displaystyle \alpha n_{I}}$

Trajectoires telles que données par la méthode directe de l'épidémie SIR (lignes oranges), avec une solution numérique (lignes noires) superposées.

${\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}}$.

Ce système n'a pas de solution analytique. Une approche pourrait consister à utiliser l'algorithme de Gillespie, à simuler le système plusieurs fois et à utiliser une technique de régression telle que les moindres carrés pour ajuster un polynôme sur toutes les trajectoires. Au fur et à mesure que le nombre de trajectoires augmente, une telle régression polynomiale se comportera asymptotiquement comme la solution numérique (lignes noires). En plus d'estimer la solution d'un problème insoluble comme l'épidémie de SIR, la nature stochastique de chaque trajectoire permet de calculer des statistiques autres que . Lors de l'exécution du script ci-dessous, des réalisations de l'épidémie de SIR sont parfois obtenues qui diffèrent considérablement de la solution numérique. Par exemple, lorsque toutes les personnes sont guéries (par hasard) très tôt, ou très tard. ${\displaystyle \mathrm {E} [n|t]}$

Les trajectoires présentées dans la figure ci-dessus ont été obtenues avec l'implémentation Python suivante de la méthode directe, ainsi qu'une classe de modèle dont les membres interagissent avec la machinerie de la méthode directe pour effectuer des simulations générales avec les théories sous-jacentes à l'algorithme de Gillespie. Ceux-ci sont présentés ci-dessus. De plus, un solveur ODE de SciPy est invoqué pour obtenir la solution numérique du système d'équations différentielles, c'est-à-dire une représentation 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()

Lectures complémentaires

• Gillespie, Daniel T. (1977). « Simulation stochastique exacte de réactions chimiques couplées ». Le Journal de Chimie Physique . 81 (25) : 2340–2361. CiteSeerX  10.1.1.704.7634 . doi : 10.1021/j100540a008 .
• Gillespie, Daniel T. (1976). « Une méthode générale pour simuler numériquement l'évolution temporelle stochastique des réactions chimiques couplées ». Journal de physique computationnelle . 22 (4) : 403-434. Bibcode : 1976JCoPh..22..403G . doi : 10.1016/0021-9991 (76) 90041-3 .
• Gibson, Michael A.; Bruck, Josué (2000). « Simulation stochastique exacte efficace de systèmes chimiques avec de nombreuses espèces et de nombreux canaux » (PDF) . Journal de chimie physique A . 104 (9) : 1876-1889. Bibcode : 2000JPCA..104.1876G . doi : 10.1021/jp993732q .
• Doob, Jacob L. (1942). « Sujets dans la théorie des chaînes de Markoff » . Transactions de l'American Mathematical Society . 52 (1) : 37-64. doi : 10.1090/S0002-9947-1942-0006633-7 . JSTOR  1990152 .
• Doob, Jacob L. (1945). "Chaînes Markoff – Cas dénombrable". Transactions de l'American Mathematical Society . 58 (3) : 455-473. doi : 10.2307/1990339 . JSTOR  1990339 .
• Presse, William H. ; Teukolsky, Saul A.; Vetterling, William T. ; Flannery, Brian P. (2007). "Section 17.7. Simulation stochastique des réseaux de réaction chimique" . Recettes numériques : L'art de l'informatique scientifique (3e éd.). New York, NY : Cambridge University Press. ISBN 978-0-521-88068-8.
• Kolmogorov, Andrey N. (1931). "Über die analytischen Methoden in der Wahrscheinlichkeitsrechnung" [Sur les méthodes analytiques dans la théorie des probabilités]. Mathematische Annalen . 104 : 415-458. doi : 10.1007/BF01457949 . S2CID  119439925 .
• Feller, Willy (1940). "Sur les équations intégro-différentielles des processus de marquage purement discontinus" . Transactions de l'American Mathematical Society . 48 (3) : 4885-15. doi : 10.2307/1990095 . JSTOR  1970064 .
• Kendall, David G. (1950). "Une réalisation artificielle d'un processus simple de "naissance et mort"". Journal de la Royal Statistical Society, série B . 12 (1) : 116-119. JSTOR  2983837 .
• Bartlett, Maurice S. (1953). "Les processus stochastiques ou les statistiques du changement". Journal de la Royal Statistical Society, série C . 2 (1) : 44-64. JSTOR  2985327 .
• Rathinam, Muruhan; Petzold, Linda R. ; Cao, Yang ; Gillespie, Daniel T. (2003). « Raideur dans les systèmes à réaction chimique stochastique : la méthode implicite de saut de tau ». Journal de physique chimique . 119 (24) : 12784-12794. Bibcode : 2003JChPh.11912784R . doi : 10.1063/1.1627296 .
• Sinitsyne, Nikolaï A. ; Hengartner, Nicolas ; Nemenman, Ilya (2009). "Gros grain adiabatique et simulations de réseaux biochimiques stochastiques" . Actes de l'Académie nationale des sciences des États-Unis d'Amérique . 106 (20) : 10546–10551. Bibcode : 2009PNAS..10610546S . doi : 10.1073/pnas.0809340106 . PMC  2705573 . PMID  19525397 .
• Salis, Howard ; Kaznessis, Yiannis N. (2005). « Simulation stochastique hybride précise d'un système de réactions chimiques ou biochimiques couplées ». Journal de physique chimique . 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 algorithme de Monte Carlo cinétique à temps constant pour la simulation de grands réseaux de réactions biochimiques ». Journal de physique chimique . 128 (20): 205101. bibcode : 2008JChPh.128t5101S . doi : 10.1063/1.2919546 . PMID  18513044 .
• (Bratsun et al. 2005) : Bratsun, Dmitri ; Volfson, Dmitri ; Hâte, Jeff ; Tsimring, Lev S. (2005). « Oscillations stochastiques induites par le retard dans la régulation des gènes » . Actes de l'Académie nationale des sciences des États-Unis d'Amérique . 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). " Régulation oscillatoire de hes1 : Modélisation et simulation de retard stochastique discret " . Biologie computationnelle PLOS . 2 (9): 1017. bibcode : 2006PLSCB ... 2..117B . doi : 10.1371/journal.pcbi.0020117 . PMC  1560403 . PMID  16965175 .
• (Cai 2007) : Cai, Xiaodong (2007). « Simulation stochastique exacte de réactions chimiques couplées avec des retards ». Journal de physique chimique . 126 (12): 124108. bibcode : 2007JChPh.126l4108C . doi : 10.1063/1.2710253 . PMID  17411109 .
• (Barnes Chu 2010) : Barnes, David J. ; Chu, Dominique (2010). Introduction à la modélisation pour les biosciences . Springer Verlag.
• (Ramaswamy González-Segredo Sbalzarini 2009) : Ramaswamy, Rajesh ; González-Segredo, Nélido; Sbalzarini, Ivo F. (2009). « Une nouvelle classe d'algorithmes de simulation stochastique exacte très efficace pour les réseaux de réaction chimique ». Journal de physique chimique . 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). « Une variante de propension partielle de l'algorithme de simulation stochastique de rejet de composition pour les réseaux de réaction chimique » (PDF) . Journal de physique chimique . 132 (4): 044102. bibcode : 2010JChPh.132d4102R . doi : 10.1063/1.3297948 . PMID  20113014 .
• (Indurkhya Beal 2010) : Indurkhya, Sagar ; Beal, Jacob S. (2005). Isalan, Mark (éd.). "La factorisation de réaction et les graphiques de mise à jour bipartite accélèrent l'algorithme de Gillespie pour les systèmes biochimiques à grande échelle" . PLOS UN . 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). « Une formulation de propension partielle de l'algorithme de simulation stochastique pour les réseaux de réaction chimique avec des retards » (PDF) . Journal de physique chimique . 134 (1): 014106. bibcode : 2011JChPh.134a4106R . doi : 10.1063/1.3521496 . PMID  21218996 .
• (Yates Klingbeil 2013) : Yates, Christian A. ; Klingbeil, Guido (2013). "Recyclage des nombres aléatoires dans l'algorithme de simulation stochastique" . Revue annuelle de chimie physique . 58 (9): 094103. bibcode : 2013JChPh.138i4103Y . doi : 10.1063/1.4792207 . PMID  23485273 .
• Gillespie, Daniel T. (2007). "Simulation Stochastique de Cinétique Chimique". Le Journal de Chimie Physique . 58 : 35-55. doi : 10.1146/annurev.physchem.58.032806.104637 . PMID  17037977 .

Liens externes

• Un référentiel contenant tout le code source utilisé dans cet article