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)

How to do a simple parameter scan using Tellurium

Originally Posted on August 27, 2018 by hsauro

A common task in modeling is to see how a parameter influences a model’s dynamics. For example, consider a simple two reaction pathway:

-> S1 ->

where the first reaction has a fixed input of vo and the second reaction a first-order rate laws k1*S1. The task is to investigate how the time course of S1 is influenced by vo.

The script below defines the model, then changes vo in increments and plots the effect on the pathway via a time course simulation. To do the parameter scan we exploit plotArray. This initially prevents the plots from being shown using show=False. We make sure that each plot gets a different color using resetColorCycle=False, finally, we show the plot using show(). To make things more interesting we also add a legend entry for each plot.

Note we call reset each time we run a simulation to ensure that S1 is reset back to its initial condition.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

import tellurium as te import numpy as np   r = te.loada("""         J1: -> S1; vo;         S1 -> ; k1*S1         k1 = 0.5;         vo = 4; """)   for vo in np.arange (0, 5, 0.5):     r.reset()     r.vo = vo;     m = r.simulate(0, 10, 100)     te.plotArray(m, labels=['v0 = '+str(vo)], show=False, resetColorCycle=False) te.show()

The following figure shows the resulting plot:

Thanks to Kiri Choi for pointing out how to use plotArray in this way.