It’s been a while since I blogged (grant writing etc gets int he way of actual thinking and doing) so here is a quick post that uses Python to model a simple reaction A -> B using the Gillespie next reaction method.
I know importing pylab is discouraged but I can never remember how to import matplotlib. pylab is an easy way to get matplotlib. It also imports scipy and numpy but I import numpy separately.
I use lists instead of array to accumulate the results because they are faster to extend as the data comes in. Note that with a Gillespie algorithm you don’t know beforehand how many data points you’ll get.
I’ll also make a shameless plug to my book Enzyme Kinetics for Systems Biology which explains the algorithm.
import random, math, pylab
import numpy as np
# A -> B
Na = 100 # Number of molecules of A
Nb = 0 # Number of molecules of B
p = 0.99 # Probability of a reaction
timeEnd = 1000
t = 0
timeData = 
moleculesA = 
moleculesB = 
if Na == 0:
tau = -(1/(p*Na))*math.log (random.random())
Na -= 1
Nb += 1
t += tau
if t > timeEnd:
pylab.plot (timeData, moleculesA)
pylab.plot (timeData, moleculesB)