Computing the steady-state using Tellurium

 If you're building a model and you want to quickly find the model's steady state, you can call the command steadyState. Let's illustrate this with an example:

import tellurium as te

r = te.loada ('''
     # Define a simple linear chain of reactions
     $Xo -> S1; k1*Xo;
      S1 -> S2; k2*S1;
      S2 -> S3; k3*S2;
      S3 -> $X1; k4*S3;

  # Initialize rate constants 
  k1 = 0.2; k2 = 0.7; k3 = 0.15; k4 = 1.3;
  Xo = 10
''')

print r.getSteadyStateValues() 

Running this script by clicking in the green arrow in the tool bar will output the following steady state levels:

[2.85714286  13.33333333   1.53846154]

But which steady state value refers to which species in my model? To find that out call getFloatingSpeciesIds:

print r.getFloatingSpeciesIds()
['S1', 'S2', 'S3']

The order of the speces Ids match the order of steady state values. In other words, S1 = 2.857, S2 = 13.333, and S3 = 1.538.

If you want to check that these values really are the steady state values you can compute the rates of change:

print r.getRatesOfChange()
array([  0.00000000e+00,  -4.44089210e-16,   4.44089210e-16])

If you look carefully all the rates of change are effectively zero, thus confirming we're at steady state.

What about the stability of the steady state? That is, is it stable to small perturbations in S1, S2 or S3? To find this out we need to compute the eigenvalues of the system Jacobian matrix. If all the eigenvalues are negative this means small pertubrations are damped so that the system returns to the steady state. There are two ways to do this, we can use the getReducedEigenValues call, or we get compute the Jacobian first and then compute the eigenvalues of the Jacobian. Both these approaches are given below. Its probably simplest to call just getReducedEigenValues.

print r.getReducedEigenValues()
[-0.7  -0.15 -1.3 ]

Notice that all the eigenvalues are negative indicating that small perturbations are stable.

This is the alternative approach which computes the Jacobian first and then the eigenvalues of the Jacobian:

print te.getEigenvalues(r.getFullJacobian())
[-0.7  -0.15 -1.3 ]

Note that we're using a utility method from the tellurium library, getEigenValues to compute the eigenvalues.