Chemical Bistables

A bistable system is a dynamic system that has two stable equilibrium states. The following examples can be used to teach and demonstrate different aspects of bistable systems or to learn how to model them using moose. Each example contains a short description, the model's code, and the output with default settings.

Each example can be found as a python file within the main moose folder under

(...)/moose/moose-examples/tutorials/ChemicalBistables

In order to run the example, run the script

python filename.py

in command line, where filename.py is the name of the python file you would like to run. The filenames of each example are written in bold at the beginning of their respective sections, and the files themselves can be found in the aformentioned directory.

In chemical bistable models that use solvers, there are optional arguments that allow you to specify which solver you would like to use.

python filename.py [gsl | gssa | ee]

Where:

  • gsl: This is the Runge-Kutta-Fehlberg implementation from the GNU Scientific Library (GSL). It is a fifth order variable timestep explicit method. Works well for most reaction systems except if they have very stiff reactions.
  • gssl: Optimized Gillespie stochastic systems algorithm, custom implementation. This uses variable timesteps internally. Note that it slows down with increasing numbers of molecules in each pool. It also slows down, but not so badly, if the number of reactions goes up.
  • Exponential Euler:This methods computes the solution of partial and ordinary differential equations.

All the following examples can be run with either of the three solvers, each of which has different advantages and disadvantages and each of which might produce a slightly different outcome.

Simply running the file without the optional argument will by default use the gsl solver. These gsl outputs are the ones shown below.



Simple Bistables

Filename: simpleBis.py

This example shows the key property of a chemical bistable system: it has two stable states. Here we start out with the system settling rather quickly to the first stable state, where molecule A is high (blue) and the complementary molecule B (green) is low. At t = 100s, we deliver a perturbation, which is to move 90% of the A molecules into B. This triggers a state flip, which settles into a distinct stable state where there is more of B than of A. At t = 200s we reverse the flip by moving 99% of B molecules back to A.

If we run the simulation with the gssa option python simpleBis.py gssa

we see exactly the same sequence of events, except now the switch is noisy. The calculations are now run with the Gillespie Stochastic Systems Algorithm (gssa) which incorporates probabilistic reaction events. The switch still switches but one can see that it might flip spontaneously once in a while.

Things to do:

1. Open a copy of the script file in an editor, and around line 124 and 129 you will see how the state flip is implemented while maintaining mass conservation. What happens if you flip over fewer molecules? What is the threshold for a successful flip? Why are these thresholds different for the different states?

  1. Try different volumes in line 31, and rerun using the gssa. Will you see more or less noise if you increase the volume to 1e-20 m^3?

Code:

Show/Hide code

Output:

../../../_images/simpleB.png

Scale Volumes

File name: scaleVolumes.py

This script runs exactly the same model as in simpleBis.py, but it automatically scales the volumes from 1e-19 down to smaller values.

Note how the simulation successively becomes noisier, until at very small volumes there are spontaneous state transitions.

Code:

Show/Hide code

Output:

vol = 1e-19, a.concInit = 1.0, a.nInit = 60221.415
../../../_images/sV1.png
vol = 1e-20, a.concInit = 1.0, a.nInit = 6022.1415
../../../_images/sV2.png
vol = 1e-21, a.concInit = 1.0, a.nInit = 602.21415
../../../_images/sV3.png
vol = 3e-22, a.concInit = 1.0, a.nInit = 180.664245
../../../_images/sV4.png
vol = 1e-22, a.concInit = 1.0, a.nInit = 60.221415
../../../_images/sV5.png
vol = 3e-23, a.concInit = 1.0, a.nInit = 18.0664245
../../../_images/sV6.png
vol = 1e-23, a.concInit = 1.0, a.nInit = 6.0221415
../../../_images/sV7.png

Strong Bistable System

File name: strongBis.py

This example illustrates a particularly strong, that is, parametrically robust bistable system. The model topology is symmetric between molecules b and c. We have both positive feedback of molecules b and c onto themselves, and also inhibition of b by c and vice versa.

../../../_images/strongBis.png

Open the python file to see what is happening. The simulation starts at a symmetric point and the model settles at precisely the balance point where a, b, and c are at the same concentration. At t= 100 we apply a small molecular 'tap' to push it over to a state where c is larger. This is stable. At t = 210 we apply a moderate push to show that it is now very stably in this state, and the system rebounds to its original levels. At t = 320 we apply a strong push to take it over to a state where b is larger. At t = 430 we give it a strong push to take it back to the c dominant state.

Code:

Show/Hide code

Output:

../../../_images/strongB.png

MAPK Feedback Model

File name: mapkFB.py

This example illustrates loading, and running a kinetic model for a much more complex bistable positive feedback system, defined in kkit format. This is based on Bhalla, Ram and Iyengar, Science 2002.

The core of this model is a positive feedback loop comprising of the MAPK cascade, PLA2, and PKC. It receives PDGF and Ca2+ as inputs.

../../../_images/mapkFB.png

This model is quite a large one and due to some stiffness in its equations, it takes about 30 seconds to execute. Note that this is still 200 times faster than the events it models.

The simulation illustrated here shows how the model starts out in a state of low activity. It is induced to 'turn on' when a a PDGF stimulus is given for 400 seconds, starting at t = 500s. After it has settled to the new 'on' state, the model is made to 'turn off' by setting the system calcium levels to zero. This inhibition starts at t = 2900 and goes on for 500 s.

