## A look at the Euler’s number: e

I’ve never particularly liked the way , Euler’s number, is introduced in textbooks. Most approaches give me a very limited intuitive feel for what 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 . These approaches include computing compound interest (I think this is one of the worst), integrating and noting that the area from 1 to the value of e is one (amazing but so what?), looking at the slope of at , the slope is one when (yes ok, and…) and finally the one I like is starting with and finding the power series solution and then using this to define .

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

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

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

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

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

We now differentiate this with respect to to give:

To find the function that also equals we need to discover the values for the coefficients, , etc. To make it easier, let’s decide that when , the value of is one, i.e: . This will mean that . What we’ll now do is match up the pairs of terms in and , that is set: and so on.

This allows us to state that: , , , and so on. Working back by subsituting into the equation and into the equation, and so on leads to the result that:

We therefore conclude that the equation, that satisfies is:

You can easily check by differentiating , you’ll get again. The solution to must therefore also be:

Not to labor the point too much, but differentiate and you’ll get the equivalence .

Let us define the value of the series when to be:

For convenience we will call this value :

The next question is, if then what does equal?

I’m going to make a jump here and propose that . From first principles, we can work out the derivative of and remarkably it is ! 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 . The Maclaurin series is a Taylor series centered on zero, that is:

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

But this is just the solution to , — see equation (1) — therefore we conclude that:

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

In other words the derivative of is the special case where . One reason why 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 ? What about its rate of increase? We know that the rate of increase is itself, , but how fast is that? One way to look at this is to derive the relative increase in . The relative growth rate is defined by:

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

If we now do the same for we find:

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

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

## 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 (pathwaydesigner.org) 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:

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

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

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:

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.

This code draw the following figure:

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

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]}-

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

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

## 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
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

// 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()