Smallest Chemical Reactions Systems that is Bistable

Originally Posted on May 15, 2018 by hsauro

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.

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

import tellurium as te import roadrunner   r = te.loada('''     // Reactions:     J0:     Y -> X + X;    k1*S*Y     J1: X + X -> X + Y;    k2*X^2     J3: X + Y -> Y + $P;   k3*X*Y     J4: X     -> $P   ;    k4*X       // Species Init:     S = 0.8; X = 3; Y = 3; P = 1;       // Variable Init:     k1 = 8; k2 = 1; k3 = 1; k4 = 1.5; ''') m = r.simulate (0, 10, 100)

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:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

import tellurium as te import roadrunner import matplotlib.pyplot as plt from rrplugins import *   r = te.loada ('''     // Reactions:     J0:     Y -> X + X;    k1*S*Y     J1: X + X -> X + Y;    k2*X^2     J3: X + Y -> Y + $P;   k3*X*Y     J4: X     -> $P   ;    k4*X     J5: -> X; 0.6;  // Comment this out to get the original plot       // Species Init:     S = 0.8; X = 3; Y = 3; P = 1;       // Variable Init:     k1 = 8; k2 = 1; k3 = 1; k4 = 1.5;   ''') m = r.simulate (0, 10, 100)   r.steadyState()   # Load the bifurcxation plugin auto = Plugin("tel_auto2000")   # Setup properties auto.setProperty("SBML", r.getCurrentSBML()) auto.setProperty("ScanDirection", "Negative") auto.setProperty("PrincipalContinuationParameter", "S") auto.setProperty("PreSimulation", "True") auto.setProperty("PreSimulationDuration", 1.0) auto.setProperty("RL0", 0.2) auto.setProperty("RL1", 2)   # execute the plugin auto.execute()   # plot Bifurcation diagram pts     = auto.BifurcationPoints lbls    = auto.BifurcationLabels biData  = auto.BifurcationData biData.plotBifurcationDiagram(pts, lbls)   # If you want the raw data use: rawData = auto.BifurcationData.toNumpy   # If you want a summary of the bifurcation analysis use: print auto.BifurcationSummary