Note that this is a somewhat unphysiological manipulation! Following this the model settles back to the same 'off' state it was in originally.

Code:

Show/Hide code

Output:

../../../_images/mapkFB2.png

Propogation of a Bistable System

File name: propagationBis.py

All the above models have been well-mixed, that is point or non-spatial models. Bistables do interesting things when they are dispersed in space. This is illustrated in this example. Here we have a tapering cylinder, that is a pseudo 1-dimensional reaction-diffusion system. Every point in this cylinder has the bistable system from the strongBis example.

../../../_images/propBis.png

The example has two stages. First it starts out with the model in the unstable transition point, and introduces a small symmetry-breaking perturbation at one end. This rapidly propagates through the entire length model, leaving molecule b at a higher value than c.

At t = 100 we carry out a different manipulation. We flip the concentrations of molecules b and c for the left half of the model, and then just let it run. Now we have opposing bistable states on either half. In the middle, the two systems battle it out. Molecule c from the left side diffuses over to the right, and tries to inhibit b, and vice versa. However we have a small asymmetry due to the tapering of the cylinder. As there is a slightly larger volume on the left, the transition point gradually advances to the right, as molecule b yields to the slightly larger amounts of molecule c.

Code:

Show/Hide code

Output:

../../../_images/propBis.gif

Steady-state Finder

File name: findSteadyState

This is an example of how to use an internal MOOSE solver to find steady states of a system very rapidly. The method starts from a random position in state space that obeys mass conservation. It then finds the nearest steady state and reports it. If it does this enough times it should find all the steady states.

We illustrate this process for 50 attempts to find the steady states. It does find all of them. Each time it plots and prints the values, though the plotting is not necessary.

The printout shows the concentrations of all molecules in the first 5 columns. Then it prints the type of solution, and the numbers of negative and positive eigenvalues. In all cases the calculations are successful, though it takes different numbers of iterations to arrive at the steady state. In some models it would be necessary to put a cap on the number of iterations, if the system is not able to find a steady state.

In this example we run the bistable model using the ODE solver right at the end, and manually enforce transitions to show where the target steady states are.

For more information on the algorithm used, look in the comments within the main method of the code below.

Code:

Show/Hide code

Output:

0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=16
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=29
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=10
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=26
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=27
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=30
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=12
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=29
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=12
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=41
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=29
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=18
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=27
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=14
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=12
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=19
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter= 6
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=14
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=23
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=25
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=16
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter= 5
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=43
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter= 9
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=43
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=29
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=27
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter= 9
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=12
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=24
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=26
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=14
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=14
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=10
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=13
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=26
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=21
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=26
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=24
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=24
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=18
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=26
0.18 0.75 0.00 0.03 0.01 Type=5 NegEig=4 PosEig=0 status=0 success Iter=13
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=23
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=24
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter= 8
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=0 status=0 success Iter=18
0.18 0.75 0.00 0.03 0.01 Type=0 NegEig=3 PosEig=1 status=0 success Iter=21
0.99 0.00 0.01 0.00 0.00 Type=0 NegEig=3 PosEig=0 status=0 success Iter=15
0.92 0.05 0.00 0.01 0.01 Type=2 NegEig=2 PosEig=1 status=0 success Iter=29
../../../_images/findS.png

Dose Response (Under construction)

File name: doseResponse.py

This example generates a doseResponse plot for a bistable system, against a control parameter (dose) that takes the system in and out again from the bistable regime. Like the previous example, it uses the steady-state solver to find the stable points for each value of the control parameter. Unfortunately it doesn't work right now. Seems like the kcat scaling isn't being registered.

Code:

Show/Hide code

Output:

scale=0.010 kcat=0.004
scale=0.013 kcat=0.005
scale=0.016 kcat=0.006
scale=0.020 kcat=0.007
scale=0.025 kcat=0.009
scale=0.032 kcat=0.011
scale=0.040 kcat=0.014
scale=0.050 kcat=0.018
scale=0.063 kcat=0.023
scale=0.079 kcat=0.029
scale=0.100 kcat=0.036
scale=0.126 kcat=0.045
scale=0.158 kcat=0.057
scale=0.200 kcat=0.072
scale=0.251 kcat=0.091
scale=0.316 kcat=0.114
scale=0.398 kcat=0.144
scale=0.501 kcat=0.181
scale=0.631 kcat=0.228
scale=0.794 kcat=0.287
scale=1.000 kcat=0.361
scale=1.259 kcat=0.454
scale=1.585 kcat=0.572
scale=1.995 kcat=0.720
scale=2.512 kcat=0.907
scale=3.162 kcat=1.142
scale=3.981 kcat=1.437
scale=5.012 kcat=1.809
scale=6.310 kcat=2.278
scale=7.943 kcat=2.868
scale=10.000        kcat=3.610
scale=12.589        kcat=4.545
scale=15.849        kcat=5.722
scale=19.953        kcat=7.203
scale=25.119        kcat=9.068
scale=31.623        kcat=11.416
scale=39.811        kcat=14.372
scale=50.119        kcat=18.093
scale=63.096        kcat=22.778
scale=79.433        kcat=28.676
scale=100.000       kcat=36.101
../../../_images/doseR.png