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:

 

How do I change the simulation tolerances in Tellurium?

November 4, 2016 8:59 am

For very complicated and large models it may be necessary to adjust the simulator tolerances in order to get the correct simulation results. Sometimes the simulator will terminate a simulation because it was unable to proceed due to numerical errors. In many cases this is due to a bad model and the user must investgate the model to determine what the issue might be. If the model is assumed to be correct then the other option is to change the simulator tolerances. The current option state of the simulator is obtained using the getInfo call, for example:

r = roadrunner.RoadRunner ('mymodel.xml')
r.getInfo()
roadrunner.RoadRunner() { 
'this' : 22F59350
'modelLoaded' : true
'modelName' : __main
'libSBMLVersion' : LibSBML Version: 5.13.0
'jacobianStepSize' : 1e-005
'conservedMoietyAnalysis' : false
'simulateOptions' : 
< roadrunner.SimulateOptions() 
{ 
'this' : 1B89AF28, 
'reset' : 0,
'structuredResult' : 0,
'copyResult' : 1,
'steps' : 50,
'start' : 0,
'duration' : 20
}, 
'integrator' : 
< roadrunner.Integrator() >
  name: cvode
  settings:
      relative_tolerance: 0.00001
      absolute_tolerance: 0.0000000001
                   stiff: true
       maximum_bdf_order: 5
     maximum_adams_order: 12
       maximum_num_steps: 20000
       maximum_time_step: 0
       minimum_time_step: 0
       initial_time_step: 0
          multiple_steps: false
      variable_step_size: false

}>

There are a variety of tuning parameters that can be changed in the simulator. Of interest are the relative and absolute tolerances, the maximum number of steps, and the initial time step.

The smaller the relative tolerance the more accurate the solution, however too small a value will result in either excessive runtimes or more likely roundoff errors. A relative tolerance of 1E-4 means that errors are controlled to 0.01%. An optimal value is roughly 1E-6. The absolute tolerance is used when a variable gets so small that the relative tolerance doesn't make much sense to apply. In these situations, absolute error tolerance is used to control the error. A small value for the absolute tolerance is often desirable, such as 1E-12, we do not recommend going below 1E-15 for either tolerance.

To set the tolerances use the statements:

r.integrator.absolute_tolerance = 5e-10
r.integrator.relative_tolerance = 1e-3

Another parameter worth changing if the simulations are not working well is to change the initial time step. This is often set by the integrator to be a relatively large value which means that the integrator will try to reduce this value if there are problems. Sometimes it is better to provide a small initial step size to help the integrator get started, for example, 1E-5.

r.integrator.initial_time_step = 0.00001

The reader is referred to the CVODE documentation for more details.

How do I plot phase plots using Tellurium?

 Phase plots are a common way to visualize the dynamics of models where time courses are generated and one variable is plotted against the other. For example, consider the following model that can show oscillations:

  v1: $Xo -> S1; k1*Xo;
  v2:  S1 -> S2; k2*S1*S2^h/(10 + S2^h) + k3*S1;
  v3:  S2 -> ;   k4*S2;      

In this model S2 positively activates reaction v2 thus forming a positive feedback loop. The rate equation for v2 include a Hill like coefficient term, S2^h, which determines the strength of the positive feedback. The oscillations originate from an interaction between the positive feedback and a non-obvious negative feedback loop at S1 between v1 and v2.

Let us assign suitable parameter values to this model, run a simulation and plot S1 versus S2.

import tellurium as te
# Import pylab to access subplot plotting feature.
import pylab as plt

r = te.loada ('''
  v1: $Xo -> S1; k1*Xo;
  v2:  S1 -> S2; k2*S1*S2^h/(10 + S2^h) + k3*S1;
  v3:  S2 -> $w; k4*S2;      

  # Initialize
  h  = 2; # Hill coefficient
  k1 = 1; k2 = 2; Xo = 1;
  k3 = 0.02; k4 = 1;
  S1 = 6; S2 = 2
''')

m = r.simulate (0, 80, 500, ['S1', 'S2'])
r.plot()

Running this script by clicking the green button in the toolbar yields the following plot:

