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 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 includng oscill8), AUTO2000 in Telluirum does not require a separate compiler to compile the model for AUTO2000 to compute the differential equations. This makes it easier to delpoy 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 text book “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 byTyson (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:
import tellurium as te
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, teplugins to access AUTO2000 and pylab to gain acess to matplotlib. Once we have the model loaded we can get a handle on AUTO2000 by calling teplugins.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 propoerties 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.
auto = teplugins.Plugin("tel_auto2000")
To run the bifurcation analysis we use the Python code:
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:
pts = auto.BifurcationPoints
lbls = auto.BifurcationLabels
biData = auto.BifurcationData
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 utilties. The import important of these is biData.plotBifurcationDiagram(pts, lbls) which takes pts and lbls as arguments. Running this code will generate the bifurcation plot shown below.
We can also print out a text summary of the computation using the command, auto.BifurcationSummary, which returns a summary of tjhe fidnnngs.
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:
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:
plt.plot(x[:,0], x[:,1], linewidth=2)
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 thee Network
One other things we can do is view the model as a network. Telluirum comes witha simple network viewer in the package nwed. import the viewer using
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 trhe ipython console.
The viewer should now display something like:
Note that every view will be different and depends on the layout algorithm.