Computing the steady-state using Tellurium

 If you're building a model and you want to quickly find the model's steady state, you can call the command steadyState. Let's illustrate this with an example:

import tellurium as te

r = te.loada ('''
     # Define a simple linear chain of reactions
     $Xo -> S1; k1*Xo;
      S1 -> S2; k2*S1;
      S2 -> S3; k3*S2;
      S3 -> $X1; k4*S3;

  # Initialize rate constants 
  k1 = 0.2; k2 = 0.7; k3 = 0.15; k4 = 1.3;
  Xo = 10
''')

print r.getSteadyStateValues() 

Running this script by clicking in the green arrow in the tool bar will output the following steady state levels:

[2.85714286  13.33333333   1.53846154]

But which steady state value refers to which species in my model? To find that out call getFloatingSpeciesIds:

print r.getFloatingSpeciesIds()
['S1', 'S2', 'S3']

The order of the speces Ids match the order of steady state values. In other words, S1 = 2.857, S2 = 13.333, and S3 = 1.538.

If you want to check that these values really are the steady state values you can compute the rates of change:

print r.getRatesOfChange()
array([  0.00000000e+00,  -4.44089210e-16,   4.44089210e-16])

If you look carefully all the rates of change are effectively zero, thus confirming we're at steady state.

What about the stability of the steady state? That is, is it stable to small perturbations in S1, S2 or S3? To find this out we need to compute the eigenvalues of the system Jacobian matrix. If all the eigenvalues are negative this means small pertubrations are damped so that the system returns to the steady state. There are two ways to do this, we can use the getReducedEigenValues call, or we get compute the Jacobian first and then compute the eigenvalues of the Jacobian. Both these approaches are given below. Its probably simplest to call just getReducedEigenValues.

print r.getReducedEigenValues()
[-0.7  -0.15 -1.3 ]

Notice that all the eigenvalues are negative indicating that small perturbations are stable.

This is the alternative approach which computes the Jacobian first and then the eigenvalues of the Jacobian:

print te.getEigenvalues(r.getFullJacobian())
[-0.7  -0.15 -1.3 ]

Note that we're using a utility method from the tellurium library, getEigenValues to compute the eigenvalues.

Bifurcation Analysis with Tellurium

Originally Posted on October 28, 2016 by hsauro

I thought I’d try and write a series of HowTos on Tellurium, our python-based tool for the construction, simulation and analysis of biochemical models. Details on this tool can be found here.

One of the unique features of Tellurium is that it comes with the AUTO2000 package. This is a well-established software library for performing a bifurcation analysis.

Unlike other implementations (not including oscill8), AUTO2000 in Tellurium does not require a separate compiler to compile the model for AUTO2000 to compute the differential equations. This makes it easier to deploy and users don’t have to worry about such details. AUTO2000 uses libRoadRunner to access and compute the model.

Let’s begin with an example from my textbook “Systems Biology: Introduction to Pathway Modeling”, Figure 12.20, page 279 in revision 1.1 of the book. The model in question is a modified model from the ‘Mutual activation’ model in the review by Tyson (Current Opinion in Cell Biology 15:221–231, Figure 1e). In this example increasing the signal results in the system switching to the high state at around 2.0. If we reduce the signal from a high level, we traverse a different steady-state. If we assume the signal can never be negative, we will remain at the high steady-state even if the signal is reduced to zero. The bifurcation plot in the negative quadrant of the graph is physically inaccessible. This means it is not possible to transition to the low steady-state by decreasing signal. As a result, the bistable system is irreversible, that is, once it is switched on, it will always remain on. To compute the bifurcation diagram we first define the model:

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

import tellurium as te import rrplugins import pylab as plt   r = te.loada ('''     $X ->; R1;  k1*EP + k2*Signal;     R1 ->; $w;  k3*R1;     EP ->; E;   Vm1*EP/(Km + EP);     E -> EP;   ((Vm2+R1)*E)/(Km + E);     Vm1 = 12; Vm2 = 6;     Km = 0.6;     k1 = 1.6; k2 = 4;     E = 5; EP = 15;     k3 = 3; Signal = 0.1; ''')