What if we'd like to investigate how the oscillations are affected by the parameters of the model. For example, how does the model behave when we change k1? One way to do this is to plot simulations at different k1 values onto the same plot. In this case, however, this will create a difficult to read graph. Instead, let us create a grid of subplots where each subplot represents a different simulation.

import tellurium as te
import pylab as plt

r = te.loada ('''
  v1: $Xo -> S1; k1*Xo;
  v2:  S1 -> S2; k2*S1*S2^h/(10 + S2^h) + k3*S1;
  v3:  S2 -> $w; k4*S2;      

  # Initialize
  h  = 2; # Hill coefficient
  k1 = 1; k2 = 2; Xo = 1;
  k3 = 0.02; k4 = 1;
  S1 = 6; S2 = 2
''')

plt.subplots(3,3,figsize=(12,6)) 
for i in range (9):
    r.reset()
    m = r.simulate (0, 80, 500, ['S1', 'S2'])
    plt.subplot (3,3,i+1)
    plt.plot (m[:,0], m[:,1], label="k1=" + str (r.k1))
    plt.legend()
    r.k1 = r.k1 + 0.2

Here we create a 3 by 3 subplot grid, start a loop that changes k1, and each time round the loop it plots the simulation onto one of the subplots. Running this script results in the following output

 

How to get the Stoichiometry Matrix using Tellurium

 Here is a simple need, given a reaction model how do we get hold of the stoichiometry matrix?

Consider the following simple model:

import tellurium as te
import roadrunner

r = te.loada("""
     $Xo -> S1; k1*Xo;
      S1 -> S2; k2*S1;
      S2 -> S3; k3*S2;
      S3 -> $X1; k4*S3;
      
      k1 = 0.1; k2 = 0.4;
      k3 = 0.5; k4 = 0.6;
      Xo = 1;
""")

print r.getFullStoichiometryMatrix()

Running this script by clicking on the green arrow in the toolbar will yield:

      _J0, _J1, _J2, _J3
S1 [[   1,  -1,   0,   0],
S2  [   0,   1,  -1,   0],
S3  [   0,   0,   1,  -1]]

The nice thing about this output is that the columns and rows are labeled so you know exact what is what.

What about much larger models? For example the iAF1260.xml model from the Bigg database (http://bigg.ucsd.edu:8888/models/iAF1260). This is a model of E. Coli that includes 1668 metabolites and 2382 reactions. We can download the iAF1260.xml file and load it into libRoadRunner using:

 
r = roadrunner.RoadRunner ('iAF1260.xml')

This might take up to a minute to load depending on how fast your computer is. We are assuming here that the file is located in the current directory (os.getcwd()). If not, move the file, change the current directory (using os.chdir), or use the appropriate path in the call.

Rather than print out the stoichiometry matrix (don't even try) to the screen, we'll save it to a file. Because the stoichiometry matrix is so large we will use numpy to write the matrix out as a text file:

import numpy as np
r = roadrunner.RoadRunner ('iAF1260.xml')
st = r.getFullStoichiometryMatrix()
print "Number of metabolites = ", r.getNumFloatingSpecies()
print "Number of reactions = ", r.getNumReactions()
np.savetxt ('stoich.txt', st)

Number of metabolites = 1668
Number of reactions = 2382

One can change the formatting of the output using savetxt, for example, the following will output the individual stoichiometry coefficient using 3 decimal places, 5 characters minimum, and separated by a comma.

np.savetxt ('c:\\tmp\\st.txt', st, delimiter=',', fmt='%5.3f',)

You can get the labels for the rows and columns by calling r.getFloatingSpeciesIds() and r.getReactionIds() respectively.

How to do a simple simulation using Tellurium

November 1, 2016 1:10 pm

The most common requirement is the ability to carry out a simple time course simulation of a model. Consider the model:

$$ S_1 \rightarrow S_2 \rightarrow S_2 $$

Two reactions and three metabolites, S1, S2 and S3. We can describe this system using an Antimony string:

import tellurium as te
import roadrunner

r = te.loada ('''
   S1 -> S2; k1*S1;
   S2 -> S3; k2*S2;

   k1= 0.4; k2 = 0.45
   S1 = 5; S2 = 0; S3 = 0 
''')

m = r.simulate (0, 20, 100)
r.plot()

Run the script by clicking on the green arrow in the tool bar to yield: