Simple Stochastic Code in Python

It’s been a while since I blogged (grant writing etc gets int he way of actual thinking and doing) so here is a quick post that uses Python to model a simple reaction A -> B using the Gillespie next reaction method.

I know importing pylab is discouraged but I can never remember how to import matplotlib. pylab is an easy way to get matplotlib. It also imports scipy and numpy but I import numpy separately.

I use lists instead of array to accumulate the results because they are faster to extend as the data comes in. Note that with a Gillespie algorithm you don’t know beforehand how many data points you’ll get.

I’ll also make a shameless plug to my book Enzyme Kinetics for Systems Biology which explains the algorithm.

Posted in General Science Interest, Modeling, Programming, Python, Software | Leave a comment

New Textbook on Metabolic Control Analysis

New Textbook Published

I am pleased to announce the publication of my new textbook:

Introduction to Metabolic Control Analysis


“This book is an introduction to control in biochemical pathways. It introduces students to some of the most important concepts in modern metabolic control. It covers the basics of metabolic control analysis that helps us think about how biochemical networks operate. The book should be suitable for undergraduates in their early to mid years at college.”


Available at Amazon for $49.95 or directly from me for only $29.95 at


The book is printed in full color with 275 pages, 118 Illustrations, 71 Exercises.


1, Traditional Concepts in Metabolic Regulation
2. Elasticities
3. Introduction to Biochemical Control
4. Linking the Parts to the Whole
5. Experimental Methods
6. Linear Pathways
7. Branched and Cyclic Systems
8. Negative Feedback
9. Stability
10. Stability of Negative Feedback Systems
11. Moiety Conserved Cycles
12. Moiety Conserved Cycles

Appendix A: List of Symbols and Abbreviations
Appendix B: Control Equations

Posted in Metabolic Control Analysis, Modeling, Pathways, Publishing, Systems Theory, Textbooks | Leave a comment

Repeatable, Reproducible [and Replictable]

There appears to be great confusion in the scientific and social sciences communities on the meaning of words related to certain aspects of the scientific method. The ArXiv paper by Lorena Barba “Terminologies for Reproducible Research” highlights the confused state that has appeared over the last 20 years. The words in question include:


I will dispense with replication simply because it’s too hard to say quickly (especially replicability) but see below for a more serious reason. The contention appears in the meaning of reproducibility or to reproduce. As Barba points out there are at least three ‘camps’ in this community which she labels, A, B1, and B2. The A camp makes no distinction, so we’ll forget about those. To describe B1 and B2 we must look at two extreme scenarios with respect to an experiment:

1. An experiment is carried out and is done again by the same author, using the same equipment, same methods, basically the same everything.

2. The experiment is carried out by a third-party using different equipment, different methods, etc. Basically, everything is different

In between these two extremes are variants, For example, the third-party could use the same methods but implement them independently of the original author, usually by reading the description given in the original published paper.

Given these descriptions, the B1 group calls the first scenario, ‘to reproduce the experiment’ while the second group, B2, calls the first scenario, ‘to replicate the experiment’ and there lies the contention.

Personally, I don’t like either of these terms as used here. As I mentioned before, replicability is a hard word to say. But not only that, from a dictionary perspective, it means the same thing as reproducibility. The Oxford English Dictionary describes replicability as “The quality of being able to be exactly copied or reproduced.” So why use two words, for two quite different things, where the two words have essentially the meaning?

My personal choice are the following two words, Rather than use the word replicability which seems redundant, I choose repeatability, hence:

Repeatability: means ‘to repeat the experiment again’, the word implies that the experiment was done exactly as before – Scenario 1

Reproducibility: means: ‘to recreate the experiment anew; reproduce implies creating a new thing, independently of the old – Scenario 2.

Of course one can get much more fine-grained, especially when it comes to computational experiments. But the fine graining can be included as levels within the class reproducibility.

Other than a change in wording from replicability to reproducibility, I appear to belong to camp B2. I should list others in camp B2, these include FASEB, NIST, 6 sigma, ACM, and Wikipedia and this Wikpedia page and The Physiome Project. I am sure there are others. For example, ACM writes:

