How do I run a stochastic simulation using Tellurium/libRoadRunner?

 In this post I will show you how to run a stochastic simulation using our Tellurium application. Tellurium is a set of libraries that can be used via Python. One of those libraries is libRoadRunner which is our very fast simulator. It can simulate both stochastic and deterministic models. Let's illustrate a stochastic simulation using the following simple model:

import tellurium as te
import numpy as np

r = te.loada('''
    J1: S1 -> S2;  k1*S1;
    J2: S2 -> S3;  k2*S2
    J3: S3 -> S4;  k3*S3;

    k1 = 0.1; k2 = 0.5; k3 = 0.5;
    S1 = 100;
''')

We've set up the number of initial molecules of S1 to be 100 molecules. The easiest way to run a stochastic simulation is to call the gillespie method on roadrunner. This is shown in the code below:

m = r.gillespie (0, 40, steps=100)
r.plot()

Running this by clicking on the green button in the tool bar will give you the plot:

What if you wanted to run a lot of gillespie simulations to get an idea of the distribution of trajectories? TO do that just just need to repeat the simulation many times and plot all the results on the same graph:

import tellurium as te
import numpy as np

r = te.loada('''
    J1: S1 -> S2; k1*S1;
    J2: S2 -> S3; k2*S2
    J3: S3 -> S4; k3*S3;

    k1 = 0.1; k2 = 0.5; k3 = 0.5;
    S1 = 100;
''')

# run repeated simulation
numSimulations = 50
points = 101
for k in range(numSimulations ):
    r.resetToOrigin()
    s = r.gillespie(0, 50)
    # No legend, do not show
    r.plot(s, show=False, loc=None, alpha=0.4, linewidth=1.0)

This script will yield the plot:

We can do one other thing and compute the average trajectory and overlay the plot with the average line. The one thing we have to watch out for is that we must set the integrator property variable_step_size = False to false. This will ensure that time points are equally spaced and that all trajectories end at the same point in time.

import tellurium as te
import numpy as np

r = te.loada('''
    J1: S1 -> S2; k1*S1;
    J2: S2 -> S3; k2*S2
    J3: S3 -> S4; k3*S3;

    k1 = 0.1; k2 = 0.5; k3 = 0.5;
    S1 = 100;
''')

# run repeated simulation
numSimulations = 50
#points = 101
# Set the physical size of the plot (units are in inches)
plt.figure(figsize=(10,5))
r.setIntegrator ('gillespie')
# Make sure we do this so that all trajectories 
# are the same length and spacings
r.getIntegrator().variable_step_size = False
s_sum = np.array(0.)
for k in range(numSimulations):
    r.resetToOrigin()
    s = r.simulate(0, 100, steps=50)
    s_sum = np.add (s_sum, s)
    # no legend, do not show
    r.plot(s, show=False, loc=None, alpha=0.4, linewidth=1.0)

# add mean curve, legend, show everything and set labels, titels, ...
s_mean = s_sum/numSimulations
r.plot(s_mean, loc='upper right', show=True, linewidth=3.0,
       title="Stochastic simulation", xlabel="time", 
       ylabel="concentration", grid=True); 

This will give us the following plot: