Steady state solutions

Steady state solutions#

Solving steady state solutions for an ODE system is to find a combination of state variables such that their derivatives are all zeroes. Their are two ways:

  • Running the ODE solver until a steady state is reached (DynamicSS() and DifferentialEquations.jl)

  • Using a root-finding algorithm to find a steady state (SSRootfind() and NonlinearSolve.jl)

Defining a steady state problem#

using SteadyStateDiffEq
using OrdinaryDiffEq

model(u, p, t) = 1 - p * u

p = 1.0
u0 = 0.0
prob = SteadyStateProblem(model, u0, p)
alg = DynamicSS(Tsit5())

SteadyStateDiffEq.DynamicSS{OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Float64}(OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}(OrdinaryDiffEqCore.trivial_limiter!, OrdinaryDiffEqCore.trivial_limiter!, static(false)), Inf)

Solve the problem. The result should be close to 1.0.

sol = solve(prob, alg)

retcode: Success
u: 0.9999994687790571

Modelingtoolkit solving nonlinear systems#

Use NonlinearSolve.jl and NonlinearSystem().

using ModelingToolkit
using NonlinearSolve

@variables x y z
@parameters σ ρ β

eqs = [
    0 ~ σ * (y - x),
    0 ~ x * (ρ - z) - y,
    0 ~ x * y - β * z
]

@mtkbuild ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])

guess = [x => 1.0, y => 0.0, z => 0.0]
ps = [σ => 10.0, ρ => 26.0, β => 8 / 3]
prob = NonlinearProblem(ns, guess, ps)
sol = solve(prob, NewtonRaphson()) ## The results should be all zeroes

retcode: Success
u: 2-element Vector{Float64}:
 -2.129924444096732e-29
 -2.398137151871876e-28

Another nonlinear system example with structural_simplify().

@independent_variables t
@variables u1(t) u2(t) u3(t) u4(t) u5(t)

eqs = [
    0 ~ u1 - sin(u5)
    0 ~ u2 - cos(u1)
    0 ~ u3 - hypot(u1, u2)
    0 ~ u4 - hypot(u2, u3)
    0 ~ u5 - hypot(u4, u1)
]

@mtkbuild sys = NonlinearSystem(eqs, [u1, u2, u3, u4, u5], [])
\[ \begin{align} 0 &= \mathtt{u5}\left( t \right) - hypot\left( \mathtt{u4}\left( t \right), \mathtt{u1}\left( t \right) \right) \end{align} \]

You can simplify the problem using structural_simplify(). There will be only one state variable left. The solve can solve the problem faster.

prob = NonlinearProblem(sys, [u5 => 0.0])
sol = solve(prob, NewtonRaphson())

retcode: Success
u: 1-element Vector{Float64}:
 1.6069926947050055

The answer should be 1.6 and 1.0.

@show sol[u5] sol[u1];

sol[u5] = 1.6069926947050055
sol[u1] = 0.9993449829954304

Finding Steady States through Homotopy Continuation#

For systems of rational polynomials. e.g., mass-action reactions and reaction rates with integer Hill coefficients. Source: https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/homotopy_continuation/#homotopy_continuation

using Catalyst
import HomotopyContinuation

wilhelm_2009_model = @reaction_network begin
    k1, Y --> 2X
    k2, 2X --> X + Y
    k3, X + Y --> Y
    k4, X --> 0
end
\[\begin{split} \begin{align*} \mathrm{Y} &\xrightarrow{\mathtt{k1}} 2 \mathrm{X} \\ 2 \mathrm{X} &\xrightarrow{\mathtt{k2}} \mathrm{X} + \mathrm{Y} \\ \mathrm{X} + \mathrm{Y} &\xrightarrow{\mathtt{k3}} \mathrm{Y} \\ \mathrm{X} &\xrightarrow{\mathtt{k4}} \varnothing \end{align*} \end{split}\]

There are three steady states

ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5]
steady_states = hc_steady_states(wilhelm_2009_model, ps)

Computing mixed cells... 2    Time: 0:00:00
  mixed_volume:  3



Computing mixed cells... 2    Time: 0:00:00
  mixed_volume:  3

Tracking 3 paths...  67%|████████████████████▋          |  ETA: 0:00:09
  # paths tracked:                  2
  # non-singular solutions (real):  2 (2)
  # singular endpoints (real):      0 (0)
  # total solutions (real):         2 (2)












Tracking 3 paths... 100%|███████████████████████████████| Time: 0:00:17
  # paths tracked:                  3
  # non-singular solutions (real):  3 (3)
  # singular endpoints (real):      0 (0)
  # total solutions (real):         3 (3)

3-element Vector{Vector{Float64}}:
 [0.5, 2.0]
 [0.0, 0.0]
 [4.500000000000001, 6.000000000000001]

Steady state stability computation#

Catalyst.steady_state_stability() shows there are two stable and one unstable steady states.

using Catalyst
[Catalyst.steady_state_stability(sstate, wilhelm_2009_model, ps) for sstate in steady_states]

3-element Vector{Bool}:
 0
 1
 1

This notebook was generated using Literate.jl.