Repeatability (Same team, same experimental setup)

The measurement can be obtained with stated precision by the same team using the same measurement procedure, the same measuring system, under the same operating conditions, in the same location on multiple trials. For computational experiments, this means that a researcher can reliably repeat her own computation.

Reproducibility (Different team, different experimental setup)

The measurement can be obtained with stated precision by a different team, a different measuring system, in a different location on multiple trials. For computational experiments, this means that an independent group can obtain the same result using artifacts which they develop completely independently.

Essentially the same definitions I gave above.

The National Academies recently studied this issue closely and is soon coming out with a report. It is apparently in favor of the more confusing option.

Posted in General Science Interest, Publishing | Leave a comment

How to do a simple parameter scan using Tellurium

A common task in modeling is to see how a parameter influences a model’s dynamics. For example, consider a simple two reaction pathway:

-> S1 ->

where the first reaction has a fixed input of vo and the second reaction a first-order rate laws k1*S1. The task is to investigate how the time course of S1 is influenced by vo.

The script below defines the model, then changes vo in increments and plots the effect on the pathway via a time course simulation. To do the parameter scan we exploit plotArray. This initially prevents the plots from being shown using show=False. We make sure that each plot gets a different color using resetColorCycle=False, finally, we show the plot using show(). To make things more interesting we also add a legend entry for each plot.

Note we call reset each time we run a simulation to ensure that S1 is reset back to its initial condition.

The following figure shows the resulting plot:

Thanks to Kiri Choi for pointing out how to use plotArray in this way.

Posted in General Science Interest, Modeling, Pathways, Programming, Python, Software, Tellurium | Leave a comment

How to plot a grid of phase plots using Tellurium

Let’s say we have a chemical network model where the species oscillate and you’d like to plot every combination of these on a grid. If so, then this code might be of help.

This code will generate:

Here is an alternative that has removed the internal ticks and labels and arranges the column and row labels as contiguous. This is probably more what one might expect. This version plots all combinations including the transpose combinations.

The following shows the steady state oscillations. This was done by simulating the model twice and using the plotting the second set of simulation results.

Posted in Enzyme Kinetics, Modeling, Pathways, Programming, Python, Software, Tellurium | Leave a comment

A look at the Euler’s number: e

I’ve never particularly liked the way e, Euler’s number, is introduced in textbooks. Most approaches give me a very limited intuitive feel for what e actually is. Modern textbooks appear to use one of four common ways to introduce e, and only one of them gives me a semblance of an intuitive feel for e. These approaches include computing compound interest (I think this is one of the worst), integrating 1/x and noting that the area from 1 to the value of e is one (amazing but so what?), looking at the slope of a^x at x = 0, the slope is one when a = e (yes ok, and…) and finally the one I like is starting with dy/dt = y and finding the power series solution and then using this to define e.

Before continuing let’s state what e is numerically equal to:

e = 2.71828….

For normal algebra, all we need are addition, and multiplication (subtraction and division are just alternative forms of these). For convenience, we also introduce a power notation such as x^2 and x^n, but these are just a short-hand notation for doing lots of multiplications at once. With the introduction of trigonometry which involves relationships between sides and angles of triangles, the basic algebra becomes cumbersome because it involves the use of infinite series. Rather than writing the series down all the time, we define short-cut names such as sine, cosine etc. Calculus brings us another special type of series which involves the solution to differential equations. The short-cut notation for this series is e^x.

Note that in all these cases we are still only doing addition and multiplication. The functions sine and e^x are just short-hand for particular combinations of addition and multiplication that happen to be in the form of an infinite series.

Given that e pops up in the form of e^x when solving differential equations, I think this is the place to start. Let’s consider the simplest possible non-trivial differential equation:

    \[ \frac{dy}{dx} = y \]

