Bifurcation Analysis with Tellurium

Originally Posted on October 28, 2016 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:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

import tellurium as te import rrplugins import pylab as plt   r = te.loada ('''     $X ->; R1;  k1*EP + k2*Signal;     R1 ->; $w;  k3*R1;     EP ->; E;   Vm1*EP/(Km + EP);     E -> EP;   ((Vm2+R1)*E)/(Km + E);     Vm1 = 12; Vm2 = 6;     Km = 0.6;     k1 = 1.6; k2 = 4;     E = 5; EP = 15;     k3 = 3; Signal = 0.1; ''')

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.

1 2 3 4 5 6 7 8 9

auto = rrplugins.Plugin("tel_auto2000")   auto.setProperty("SBML", r.getCurrentSBML()) auto.setProperty("ScanDirection", "Negative") auto.setProperty("PrincipalContinuationParameter", "Signal") auto.setProperty("PreSimulation", "True") auto.setProperty("PreSimulationDuration", 1.0) auto.setProperty("RL0", -2.0) auto.setProperty("RL1", 3.0)

To run the bifurcation analysis we use the Python code:

1	auto.execute()

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:

1 2 3 4

pts     = auto.BifurcationPoints lbls    = auto.BifurcationLabels biData  = auto.BifurcationData biData.plotBifurcationDiagram(pts, lbls)

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.

1 2 3 4 5 6 7 8 9 10 11 12

Summary : BR PT TY LAB PAR(0) L2-NORM U(1) U(2) 1 1 EP 1 3.00000E+000 2.38908E+001 1.42336E+001 1.91879E+001 1 50 2 5.48632E-001 2.14502E+001 1.06592E+001 1.86144E+001 1 95 LP 3 -9.31897E-001 1.79088E+001 7.44448E+000 1.62881E+001 1 100 4 -9.13592E-001 1.74094E+001 7.22867E+000 1.58377E+001 1 150 5 1.74114E-001 1.26908E+001 6.15211E+000 1.10999E+001 1 200 6 1.54579E+000 8.36095E+000 5.44499E+000 6.34488E+000 1 234 LP 7 2.06965E+000 5.51178E+000 4.47540E+000 3.21723E+000 1 250 8 1.84381E+000 4.03137E+000 3.51306E+000 1.97747E+000 1 300 9 -5.77433E-001 7.02653E-001 -5.14930E-001 4.78089E-001 1 325 EP 10 -2.00625E+000 2.56175E+000 -2.55121E+000 2.32100E-001

We can manually plot the data by gaining access to the numpy version of the data. To do this we use:

1	pltData = biData.toNumpy

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:

1 2 3 4

plt.plot(x[:,0], x[:,1], linewidth=2) plt.axvline(0, color='black') plt.xlabel ("Signal") plt.ylabel ("R1")

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:

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.

Tikz Code for Drawing Metabolic Feedback Loops

Originally Posted on October 15, 2016 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} 

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

Bode Plots using Python

 I needed a quick way to plot some Bode plots for a second-order system. I didn't have access to Matlab, instead, I searched for a solution using Python, and I found one. Documentation is a bit sparse so this example might be helpful.

The signals package supports the signal.bode method which turned out to be quite easy to use. Signal is part of the scipy package and is something we bundle with our Tellurium platform.

There is a more comprehensive discussion of Python and Control Theory here.

[code lang="python" fontsize="10"]
#
import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

# Coefficients in numerator of transfer function
num = [1]
# Coefficients in denominator of transfer function
# High order to low order, eg 1*s^2 + 0.1*s + 1
den = [1, 0.1, 1]

# Scan over zeta, a parameter for a second-order system
zetaRange = np.arange(0.1,1.1,0.1)

f1 = plt.figure()
for i in range(0,9):
    den = [1, 2*zetaRange[i], 1]
    print den
    s1 = signal.lti(num, den)
    # Specify our own frequency range: np.arange(0.1, 5, 0.01)
    w, mag, phase = signal.bode(s1, np.arange(0.1, 5, 0.01).tolist())
    plt.semilogx (w, mag, color="blue", linewidth="1")
plt.xlabel ("Frequency")
plt.ylabel ("Magnitude")
plt.savefig ("c:\\mag.png", dpi=300, format="png")

plt.figure()

