## Plotting Bar graph of Species Concentations in Tellurium

I had a model with 27 flaoting species and I wanted to plot the steady state concentrations on a histogram where the labels were the names of the different species. Here is a general purpose script that will do that:

def plotFloatingSpecies (r, width=12, height=6):
import matplotlib.pyplot as plt
xlabels = r.getFloatingSpeciesIds()
concs = r.getFloatingSpeciesConcentrations()
plt.figure(figsize=(width,height))
plt.bar(xlabels, concs, label=xlabels)
plt.xticks(range (len (xlabels)), xlabels, ha='right', rotation=45)


The function has a default width and height for the resulting plot. The default widens the usual size so that the labels don’t collide. The first argument is the libroadrunner instance.  To given a contrived example, the following is a model of 8 floating species that form a linear chain governed by simple irreversible mass-action kinetics.

import tellurium as te

def plotFloatingSpecies (r, width=12, height=6):
import matplotlib.pyplot as plt
xlabels = r.getFloatingSpeciesIds()
concs = r.getFloatingSpeciesConcentrations()
plt.figure(figsize=(width,height))
plt.bar(xlabels, concs, label=xlabels)
plt.xticks(range (len (xlabels)), xlabels,  ha='right', rotation=45)

$Xo -> S1; k1*Xo; S1 -> S2; k2*S1; S2 -> S3; k3*S2; S3 -> S4; k3*S3; S4 -> S5; k4*S4; S5 -> S6; k5*S5; S6 -> S7; k4*S6; S7 -> S8; k3*S7; S8 -> ; k4*S8; k1 = 0.3; k2 = 0.5; k3 = 0.27; k4 = 0.9; k5 = 0.14 Xo = 10; """) r.steadyState() plotFloatingSpecies (r, width=6, height=3)  The resulting plot is shown below: Posted in Modeling, Pathways, Programming, Python, Software, Tellurium | Leave a comment ## Multilayered Cascade using TikZ A summer student working in my lab, Ming Hong Lui from Hong Kong University (HKUST), worked on the perturbation analysis of signaling cascades and in his writeup he use TikZ to draw a nice cascade diagram which I present here. \documentclass{article} \usepackage{tikz} \usetikzlibrary{arrows} \begin{document} \begin{tikzpicture}[>=latex',node distance = 2cm] \node (S1) {$S_1$}; \node [right of = S1] (S2) {$S_2$}; \node [right of = S1, node distance = 1.8cm] (X1) {}; \node [below of = X1] (S3) {$S_3$}; \node [right of = S3] (S4) {$S_4$}; \draw [->,very thick] (S1) to [bend left=40] node[above] {$v_1$} (S2); \draw [->,very thick] (S2) to [bend left=40] node[below] {$v_2$} (S1); \draw [->,very thick] (S3) to [bend left=40] node[above] {$v_3$} (S4); \draw [->,very thick] (S4) to [bend left=40] node[below] {$v_4$} (S3); \node [right of= S3, node distance = 1cm] (X2){}; \node [above of= X2, node distance = 1cm] (X3){}; \draw [->,very thick] (S2) to [bend left=20] (X3); \node [right of= S4, node distance = 0.8cm] (X4){}; \node [below of= X4, node distance = 1cm] (X5){}; \draw [->,very thick] (S4) to [bend left=20] (X5); \node [right of = S3, node distance = 1.8cm] (X3_5) {}; \node [below of = X3_5] (S5) {$S_5$}; \node [right of = S5] (S6) {$S_6$}; \node [right of = S5, node distance = 1.8cm] (X5) {}; \node [below of = X5] (S7) {$S_{2n-1}$}; \node [right of = S7] (S8) {$S_{2n}$}; \draw [->,very thick] (S5) to [bend left=40] node[above] {$v_5$} (S6); \draw [->,very thick] (S6) to [bend left=40] node[below] {$v_6$} (S5); \draw [->,very thick] (S7) to [bend left=40] node[above] {$v_{2n-1}$} (S8); \draw [->,very thick] (S8) to [bend left=40] node[below] {$v_{2n}$} (S7); \node [right of= S7, node distance = 1cm] (X6){}; \node [above of= X6, node distance = 1cm] (X7){}; %\node at ($(S6)!.5!(X7)$) {$\ddots$}; \draw [dashed, ->,very thick] (S6) to [bend left=20] (X7); \node [right of = S1, node distance = 1cm] (T1) {$T_1, J_1$}; \node [right of = S3, node distance = 1cm] (T2) {$T_2, J_2$}; \node [right of = S5, node distance = 1cm] (T3) {$T_3, J_3$}; \node [right of = S7, node distance = 1cm] (T4) {$T_n, J_n$}; \end{tikzpicture} \end{document}  This code draw the following figure: | Leave a comment ## Another Inhibition Pathway Diagram using TikZ Here is another pathway diagram I needed to draw using TikZ. IN this case I needed an inhibited step. This was more tricky because I needed the inhibition line to point midway to a reaction but without touching the reaction itsef. To solve this I had to ask stackoverflow and within 2 hours I had an answer. The linear portion of the pathway is straightforward:  \begin{tikzpicture}[>=latex',node distance = 2cm] \node (S) {\scalebox{1.4}{$S$}}; \node [left of = S] (Xo) {}; \draw[->,line width = 1.2pt] (Xo) to node[above] {\scalebox{1.3}{$v_1$}} (S); \node [right of = S] (P) {$P$}; \draw[->,line width = 1.2pt] (S) to node[below] {\scalebox{1.3}{$v_2$}} (P); \node [right of = P] (X1) {}; \draw[->,line width = 1.2pt] (P) to node[above] {\scalebox{1.3}{$v_3$}} (X1); \end{tikzpicture}  The question was how to draw the inhibition line from inhibitior I to reaction v2. The trick is to use the tikz calc extension package which can be used to do more more advanced coordinate calculations. \usetikzlibrary{calc} The line to add that that will draw the vertical inhibition line is: \draw[|-,line width=1.2pt] ([yshift=4pt]$(S)!.5!(P)$) –++(0,1.5cm)node[above]{$I$}; Let’s decompose this. The first two arguments in the \draw command [[|-,line width=1.2pt] simply set the arrow style (blunt end) and the line width. The second part specifies the coordinates of the line itself and the remainder node[above]{$I$}; indicates what text to draw and where to draw it relative to the line. The real meat is in the drawing line section, that is: ([yshift=4pt]$(S)!.5!(P)$) –++(0,1.5cm) The basic syntax for a line coordinate is (x1,y1) — (x2,y2). In the example there is a modification to this where the coordinate is specified as (x1,y1) — ++(x2,y2). The ++ means that the coordinate at x2,y2 is computed relative to x1,y1, that is the coordinate at x2,y2 is actually (x1,y1)+(x2,y2). The (0,1.5cm) then means that the coordinate for the end point has the same x coordinate but the y coordinate is displaced upwards by 1.5cm. Note that the tikz y coordinate is like a normal graph plotting axis. The trickest bit is computing the starting point for the vertical line. Note that this has to be mid point between nodes S and P. This requires some coordinate calculations. The x,y coordinate is specified by$(S)!.5!(P)$). The bit in front, ([yshift=4pt], just moves the computed y coordinates up 4pts. The text inside the$$means that we are doing a coordinate calculation,ie$….$represents a math calculation. (S) and (P) represent the coordinates of the nodes S and P respectively. The important bit is !.5!. The explanation point is a pathway modifier and can be put between two coordinates. (S)!.5!(P) means compute the coordinate that is half way between (S) and (P). The last thing to add is that I nocticed that the blunt end was a bit too sort. To widen the blunt end I used the arrow.meta extension package. This allows one to modify the size of the arrow heads, in this case I used the following to widen the blunt end to 5mm. {|[width=5mm]}- | Leave a comment ## Drawing a Pathway Fragment with Inhibition using Tikz Here is another quick pathway fragment I needed today. This won’t scale well because I’ve used some fixed dimensions eg the width of the lines and the size of the text. But these are easily adjusted if you want to size the figure differently. The node distance = 2.5cm gives the sizes for the arrows except for v3. One might be able to calculate the some of the fixed sizes based on the node distance but I didn’t have time to investigate this.   \documentclass{article} \usepackage{tikz} \usetikzlibrary{arrows} \begin{document} \begin{tikzpicture}[>=latex',node distance = 2.5cm] \node (S) {\scalebox{1.4}{$S$}}; \node [left of = S] (Xo) {}; \draw[->,line width = 1.2pt] (Xo) to node[above] {\scalebox{1.3}{$v_1$}} (S); \node [right of = S] (X1) {}; \draw[->,line width = 1.2pt] (S) to node[above] {\scalebox{1.3}{$v_2$}} (X1); \node[above of = S, node distance = 1.5cm] (V3) {}; \draw[-|,line width = 1.2pt] (S) to (V3) {}; \node[left of = V3, node distance = 1.25cm] (Xa) {}; \node[right of = V3, node distance = 1.25cm] (Xb) {}; \draw[->,line width = 1.2pt] (Xa) to node[above] {\scalebox{1.3}{$v_3$}} (Xb); \end{tikzpicture} \end{document}   | Leave a comment ## Drawing a ‘Futile’ Cycle using Tikz I was in need of a diagram of a futile cycle. I tried illustator but I didn’t have the right LaTeX fonts so the figure did blend well with the rest of the document. I decided to make one using TiKz. And here it is.  \documentclass{article} \usepackage{tikz} \usetikzlibrary{arrows} \begin{document} \begin{tikzpicture}[>=latex',node distance = 2.5cm] \node (S2) {$S_2$}; \node [left of = S2] (S1) {$S_1$}; \node [left of = S1, node distance = 1.5cm] (Xo) {}; \draw[->,thick] (Xo) to node[above] {$v_1$} (S1); \node [right of = S2, node distance = 1.5cm] (X1) {}; \draw[->,thick] (S2) to node[above] {$v_4$} (X1); \draw [->,thick] (S2) to [bend left=40] node[below] {$v_3$} (S1); \draw [->,thick] (S1) to [bend left=40] node[above] {$v_3$} (S2); \end{tikzpicture} \end{document}  Posted in General Interest, LaTeX, Pathways, Publishing, Textbooks | Leave a comment ## Plotting 3D graphs using Python and Tellurium As an example I wanted to show how one could plot a 3D phase plot. A great example to use for this is the Lorenz Attractor. This system is interesting because it displays chaotic behavior. The differential equations for the system are given by the following three: Different values for the parameters, sigma, rho and beta, lead to different behaviors. The Python script below that uses Tellurium plots a 3D phase plot: # The Lorenz Model # The Lorenz system of equations is the classic set of # equations which exhibit chaotic dynamics. The equations # represent a highly simplified model of 2D thermal # convection known as Rayleigh-Benard convection; # The three variables, x, y and z represent the # convective overturning, the horizontal temperature # and vertical temperature variation respectively; ## The parameters a, b and c represent approximately # the energy losses within the fluid due to viscosity, # the temperature difference between the two plates # which separate the fluid and the ratio of the vertical # height of the fluid layer to the horizontal extent # of the convective rolls respectively; import tellurium as te import roadrunner from mpl_toolkits.mplot3d import Axes3D import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np r = te.loada (''' // Define the model in terms of three differential equations;$s -> x; -a*(x - y);
$s -> y; -x*z + b*x - y;$s -> z; x*y - c*z;

a = 10; b = 28; c = 2.67;
''')

# Uncomment/comment the following lines to obtain behaviour you wish to observe;

# Non-chaotic oscillations;
#r.x = 5.43027; r.y = 26.74167; r.z = 46.0098;
#r.a = 4; r.b = 91.9; r.c = 2.67;

# Single oscillation;
#r.x = 32.82341; r.y = -44.84651; r.z = 258.02389;
#r.a = 4; r.b = 191.9; r.x = 2.67;

#r.x = 0.1; r.y = 0.1; r.z = 0.1;
#r.a = 10; r.b = 10; r.c = 2.67;

# Chaos
r.x = 0.96259; r.y = 2.07272; r.z = 18.65888;

m = r.simulate (0, 40, 3000, ['x', 'y', 'z']);

fig = plt.figure(figsize=(8,8))
ax = fig.gca(projection='3d')

#theta = np.linspace(-4 * np.pi, 4 * np.pi, 100)
ax.plot(m[:,0], m[:,1], m[:,2], label='Lorenz Attractor')
ax.legend()

plt.show()


## How do I change the simulation tolerances in Tellurium?

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()
'this' : 22F59350
'modelName' : __main
'libSBMLVersion' : LibSBML Version: 5.13.0
'jacobianStepSize' : 1e-005
'conservedMoietyAnalysis' : false
'simulateOptions' :
{
'this' : 1B89AF28,
'reset' : 0,
'structuredResult' : 0,
'copyResult' : 1,
'steps' : 50,
'start' : 0,
'duration' : 20
},
'integrator' :
name: cvode
settings:
relative_tolerance: 0.00001
absolute_tolerance: 0.0000000001
stiff: true
maximum_bdf_order: 5
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 maximium number of steps and the initial time step.

The smaller the relative tolerance the more accuate the solution, however too small a value will result in either excessive runtimes or more likely rounadoff errors. A relative tolerance of 1E-4 means that errors are controlled to 0.01%. An optimal value is roughly 1E-6. The absoute tolerance is used when a variable gets so small that the relative tolerance don’t make much sense to apply. In these situations the 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 tolerances.

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 if refered 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, this 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 beteen 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 parmaeters of the model. For example how does the model bahave 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) plt.figure(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], "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 Posted in Modeling, Pathways, Programming, Software, Systems Theory | Leave a comment ## 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 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:




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
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.

## How to do a simple simulation using Tellurium

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

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

import tellurium as te

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: