Another way to find unstable steady states

James Glazier recently told me of a trick he uses to find unstable steady states. Consider the following model which has two stable and one unstable steady state (it's a bistable system using postiive feedback). Note that $Xo means the species Xo is fixed.

import tellurium as te

r = te.loada('''
    $Xo -> S1; (0.1 + k1*S1^4/(k2+S1^4));
    S1 ->; k3*S1;
    
    k1 = 0.9; k2 = 0.3; k3 = 0.7; S1 = 0.5;
''')

If we run a simulation of this system it evolves to one of the stable steady states, in this case 0.144635. If we set the initial conditon to S1 = 10, we can also get the other stable state at S1 = 1.3095. There is the code to do that:

import tellurium as te

r = te.loada('''
    $Xo -> S1; ( 0.1 + k1*S1^4/(k2+S1^4));
    S1 ->; k3*S1;
    
    k1 = 0.9; k2 = 0.3; k3 = 0.7; S1 = 0.5;
''')

r.steadyState()
print (r.S1)

# Find theother steady state
r.S1 = 10
r.steadyState()
print (r.S1)

But how can we find the unstable one? There as an old trick where if one integrates backwards in time, stable states became repelers and unstable states attractors. However we don't allow someone to specfiy a start time that is bigger then the eed time. Instead James Glazier realized one chould just a put minus sign in front of every rate law to mimic the same effect. For example, a simulation of the following modiified model:

import tellurium as te

r = te.loada('''
    $Xo -> S1; -(0.1 + k1*S1^4/(k2+S1^4));
    S1 ->; -k3*S1;
    
    k1 = 0.9; k2 = 0.3; k3 = 0.7; S1 = 0.5;
''')
m = r.simulate (0, 50, 100)
print (r.S1)

will yield the unstable state at S1 = 0.68256.