This equation is saying something quite interesting, that the derivative of y is the same as y. What this means is that if we were to find a solution to this differential equation, the solution would have to also equal dy/dt. We can express the above equation in the form:

    \[ f'(x) = f(x) \]

this perhaps makes it more obvious that the derivative is the same as the function f(x). Can we find an actual function that is like this? The way to find this answer is to define f(x) as a general power series:

    \[ f(x) = a_o + a_1 x + a_2 x^2 + a_3 x^3 + \ldots + a_n x^n \]

We now differentiate this with respect to x to give:

    \[ f'(x) = a_1 + 2 a_2 x + 3 a_3 x^2 + \ldots + (n + 1) a_n x^n \]

To find the function f(x) that also equals f'(x) we need to discover the values for the coefficients, a_o, a_1, etc. To make it easier, let’s decide that when x=0, the value of f(0) is one, i.e: f(0) = 1. This will mean that a_o = 1. What we’ll now do is match up the pairs of terms in f(x) and f'(x), that is set: a_o = a_1, a_1 x = 2 a_2 x, a_2 x^2 = 3 a_3 x^2 and so on.

This allows us to state that: a_1 = 1, a_2 = 1/2, a_3 = a_2/3, a_4 = a_3/4 and so on. Working back by subsituting a_3 into the a_4 equation and a_2 into the a_3 equation, and so on leads to the result that:

    \[ a_n = \frac{1}{n!} \]

We therefore conclude that the equation, f(x) that satisfies f(x) = f'(x) is:

    \[ f(x) = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} \hdots \]

You can easily check by differentiating f(x), you’ll get f(x) again. The solution to dy/dt = y must therefore also be:

    \[ y = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} \hdots \qquad\qquad (1) \]

Not to labor the point too much, but differentiate y and you’ll get the equivalence dy/dx = y.

Let us define the value of the series when x = 1 to be:

    \[ f(1) = 1 + 1 + \frac{1}{2!} + \frac{1}{3!} + \hdots + \frac{1}{n!} \]

For convenience we will call this value e:

    \[ e = f(1) = 1 + 1 + \frac{1}{2!} + \frac{1}{3!} + \hdots + \frac{1}{n!} \]

The next question is, if f(1) = e then what does f(x) equal?

I’m going to make a jump here and propose that f(x) = e^x. From first principles, we can work out the derivative of e^x and remarkably it is e^x! You’ll find a proof at the excellent Paul’s Online Math Notes