We’ve imported three packages, tellurium to load the model, rrplugins to access AUTO2000 and pylab to gain access to matplotlib. Once we have the model loaded we can get a handle on AUTO2000 by calling rrplugins.Plugin(“tel_auto2000”) and set a number of properties in AUTO2000. This includes loading the model into AUTO200, identifying the parameter we wish to modify for the bifurcation diagram (in this case signal), following by some options to carry out a pre simulation to help with the initial location of the steady-state and finally the limits for x axis for the plot, in this case -2 to 3. Details of other properties to change can be found by typing auto.viewManual(), make sure you have a pdf reader available. The alternative is to go to the intro page.

1 2 3 4 5 6 7 8 9

auto = rrplugins.Plugin("tel_auto2000")   auto.setProperty("SBML", r.getCurrentSBML()) auto.setProperty("ScanDirection", "Negative") auto.setProperty("PrincipalContinuationParameter", "Signal") auto.setProperty("PreSimulation", "True") auto.setProperty("PreSimulationDuration", 1.0) auto.setProperty("RL0", -2.0) auto.setProperty("RL1", 3.0)

To run the bifurcation analysis we use the Python code:

1	auto.execute()

If all was successful we can next plot the results. It is possible to plot the results using your own code (see below) but it might be more convenient to use the builtin facilties, for example:

1 2 3 4

pts     = auto.BifurcationPoints lbls    = auto.BifurcationLabels biData  = auto.BifurcationData biData.plotBifurcationDiagram(pts, lbls)

The pts vector contains the point coordinates where the bifurcation points are located. lbls give the labels that correspond to the pts vector and indicate what type of bifurcation point it represented. Finally a special object, here called biData contains the data together with a number of useful utilities. The import important of these is biData.plotBifurcationDiagram(pts, lbls) which takes pts and lbls as arguments. 

We can also print out a text summary of the computation using the command, auto.BifurcationSummary, which returns a summary of the findings.

1 2 3 4 5 6 7 8 9 10 11 12

Summary : BR PT TY LAB PAR(0) L2-NORM U(1) U(2) 1 1 EP 1 3.00000E+000 2.38908E+001 1.42336E+001 1.91879E+001 1 50 2 5.48632E-001 2.14502E+001 1.06592E+001 1.86144E+001 1 95 LP 3 -9.31897E-001 1.79088E+001 7.44448E+000 1.62881E+001 1 100 4 -9.13592E-001 1.74094E+001 7.22867E+000 1.58377E+001 1 150 5 1.74114E-001 1.26908E+001 6.15211E+000 1.10999E+001 1 200 6 1.54579E+000 8.36095E+000 5.44499E+000 6.34488E+000 1 234 LP 7 2.06965E+000 5.51178E+000 4.47540E+000 3.21723E+000 1 250 8 1.84381E+000 4.03137E+000 3.51306E+000 1.97747E+000 1 300 9 -5.77433E-001 7.02653E-001 -5.14930E-001 4.78089E-001 1 325 EP 10 -2.00625E+000 2.56175E+000 -2.55121E+000 2.32100E-001

We can manually plot the data by gaining access to the numpy version of the data. To do this we use:

1	pltData = biData.toNumpy

pltData is a numpy array where the first column is the bifurcation parameter and the remaning columns contain the species. For example to plot the bifurcation diagram for the first species in the model, R1 we would use:

1 2 3 4

plt.plot(x[:,0], x[:,1], linewidth=2) plt.axvline(0, color='black') plt.xlabel ("Signal") plt.ylabel ("R1")

I added a axvline command to draw a vertical line from the zero axis. I also added some axis labeling statements. These commands will result in:

What is interesting about this model is that the upper branch reaches the zero parameter value before the turning point. This means it is difficult to switch to the lower steady-state by just lowering the signal.

Viewing the Network

One other thing we can do is view the model as a network. Tellurium comes with a simple network viewer in the package nwed. import the viewer using

import nwed

at the ipython console. To view the network make sure the network viewer panel is visible, do this by going to the View menu, find panes and select, then look down the menu items, and near the bottom you'll find Network Viewer, select this option. To view the network, type the following at the ipython console.

nwed.setsbml (r.getSBML())

The viewer should now display something like:

Note that every view will be different and depends on the layout algorithm.