Here is a simple bistable system I sometimes use in class. I like it because it’s easy to explain plus the addition of one extra reaction turns it into a relaxation oscillator.

import tellurium as te import random import pylab as plt import numpy as np # ---- BISTABLE SYSTEM ----- r = te.loada (''' $Xo -> S1; 0.5 + Vmax*S1^n/(15 + S1^n); S1 -> ; k1*S1; Xo = 1; S1 = 1; n = 4; Vmax = 10; k1 = 2; ''') plt.xlabel ('Time') plt.ylabel ('Concentration') for S1 in np.arange (1.45, 1.55, 0.01): r.S1 = S1 m = r.simulate (0, 5, 100) plt.plot (m['time'], m['[S1]']) r.steadyStateSolver.allow_presimulation = False for i in range (10): r.S1 = random.random()*5 r.steadyState() print ('S1 = ', r.S1, 'eigenvalue = ', r.getFullEigenValues()[0])

If you run the above code you’ll get:

S1 = 5.145227 eigenvalue = -1.84055

S1 = 5.145227 eigenvalue = -1.84055

S1 = 0.251329 eigenvalue = -1.95768

S1 = 5.145227 eigenvalue = -1.84055

S1 = 1.492258 eigenvalue = 3.004970

S1 = 5.145227 eigenvalue = -1.84055

S1 = 1.492258 eigenvalue = 3.004970

S1 = 5.145227 eigenvalue = -1.84055

S1 = 5.145227 eigenvalue = -1.84055

S1 = 0.251329 eigenvalue = -1.957684

The above output is what I call the poor-mans way of finding multiple steady states. I scan across the inital conditions computing the steady state in each case and the corresponding eigenvalue (You can also just assign random values to the initial value for S1, which is better if you don’t know what starting conditions to use). As you can see it found three steady states, two stable (neg eignvalue) and one unstable (1.492, pos eigenvalue).

The plot above shows time course trajectories for about 10 different starting points for S1. The two steady states are clearly evident.

Out of interest here is the addition of the extra reaction at the front to turn it into a relaxation oscillator. You have to make a slight modification to the second reaction’s rate law so that it becomes a function of its reactant (which it wasn’t before)

$Xo -> So; k0*Xo; So -> S1; k1*So + Vmax*So*S1^n/(15 + S1^n); S1 -> ; k2*S1; Xo = 1; S1 = 1; k0 = 0.07; k1 = 0.04 n = 4; Vmax = 6; k2 = 0.1;