Friday, October 28, 2016

Bifurcation Analysis with Tellurium

Originally Posted on  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:

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:
bifirreversible2


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.




Saturday, October 15, 2016

Tikz Code for Drawing Metabolic Feedback Loops

Originally Posted on  by hsauro

I needed some figures that displayed a variety of different negative feedback loops so I created these using Tikz. Nothing particularly special. There are some absolute distances in the code which perhaps could be removed to make it more generic.



 

\documentclass{article}

\usepackage{amsmath}
\usepackage{tikz}
\usetikzlibrary{arrows}
\usetikzlibrary{calc}

\begin{document}

\begin{tikzpicture}[>=latex', node distance=2cm]
\node (Xo) {};
\node [right of = Xo] (S1) {\Large $x$};
\node [right of = S1] (S2) {};

\draw [->,ultra thick,blue] (Xo) -- node[above, black] {$v_1$} (S1);
\draw [->,ultra thick,blue] (S1) -- node[above, black] {$v_2$} (S2);

% Lets draw a line with a blunt end, -|
\draw [-|,ultra thick,blue]
% start in the middle of S2, and move down 2.75 mm
% ($ ... $) notation is used to add the coordinates
($ (S1) + (0mm,-2.75mm) $)

% Now draw the line down by 3mm
% -- means draw to, + means move by
-- +(0,-3mm)
% Now move back to the left of S2
% The symbol -| means draw horizontal then vertical.
% If we used -- instead the line would be drawn
% diagonally to the reaction edge.
% (S1) is the center of the node. But we want the
% blunt end to end below the S1 line and
% half way to the left. The 10mm is half the node
% distance of 2cm, and 1mm is slightly below
% the reaction line.
-| ($ (S1) - (10mm,1mm) $);
\end{tikzpicture}

\vspace{1cm}

\begin{tikzpicture}[>=latex', node distance=2cm]
\node (Xo) {};
\node [right of = Xo] (x1) {\Large $x_1$};
\node [right of = x1] (x2) {\Large $x_2$};
\node [right of = x2] (x3) {};

\draw [->,ultra thick,blue] (Xo) -- node[above, black] {$v_1$} (x1);
\draw [->,ultra thick,blue] (x1) -- node[above, black] {$v_2$} (x2);
\draw [->,ultra thick,blue] (x2) -- node[above, black] {$v_3$} (x3);

% Lets draw a line with a blunt end, -|, using the following coords
\draw [-|,ultra thick,blue]
% start in the middle of x2, and move down 2.75 mm
% ($ ... $) notation is used to add the coordinates
($ (x2) + (0mm,-2.75mm) $)

% Now draw the line down by an additional 3mm
% -- means draw to, + means move by
-- +(0,-3mm)

% Now move back to the left of x2
% The symbol -| means draw horizontal then vertical.
% If we used -- instead the line would be drawn
% diagonally to the reaction edge.
% (S1) is the center of the node. But we want the
% blunt end to end below the S1 line and
% half way to the left. The 10mm is half the node
% distance of 2cm, and 1mm is slightly below
% the reaction line.
-| ($ (x1) - (10mm,1mm) $);
\end{tikzpicture}

\vspace{1cm}

\begin{tikzpicture}[>=latex', node distance=2cm]
\node (Xo) {};
\node [right of = Xo] (x1) {\Large $x_1$};
\node [right of = x1] (x2) {\Large $x_2$};
\node [right of = x2] (x3) {\Large $x_3$};
\node [right of = x3] (x4) {};

\draw [->,ultra thick,blue] (Xo) -- node[above, black] {$v_1$} (x1);
\draw [->,ultra thick,blue] (x1) -- node[above, black] {$v_2$} (x2);
\draw [->,ultra thick,blue] (x2) -- node[above, black] {$v_3$} (x3);
\draw [->,ultra thick,blue] (x3) -- node[above, black] {$v_4$} (x4);

% Lets draw a line with a blunt end, -|
\draw [-|,ultra thick,blue]
($ (x3) + (0mm,-2.75mm) $)
-- +(0,-3mm)
-| ($ (x1) - (10mm,1mm) $);
\end{tikzpicture}

\vspace{1cm}

\begin{tikzpicture}[>=latex', node distance=2cm]
\node (Xo) {};
\node [right of = Xo] (x1) {\Large $x_1$};
\node [right of = x1] (x2) {\Large $x_2$};
\node [right of = x2] (x3) {\Large $x_3$};
\node [right of = x3] (x4) {\Large $x_4$};
\node [right of = x4] (x5) {};

\draw [->,ultra thick,blue] (Xo) -- node[above, black] {$v_1$} (x1);
\draw [->,ultra thick,blue] (x1) -- node[above, black] {$v_2$} (x2);
\draw [->,ultra thick,blue] (x2) -- node[above, black] {$v_3$} (x3);
\draw [->,ultra thick,blue] (x3) -- node[above, black] {$v_4$} (x4);
\draw [->,ultra thick,blue] (x4) -- node[above, black] {$v_5$} (x5);

% Lets draw a line with a blunt end, -|
\draw [-|,ultra thick,blue]
($ (x4) + (0mm,-2.75mm) $)
-- +(0,-3mm)
-| ($ (x1) - (10mm,1mm) $);
\end{tikzpicture}

\vspace{1cm}

\begin{tikzpicture}[>=latex', node distance=2cm]
\node (Xo) {};
\node [right of = Xo] (x1) {\Large $x_1$};
\node [right of = x1] (x2) {\Large $x_2$};
\node [right of = x2] (x3) {\Large $x_3$};
\node [right of = x3] (x4) {};

\draw [->,ultra thick,blue] (Xo) -- node[above, black] {$v_1$} (x1);
\draw [->,ultra thick,blue] (x1) -- node[above, black] {$v_2$} (x2);
\draw [->,ultra thick,blue] (x2) -- node[above, black] {$v_3$} (x3);
\draw [->,ultra thick,blue] (x3) -- node[above, black] {$v_4$} (x4);

% Lets draw a line with a blunt end, -|
\draw [-|,ultra thick,blue] (10mm, -8mm)
-- +(0,6mm);
\node (x) at (10mm,-11mm) {\Large $x$};
\end{tikzpicture}

\end{document} 

Tuesday, September 6, 2016

The Confusion of Modern Textbooks - Stylistic Sugar

August 6, 2016 11:43 am

I've been looking for a textbook on statistics for a class I'll teach in the autumn term. While the content of many textbooks might be ok the way the information is presented makes them difficult to read - at least it does for me. This applies to most undergraduate textbooks published today whether they be about statistics or other topics. There seems to be a need by publishers to embellish textbooks with so much stylistic sugar that the content is buried.

I scanned two typical pages from a second-hand stats textbook I bought from Goodwill to illustrate what I mean.



I've marked using a red star the many different styles used in the textbook on two random pages, these include:

1. Section heading with a vertical dark purple line
2. Highlighted area using three colors (black, blue and purple)
3. A notes section with blue heading
4. An Illustration section using blue font but spaced out letters and thin left-bar in dark purple
4. Up and down purple arrows in the illustration section indicating question and answer
6. Paragraph of text using back font (finally something normal)
7. Figure caption in dark purple text with green vertical line and horizonal line in black
8. Exercise box in the margin with heading in white with black background
9. Exercise box where the question in black with yellow background
10. Typewriter font for computer code with black and blue horizontal lines delimiting code

Not included in these pages are also other stylistic sugar:

1. Case study section that uses five colors and six stylistic features
2. Exercises with six stylistic features including at least eleven different symbols
3. Call-out in exercises using a script font with blue background.
3. Chapter practice tests in black font with blue background and thick blue horizontal line
4. Chapter Objective in black font, red bullet point in off yellow background with a vertical dotted line.
5. Chapter heading page, eight stylistic features with multiple fonts

And this is before a student has even started to read the content I counted at least 15 different fonts used in the text.

 

Tips for Matplotlib Users

Matplotlib is the plotting library that is commonly used with Python. However the documentation for matplotlib, though extensive is quite difficult to understand and navigate. Many of the most common things one might do when plotting are either not well explained or buried deep in the documentation. Here are some tips for the common things that I tend to do. More will be added as time goes by.

1. Saving a matplotlib plot as a PDF file:

import pylab as plt
a = [1,2,3,4]
b = [2,4,6,8]
plt.plot (a, b) 
plt.savefig ('myplot.pdf')

2. Setting the x and y axes labels and the main title:

import pylab as plt
a = [1,2,3,4]; b = [2,4,6,8] 
plt.plot (a, b)
plt.xlabel ('Time', fontsize=16)
plt.ylabel ('Variable')
plt.title ('My Experiment')
plt.show()

3. Set x and y axes limits

import pylab as plt
a = [1,2,3,4]; b = [2,4,6,8] 
plt.plot (a, b)
plt.xlim ((0, 2))
plt.ylim ((-5, 5))
plt.show()

4. Set the physical size of the plot to 8 inches by 6 inches

import pylab as plt
a = [1,2,3,4]; b = [2,4,6,8] 
plt.figure(figsize=(8,6))
plt.plot (a, b)
plt.show()

5. Setting the color and width of plotted lines

import pylab as plt
a = [1,2,3,4]; b = [2,4,6,8] 
plt.figure(figsize=(8,6))
plt.plot(a, b, color="blue", linewidth=2.5)
plt.show()

7. Set the x and y axis labels and graph title

import pylab as plt
a = [1,2,3,4]; b = [2,4,6,8] 
plt.plot(a, b, color="blue", linewidth=2.5)
plt.show()