for i in range(0,9):
    den = [1, zetaRange[i], 1]
    s1 = signal.lti(num, den)
    w, mag, phase = signal.bode(s1, np.arange(0.1, 10, 0.02).tolist())
    plt.semilogx (w, phase, color="red", linewidth="1.1")
plt.xlabel ("Frequency")
plt.ylabel ("Amplitude")
plt.savefig ("c:\\phase.png", dpi=300, format="png")

[/code]

Output from Python script:

Magnitude plot:

Phase plot:

Simple RK4 Code

January 4, 2015 10:11 pm 

Thought I'd start off the New Year with a simple freebie, a 4th Order Runge-Kutta class for Delphi, Windows and Mac OS. Should also work on Android and iOS mobile platforms and with Free Pascal. The code below shows how you'd use it. First declare the set of differential equations that need to be solved (TDoubleDynArray type can be found in Types):

 TMySystem = class (TObject)
    class procedure func (time : double; y, p, dydt : TDoubleDynArray);
 end;

class procedure TMySystem.func (time : double; y, p, dydt : TDoubleDynArray);
begin
  dydt[0] := p[0] - p[1]*y[0];
  dydt[1] := p[1]*y[0] - p[2]*y[1];
end;

Next, create a Runge-Kutta object:

// First argument = number of ODEs
// Second argument = number of parameters, if any
// Third argument = ODE Function
rk4 := TRK4.Create (2, 3, TMySystem.func);
rk4.p[0] := vo; rk4.p[1] := k1; rk4.p[2] := k2;

The above code include the option to declare and initialize some parameters that are part of the differential equations. To actually integrate the system by one step use the line:

// Return new time point, assign to startTime ready for next time
startTime := rk4.eval (startTime, y, stepSize);

The y argument will contain the updated solution. Run the line again to do the next step etc. Note that the eval method takes three arguments, that current value for the independent variable (time), an array that contains the current values for the dependent variable, and the step size to use.

Download the code here. RK4.zip The download include a Windows exe that you can try right away.

My First Mac OS Application

October 21, 2014 8:36 pm

I finished my first Mac OS application. Two screenshots are shown below, I'll be making it available download at the end of the week. This was a test application to see how easy it would be use write a portable application using Firemonkey. The application in question is a specialist tool that allows users to load systems biology models expressed in SBML models and simulate them. It also has some limited steady state functionality. It users libRoadRunner as the backend simulator (libroadrunner.org). Most of the underling UI business code is new but the graphing component is one I resurrected from an old Windows project I had. Took me a week or two to extract the code from the old VCL Delphi application and swap in Firemonkey canvas methods. The biggest issue I had was remembering that Firemonkey now uses single precision values to represent canvas coordinates (VCL uses integer values).

I developed the application first on Windows using XE6. Once it was working on Windows I compiled it for Mac OS. Other than a few very minor issues (eg forgetting to remove Windows in the uses clause) getting it running on the Mac was surprisingly uneventful. It took me about 20 mins to get it working on the Mac, I was impressed. Pserver, which is the conduit for moving a Delphi application to the Mac makes the process child's play. It copies over the relevant files and metadata and creates the bundle on the Mac's hard drive ready to be distributed.

Some gotchas to watch out for:

1. Scroll bars don't show automatically on the Mac in scrollable controls, to make them show up include the following line in your main form OnCreate event:

memo1.AniCalculations.AutoShowing:=false;

where memo1 is the name of a memo control.

2. Use the Deployment option in the Project menu to specify any Mac OS dylibs you want to distribute. pserver will handle the copying and bundling of these automatically.

3. Loading dll and dylibs. Not really a gotcha but I thought it would be tricky. This was surprisingly easy and in fact no different from Windows.

DllHandle :=LoadLibrary(PWideChar (fullpath));

Obviously the path to the libraries will change because the files extensions are different.

but I need this code to deal with any errors in the loading:

{$IFDEF POSIX}
errStr := string(dlerror());
{$ELSE}
errStr := SysErrorMessage(Winapi.Windows.GetLastError);
{$ENDIF}

Use of GetProcAddress is identical in Windows and Mac OS.

4. I couldn't find the equivalent of GetDeviceCaps on the Mac so assumed the same size for the Mac and Windows.

5. I didn't realize the current version of Delphi (even XE7) generates 32-bit Mac OS binaries which means it can't load any 64-bit libraries. This posed a bit of a problem initially as we only distribute libRoadRunner as a 64 bit library and the current guardian of libRoadrunner was unable to create a 32-bit version. Luckily one of my other students, (Kyle Medley) was savvy enough to compile a 32-bit version which saved the day.

