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:
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.
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:
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.
We can manually plot the data by gaining access to the numpy version of the data. To do this we use:
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:
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:
No comments:
Post a Comment