Let’s now use the Maclaurin series to find an approximation to e^x. The Maclaurin series is a Taylor series centered on zero, that is:

    \[ f(x)=f(0)+xf'(0)+\frac{x^2}{2 !}f''(0)+...+\frac{x^n}{n !}f^{(n)}(0)+...\]

Since the derivative of e^x is e^x, as well as the third, fourth fifth derivatives etc and noting that e^x at x = 0 is one then we can write:

    \[ f(x) = 1 + 1 \cdot x + \frac{1 \cdot x^2}{2} + \frac{1 \cdot x^3}{6} + \ldots = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \ldots = \sum_{n=0}^{\infty} \frac{x^n}{n~} \]

But this is just the solution to dy/dx= y, — see equation (1) — therefore we conclude that:

    \[ y = e^x \]

It is also worth noting that the derivative for the general exponential a^x is given by:

    \[ \frac{d a^x}{dx} = a^x \log (a) \]

In other words the derivative of e^x is the special case where \log (a) = 1. One reason why e^x is special is that its the purest exponential in the sense that the derivative has no scaling factor. It’s the canonical exponential, the simplest possible. It’s the only exponential function where the function and its derivative are identical.

What else can we say about e^x? What about its rate of increase? We know that the rate of increase is itself, e^x, but how fast is that? One way to look at this is to derive the relative increase in e^x. The relative growth rate is defined by:

    \[ (dy/y)/(dx/y) = \frac{dy}{dx} \frac{x}{y} \]

Given this definition, let’s look first at how fast a^x is increasing by computing

    \[ \frac{d a^x}{dx} \frac{x}{a^x} = a^x \log (x) \frac{x}{a^x} = \log (a) x \]

If we now do the same for e^x we find:

    \[ \frac{d e^x}{dx} \frac{x}{e^x} = e^x \log (e) \frac{x}{e^x} = x \]

In other words e^x is the only exponential function that increases at a rate equal to x. For example if x = 1, it means that a 1% increase in x leads to a 1% increase in e^x while if x = 10 then a 1% increase in x leads to a 10% increase in e^x. This is the characteristic of exponential functions, they increase at an ever increasing rate. Contrast this with a power term such as x^2. If we compute the relative increase for x^2 we find it has a fixed increase of 2:

    \[ \frac{d x^2}{dx} \frac{x}{x^2} = 2 x\frac{x}{x^2} = 2 \]

The growth of a bacterial colony follows this pattern, where the growth in the colony is a fixed percentage over time.

Posted in General Interest, Math | Leave a comment

Smallest Chemical Reactions Systems that is Bistable

A while back Thomas Wilhelm, published a paper that described the smallest chemical network that could display bistability. The paper that describes this result is:

Wilhelm, T. (2009). The smallest chemical reaction system with bistability. BMC systems biology, 3(1), 90.

This is a diagram of the network generated using pathwayDesigner:

Here is a Tellurium script that uses Antimony to define the model (Note that $P means that species P is fixed). The S term in the first reaction is supposed to represent an input signal.

Using the auto200 extension to roadrunner we can plot the bifurcation diagram for this system as a function of the signal S. If S is below 0.8 only one stable steady state exists and the values for X and Y are both zero. Above S = 0.8 we see three steady states emerge.

For the system where it shows three steady states, one is at zero concentration, and is stable but not shown in the plot. The other two are marked by the horizontal line roughly at 1.5 on the y-axis and is unstable. The other is represented by the line that moves up from the turning point at about x = 0.8. This steady state is stable. The unstable branch appears to asymptotically approach a limiting value at high S values, 1.5 for X and approximately 0.00028 for Y.

The paper also describes what happens when we add a fixed input flux to X at a rate of 0.6. This can be simply done by adding the line J5: ->X; 0.6; to the model. This change results in a more classical look for the bifurcation plot as shown below:

The Tellurium script for generating these bifurcation plots is shown below:

Posted in Modeling, Pathways, Python, SBML, Software, Systems Theory, Tellurium | Leave a comment

Smallest Chemical Reaction System with Hopf Bifurcation

A while back Wilhelm, and Heinrich published a paper that described the smallest chemical network that could display a Hopf bifurcation. That is, the chemical species oscillated. The paper that describes this result is:

Wilhelm, Thomas, and Reinhart Heinrich. “Smallest chemical reaction system with Hopf bifurcation.” Journal of mathematical chemistry 17.1 (1995): 1-14.

This is a diagram of the network taken from their paper:

Here is a Tellurium script that uses Antimony to define the model (Note that $A means that species A is fixed):

The paper gives parameter values that result in oscillations but I wanted to find some other parameter sets. One way to do this is to load the model into the SBW slider control. By changing the sliders one can observe the dynamics. Here I used pathwaydesigner ( to do the same thing. I first exported the model as SBML using:

I loaded the saved SBML file into pathwayDesigner and used the autolayout plugin to get the following network:

I started the slider plugin and varied the sliders until I saw oscillations. In fact, it didn’t take much to get oscillations, all I had to do was increase the value of k1:

I copied the new value of k1 to the Tellurium model (see above) and got the following output:

Posted in Modeling, Pathways, Python, SBML, Software, Tellurium | Leave a comment

C Based Reduce Row Echelon Code

I recently needed some code to compute the reduced row echelon of a matrix. Applications such as Matlab, Mathematics, sympy and R support this functionality out of the box. Libraries such as LAPACK do not, including the linear algebra package in Python. The linear algebra package in Python has some surprising omissions which appear to be by design.

Here I include a C based function that will compute the reduced row echelon. It uses partial pivoting to prevent numerical stability but I don’t know how it would fare when confronted with a large matrix (i.e > 500 rows or columns).

It uses a simple matrix type called TMATRIX. It shouldn’t be difficult to modify code if you use a different matrix structure.

To call rref use:

The 1e-9 argument is the setting that decides whether a value is zero or not. Numbers below tolerance are considered zeros.

The output matrix givne the input should be:

Posted in Modeling, Pathways, Programming, Software | Leave a comment

Plotting Bar graph of Species Concentrations 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:

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.

The resulting plot is shown below:


Posted in Modeling, Pathways, Programming, Python, Software, Tellurium | Leave a comment