# Stability of a Steady State

After the calculation of the steady-state concentrations of metabolites in a reaction scheme it is often relevant to determine if the particular steady state is stable. A steady state is said to be asymptotically stable if, after a minor perturbation in system concentrations, the system returns to the original set of values. In other words, as the time since perterbation approaches infinity, the system will approach the original steady state. If the perterbation causes the system to move toward another steady state, or to not reach a steady state at all, then the steady state is said to be unstable. In the

MetabolicControlAnalysis package the function Stability determines whether a set of simultaneous differential equations, defined by S, N, v, and p, is stable or not at a given steady state. It can be shown mathematically that the presence or absence of asymptotic stability is determined by the so-called eigenvalues of the Jacobian of the system of differential equations.(1) The Jacobian is defined as follows11

If the eigenvalues all have negative real parts, the steady state is said to be asymptotically stable. If at least one of the eigenvalues has a positive real part, the steady state is unstable. In the case that some of the eigenvalues have zero real parts (with the remaining eigenvalues all having negative real parts) no conclusions can been drawn about the stability. In this final case, further analysis is needed (see Heinrich and Schuster1 and references therein).

by k

 Stability^, N, v, p, Assesses whether the system of differential equations, SteadyStateConc 0 steadystate] defined by S, N, v, and p is asymptotically stable at the steady state given by the replacement rule steadystate. Also, the function returns the eigenvalues of the Jacobian of the system of differential equations defined by S, N, v, and p. Note that the inclusion of the parameter table p is optional MMatrixftf, N, v, Calculates the Jacobian of the system of p, Normalized 0 False, differential equations defined by S, N, v, SteadyStateConc 0 steadystate] and p. Note that the last three arguments are optional; however, the default value for the option Normalized is True

Assessing the stability of a steady state.

Assessing the stability of a steady state.

Q: Assess the stability of the steady state calculated in the question/answer in Section 4.5.

A: The function Stability requires the full description of the reaction scheme, and in addition it must have the full set of steady-state concentrations. Thus, the reaction list is

Clear[Subscript, Const] eqns = { {si [t] 0 s2 [t]} , {s2[t] 0 ss[t]}, {sa[t] 0 s4[t]}};

The S, N, and v matrices and the parameters are defined as they were previously in Section 4.6.

S := SMatrix[eqns, {s1[t] , s4[t]}] ; N = NMatrix[eqns, {s1[t] , s4 [t]}] ;

p = {s1 [t] , s4 [t] , Vmax,1,f , Vmax,2,f , Vmax,3,f , Km,1,s1 , Kni,2,s2 , Kni,3,s3 } ;

If we determine the steady-state vector, we can then use the Stability function to assess the stability of the particular steady state.

Stability[S, N, v, Transpose[{p, pv}], SteadyStateConc^ sspj]

Asymptotically Stable

This output verifies that the steady state is stable, since the eigenvalues returned by the Stability function are all real and negative.

Q: Calculate the Jacobian of the reaction system

assuming that both reactions are each dependent on S1 and S2. Thus, v1 = v1 (S1, S2) and V2 = V2(S1,S2).

A: The following series of functions defines the reaction scheme and then MMatrix carries out the evaluation of the Jacobian.

S = Table[Si [t] , {i, 1, 2}] ; N = IdentityMatrix[2]; V : = Table[v[i][S] , {i, 1, 2}] ;

{{V[1]((1'0)) [{3i [t], S2 [t]}], V[1] ({0'1)) [{3i [t], 32[t]}]},

{V[2]((1'0)) [{31 [t], 32 [t]}], V [ 2 ] ({0,1)) [{31 [t], 32 [ t ]}]}}

The eigenvalues of the Jacobian are determined by using the Mathematica function EigenvalueS. The sign of the real part of these eigenvalues must be negative for the system to be stable at the given steady state.

EigenvalueS[jacob]

{ 1 (v[2]({0'1}) [{Si [t] , S2 [t]} ] + v[1] ({1'0}) [{Si [t] , S2 [t]} ] -

,(v[2]({0'1}) [{Si[t], S2 [t]}]2 - 2v[2]({0,1}) [{ S1 [ t ] , S2 [t]}]

v[1]({1,0}) [{S1[t], S2 [t]}] + v[1]({1,0}) [{S1[t], S2 [t]}]2 + 4v[1]({0,1}) [{ S1 [ t ] , S2 [t]}] v [ 2 ]({1,0}) [{ S1 [ t ] , S2 [t]}])) ,

1 (v [ 2 ] ({0,1}) [{S1 [t] , S2 [t]} ] + v[1]({1'0})[{S1[t], S2 [t]} ] +

, (v [ 2 ]({0,1}) [{ S1 [ t ] , S2 [t]}]2 - 2v[2]({0,1}) [{ S1 [ t ] , S2 [t]}]

v[1]({1,0}) [{S1[t], S2 [t]}] + v[1]({1,0}) [{S1[t], S2 [t]}]2 + 4v[1]({0,1}) [{ S1 [ t ] , S2 [t]}] v [ 2 ]({1,0}) [{ S1 [ t ] , S2 [ t ]}])) =

Clearly, the sign of the eigenvalues is dependent on the expressions for the velocities and the values of their parameters, as well as the actual steady-state values themselves.

Q: What is the steady-state concentration of Si in the reaction scheme defined below (from Heinrich and Shuster,(1) Scheme 4, p. 47)? The scheme involves the supply of the metabolite, Si , and two pathways of removal, with one displaying high-substrate inhibition of the enzyme that catalyzes it.

Suppose that v1= k1, v2= k2[S1], and v3= ^n, with k1 = 0.9 mmol L-1 h 1, k2 =

A: First we need to define the appropriate matrices, reaction rates, and parameter values for this system.

Clear[Subscript];

kssi[t]

p = {ki ,k2 ,k3 ,Ki , n} ; pv = {0.9, 0.11, 0.4, 3, 4},

Then the steady state is calculated as follows:

{{si[t] 0-2.33223 - 2.7614 i}, {si[t] 0-2.33223 + 2.7614 i}, {si[t] 0 2.05637}, {s1[t] 0 3.28846}, {s1[t] 0 7.50144}}

By selecting only the real solutions we find three physically relevant steady states. realss = ss /. {a_, b_ -> c_Complex, d_} -> {} // Flatten

{s1[t] 0 2.05637, s1[t] 0 3.28846, s1[t] 0 7.50144}

We can now use a Do statement around the Stability function to assess the stability of each steady state individually.

Do[Stability[S, N, v, Transpose[{p, pv}] , SteadyStateConc-> realsspj] , {i, 3}]

Asymptotically Stable

Asymptotically Unstable

{0.113127}

Asymptotically Stable

The stability analysis reveals that the second steady state in the list is asymptotically unstable. This is illustrated in the graph from the following simulation where the initial concentration of S1 is set to the second value. (More significant figures are used than are shown above.)

{s1 [0] == 3.28846091596}, {t, 0, 800}, Transpose[{p, pv}]];

AxesLabel -> {"Time (h)", "Concentration (mM)"}, PlotRange ^ {1, 4}] ;

Concentration (mM) 4

200 400 600 800

Figure 4.5. Time course of the reaction scheme shown in Eqn [4.13].

The graph shows that the concentration is stable at ~3.3 mmol L 1 and after ~320 h it undergoes a transition to a new steady state at ~2.0 mmol L-1.