Steady state stability computation

After system steady states have been found using HomotopyContinuation.jl, NonlinearSolve.jl, or other means, their stability can be computed using Catalyst's steady_state_stability function. Systems with conservation laws will automatically have these removed, permitting stability computation on systems with singular Jacobian.

Basic examples

Let us consider the following basic example:

using Catalyst
rn = @reaction_network begin
    (p,d), 0 <--> X
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{p} \mathrm{X} \end{align*} \]

It has a single (stable) steady state at $X = p/d$. We can confirm stability using the steady_state_stability function, to which we provide the steady state, the reaction system, and the parameter values:

ps = [:p => 2.0, :d => 0.5]
steady_state = [:X => 4.0]
steady_state_stability(steady_state, rn, ps)

true

Next, let us consider the following self-activation loop:

sa_loop = @reaction_network begin
    (hill(X,v,K,n),d), 0 <--> X
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{\frac{X^{n} v}{K^{n} + X^{n}}} \mathrm{X} \end{align*} \]

For certain parameter choices, this system exhibits multi-stability. Here, we can find the steady states using homotopy continuation:

import HomotopyContinuation
ps = [:v => 2.0, :K => 0.5, :n => 3, :d => 1.0]
steady_states = hc_steady_states(sa_loop, ps)

3-element Vector{Vector{Float64}}:
 [0.26870078851261286]
 [0.0]
 [1.967716165985015]

Next, we can apply steady_state_stability to each steady state yielding a vector of their stabilities:

[steady_state_stability(sstate, sa_loop, ps) for sstate in steady_states]

3-element Vector{Bool}:
 0
 1
 1

Finally, as described above, Catalyst uses an optional argument, tol, to determine how strict to make the stability check. I.e. below we set the tolerance to 1e-6 (a larger value, that is stricter, than the default of 10*sqrt(eps()))

[steady_state_stability(sstate, sa_loop, ps; tol = 1e-6) for sstate in steady_states]

Pre-computing the Jacobian to increase performance when computing stability for many steady states

Catalyst uses the system Jacobian to compute steady state stability, and the Jacobian is computed once for each call to steady_state_stability. If you repeatedly compute stability for steady states of the same system, pre-computing the Jacobian and supplying it to the steady_state_stability function can improve performance.

In this example we use the self-activation loop from previously, pre-computes its Jacobian, and uses it to multiple steady_state_stability calls:

ss_jac = steady_state_jac(sa_loop)

ps_1 = [:v => 2.0, :K => 0.5, :n => 3, :d => 1.0]
steady_states_1 = hc_steady_states(sa_loop, ps)
stabs_1 = [steady_state_stability(st, sa_loop, ps_1; ss_jac) for st in steady_states_1]

ps_2 = [:v => 4.0, :K => 1.5, :n => 2, :d => 1.0]
steady_states_2 = hc_steady_states(sa_loop, ps)
stabs_2 = [steady_state_stability(st, sa_loop, ps_2; ss_jac) for st in steady_states_2]