6. Bugs in XE6:

a) Don't use a separate stylebook for each form because closing and freeing a second form will corrupt the buttons or other visual element on the calling form.

b) TextHeight returns incorrect an value for an empty string

7. The Firemonkey TStringGrid is severely limited or at least it seems that way. It's currently a poor cousin of this VCL partner. Not sure how to get round its limitations yet.

8. The trackar is a limited control, there is no way to alter the thumb size or shape, track width, etc.

9. Various visual artifacts on some of the styles, eg Auric has problems with thumbs in tack bars.

10. The online documentation on menu handling between Mac OS and Windows appears to be completely incorrect.

Screenshots from a Mac Pro:

What do 17 year old science students know (or not)?

Originally Posted on July 9, 2014 by hsauro

This summer I had the opportunity to host five high school students at a UW science camp. I decided to give them a project where they could build circuits that could do addition and multiplication just using resistors.

I’ve been a high school teacher in the past so I though I knew the general level of education for 17 year olds. Let me remind everyone that we’re talking about students who are one year away from entering University. Let me say that again, they are one year away from entering University. This is what I found out:

1. They’ve never seen or used a multimeter before
2. They didn’t know what the symbol was for a battery or a resistor
3. They didn’t understand what the lines were between the circuit symbols (wires!)
4. They’d never heard of a capacitor
5. They’d never seen a breadboard before and had no idea what to with it.
6. They didn’t know the difference between DC and AC (or even what they were)

What is surprising was they did know about ohms law but didn’t know the symbol for an ohm. I am certain that if I asked them to wire up a battery and a light they would have been lost.

Is this the same country that landed men on the moon?

Roman Fort at Clunderwen?

Originally Posted on August 16, 2014 by hsauro

About ten years ago, a roman road was discovered tracing a path west of Carmarthen. Currently the road is known to reach as far as Wiston where a roman fort was recently found in 2013. This was the first roman fort ever found in Pembrokeshire. Excavation at Whitland as a result of road improvements also confirmed the presence of the road (and witnessed by the writer).

While looking along the route of the roman road on Google maps I recently noticed a rectangular structure very close to the roman road as it passed south of Clunderwen. The image below was taken from Google Maps. The rectangular structure is shown inside the red square (SN12311893). The red line is the rough alignment for the roman road. North is at the top of the image. The houses to the left of the enclosure mark the southern part of Clunderwen with the railway line arcing just at the top of center of the image. The Royal Commission marks this site as a generic ‘ENCLOSURE;DEFENDED ENCLOSURE;HENGE’ but the rectangular dimensions do not match the usual shape for an iron age or post roman site. The site was only recorded in 2011 and even the Royal Commission is unsure how to interpret the site and suggests a possible roman origin. Of more interest is that the distance between Wiston where there is a confirmed roman fort and the one suggested here is between 9 and 10 miles which is a common roman distance between marching camps. The forts that extend from Edinburgh to the north of Scotland are roughly 10 miles apart. The size of the enclosure is roughly 1.3 acres. (94 meters by 54 meters). This is small for a marching camp but not unheard of. For example Knockcross in Cumbria is only 1.5 acres.  The Wiston camp is a little bigger at about 130 by 160 meters. 

The circumstantial evidence is therefore suggestive that the rectangular enclosure south of clunderwen is in fact another Pembrokeshire roman fort. Only excavation can confirm this hypothesis. 

Rhod Kemp says:

January 4, 2016 at 2:27 am

“GAER Y;FFYNNON BRODYR

Primary Reference Number (PRN) : 3729
Trust : Dyfed
Community : Clynderwen
NGR : SN12311893
Site Type (preferred type first) : Neolithic Henge / Bronze Age Henge
Legal Protection : scheduled ancient monument

Summary :
Large oval-shaped enclosure with the lowered bowl-like interior giving the impression of an amphitheatre. Defined by an exterior bank with a possible entrance on the SSE side. The field is still cultivated and the banks have been spread by ploughing but the monument remains impressive on the landscape. This site is similar to the henge at Castell Garw (1024) and is likely to be a ritual rather than domestic/defended enclosure, dating to the late Neolithic/early Bronze age. NB 2000

