Simple Stochastic Code in Python

Originally Posted on February 5, 2019 by hsauro

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.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

import random, math, pylab import numpy as np   # Model: #   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 = [] while True:    if Na == 0:        break    tau = -(1/(p*Na))*math.log (random.random())    print (tau)    Na -= 1    Nb += 1    t += tau    timeData.append (t)    moleculesA.append (Na)    moleculesB.append (Nb)    if t > timeEnd:       break    pylab.xlabel ('Time') pylab.plot (timeData, moleculesA) pylab.plot (timeData, moleculesB)