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 tool bar 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 formating of the output using savetxt, for example the following will output the individual stoichiometry coefficeint using 3 decimal places, 5 characters minimium, 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.