Description :
Character as described by Greives, & unlike any settlement enclosure I have ever seen. Lowered interior gives impression of a large amphitheatre. Similar to Fynnon Neueydd. Nantgaredig, in respect of this lowered interior, resulting drop to inside of entrance shape. Under alternating arable/pasture regime. GH Williams 1980″

from Archwilio website (http://www.cofiadurcahcymru.org.uk/arch/dyfed/english/dyfed_interface.html) 04/01/16

Skynet Edges Closer

Originally Posted on August 7, 2014 by hsauro

There is a fascinating article in science this week about the construction in hardware of a neural chip. This isn’t a new idea but scale and flexibility is novel. The chip is made by researchers from IBM and Cornell and can emulate 1,000,000 ‘neurons’. The developers claim that the approach is scalable, and efficient. Apparently chips can be cascaded to make even bigger networks.

I did a back of the envelope calculation to guestimate how many transistors it might take to emulate the brains from different organisms since we know how many transistors it took to make the chip.

This sentence was taken from Ars Technica: “The new processor, which the team is calling TrueNorth, takes a radically different approach. Its 5.4 billion transistors include over (1 million neurons) 4,000 individual cores” Assuming 100 billion neurons in a human brain that means it would take:

54 Trillion transistors to make a human brain, that is or roughly equivalent to 5,200 modern PCs, assuming 10 billion transistors per PC.

In other words not currently possible (In addition the fact that the chip only talks to 256 or neurons whereas the each brain neuron talks to about 7000). The IBM team would have to connect over 5000 chips to make a human brain, not that many it seems. Some other numbers:

Jelly Fish: ~ 4 Million transistors
Pond Snail Brain: ~ 60 Million transistors
Fruit Fly Brain: ~ 540 Million transistors
Honey Bee Brain: ~5.2 Billion transistors
Cockroach: ~5.3 Billion transistors 
Frog Brain: ~86 Billion transistors
Rat Brain: ~1 Trillion transistors
Human Brain: ~54 Trillion transistors

The single chip could just about simulate a Cockroach brain.

Brain data from http://en.wikipedia.org/wiki/List_of_an … of_neurons1 and
http://en.wikipedia.org/wiki/Neuron2

Reversibility and Product Inhibition

 Originally Posted on March 28, 2014 by hsauro

It is common in biochemistry to categorize reactions as either reversible or irreversible. What do we mean by these terms? Let’s look at a simple reaction such as:

  

where A is transformed to B. The rate of reaction for this transformation can be assumed to be first-order, that is:

  

Mathematically this is an irreversible reaction because the reaction only goes from reactant  to product . We can easily make the reaction reversible by adding the reverse reaction:

  

If we again assume first-order kinetics for the reverse process then the reaction rate is now given by:

  

Depending on the ratio of  to , if the concentration of  is high enough the reaction rate,  will actually go negative. This is the main characteristic of a reversible reaction, the potential for the forward reaction rate to go negative. If  were very small relative to  we might find that the concentration of  required to make the reaction go in reverse is unrealistically high. For example, if  and , and the concentration of  was 1 mM, then the concentration of  required to reverse the reaction would be 10 M. In a biological cell such a concentration would be highly unlikely. In this situation we would probably assume for all practical purposes that our reaction is irreversible.

Another way to look at this is to recall that the ratio of the forward to backward rate constant () is equal to the equilibrium constant, . In our example the equilibrium constant was 10,000. One may think that such equilibrium constants are quite rare but as Cornish-Bowden and Cardenas reported, the equilibrium constant for pyruvate kinase is of the order of 100,000.

An important question is how should we represent reactions that have high equilibrium constants in our computer models? Do we simply use irreversible rate laws to model such reactions? The short answer to this is probably no.

What we often forget is that the product of an enzyme catalyzed reaction can reversible bind to the enzyme and when it does it will compete with the substrate. This means that the reaction rate is likely to decrease as product accumulates simply due to the competitive effect of product. This effect is called product inhibition. For a Michaelis-Menten like rate law it is easy to adjust the equation to take into account product inhibition, we just treat the product as a competitive inhibitor:

  

where  is the maximal rate,  the substrate concentration with no product that gives half the maximal rate and  is the dissociation constant for the product. Depending on the value for  it is quite possible for the product to greatly affect the forward rate. The reaction is irreversible but the forward reaction can be reduced via product inhibition.