18 juin

Implémentation de l’algorithme de Gillespie : PseudoCode

Dans le précédent article, j’expliquais les bases de l’algorithme de Gillespie. Je vous propose maintenant une implémentation en pseudo-code. Je reprends volontairement les étapes de la même façon que Gillespie dans son article A general method for Numerically simulating the stochastic time evolution of coupled chemical reactions.

Step 0 Initialisation

Variables :

Fixer le temps : t = 0
Fixer les concentrations initiales des N espèces : X1 X2 X3... XN
Spécifier et stocker les valeurs c1 c2 c3 ... cM pour les M réactions
Calculer les M quantité c1h1 c2h2 c3h3 ... cMhM pour les M réactions

Spécifier et stocker une série de temps auxquels on va générer
 une sortie des concentrations XN.

Fixer un temps maximum tstop.

Step 1

Générer une paire aléatoire Tau et Mu, tel que présenté dans le précédent article.

r2 = uniforme(0,1)
suma = a1 + ... + an
i = 1
sumpart = a1
tant que r2 * a < sumpart :
    incrémente i
    sumpart = sumpart + ai
sinon
    mu = i - 1
r1 = uniforme(0,1)
tau = (1/suma) ln(1/r1)

Step 2

En utilisant les nombres Tau et Mu, on fait avancer t de Tau et on change les valeurs de quantités Xi impliquées dans la réaction Mu.

pour i dans 1, n n nombre de réactants de Rmu
Xi = Xi - ai avec ai le nombre stoechiométrique correspondant
pour i dans 1, m m nombre de produits de Rmu
Xi = Xi + ai
t = t + tau

Step 3

Si t vient de dépasser un des temps auxquels on doit générer une sortie, lire les quantités de toutes les concentrations Xi. Si t est supérieur à tstop ou si tous les h sont nuls (preuve qu’il n’y a plus de réactants), s’interrompre, sinon retourner à l’étape 1.