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

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

**Output:**

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

**Output:**

```
vol = 1e-19, a.concInit = 1.0, a.nInit = 60221.415
```

```
vol = 1e-20, a.concInit = 1.0, a.nInit = 6022.1415
```

```
vol = 1e-21, a.concInit = 1.0, a.nInit = 602.21415
```

```
vol = 3e-22, a.concInit = 1.0, a.nInit = 180.664245
```

```
vol = 1e-22, a.concInit = 1.0, a.nInit = 60.221415
```

```
vol = 3e-23, a.concInit = 1.0, a.nInit = 18.0664245
```

```
vol = 1e-23, a.concInit = 1.0, a.nInit = 6.0221415
```

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

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

**Output:**

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

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

**Output:**

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

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

**Output:**

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

**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
```

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

**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
```