Finding Steady States through Homotopy Continuation

The steady states of a dynamical system ${dx \over dt} = f(x)$ can be found by solving $0 = f(x)$. This is typically a hard problem, and generally, there is no method guaranteed to find all steady states for a system that has multiple ones. However, many chemical reaction networks generate polynomial systems (for example those which are purely mass action or have only have Hill functions with integer Hill exponents). The roots of these can reliably be found through a homotopy continuation algorithm^[1][2]. This is implemented in Julia through the HomotopyContinuation.jl package^[3].

Catalyst contains a special homotopy continuation extension that is loaded whenever HomotopyContinuation.jl is. This exports a single function, hc_steady_states, that can be used to find the steady states of any ReactionSystem structure.

For this tutorial, we will use the Wilhelm model (which demonstrates bistability in a small chemical reaction network). We declare the model and the parameter set for which we want to find the steady states:

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
ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5]

Here, we only run import HomotopyContinuation as we do not require any of its functions, and just need it to be present in the current scope for the extension to be activated.

Now we can find the steady states using:

hc_steady_states(wilhelm_2009_model, ps)

3-element Vector{Vector{Float64}}:
 [4.5, 6.0]
 [0.0, 0.0]
 [0.4999999999999999, 1.9999999999999998]

The order of the species in the output vectors are the same as in species(wilhelm_2009_model).

It should be noted that the steady state multivariate polynomials produced by reaction systems may have both imaginary and negative roots, which are filtered away by hc_steady_states. If you want the negative roots, you can use the hc_steady_states(wilhelm_2009_model, ps; filter_negative=false) argument.

Systems with conservation laws

Some systems are under-determined, and have an infinite number of possible steady states. These are typically systems containing a conservation law, e.g.

two_state_model = @reaction_network begin
    (k1,k2), X1 <--> X2
end

\[ \begin{align*} \mathrm{\mathtt{X1}} &\xrightleftharpoons[\mathtt{k2}]{\mathtt{k1}} \mathrm{\mathtt{X2}} \end{align*} \]

Catalyst allows the conservation laws of such systems to be computed using the conservationlaws function. By using these to reduce the dimensionality of the system, as well as specifying the initial amount of each species, homotopy continuation can again be used to find steady states. Here we do this by designating such an initial condition (which is used to compute the system's conserved quantities, in this case $X1 + X2$):

ps = [:k1 => 2.0, :k2 => 1.0]
u0 = [:X1 => 1.0, :X2 => 1.0]
hc_steady_states(two_state_model, ps; u0)

1-element Vector{Vector{Float64}}:
 [0.6666666666666666, 1.3333333333333335]

Final notes

• hc_steady_states support any systems where all rates are systems of rational polynomials (such as Hill functions with integer Hill coefficients).
• When providing initial conditions to compute conservation laws, values are only required for those species that are part of conserved quantities. If this set of species is unknown, it is recommended to provide initial conditions for all species.
• Additional arguments provided to hc_steady_states are automatically passed to HomotopyContinuation's solve command. Use e.g. show_progress = false to disable the progress bar.

Citation

If you use this functionality in your research, please cite the following paper to support the authors of the HomotopyContinuation package:

@inproceedings{HomotopyContinuation.jl,
  title={{H}omotopy{C}ontinuation.jl: {A} {P}ackage for {H}omotopy {C}ontinuation in {J}ulia},
  author={Breiding, Paul and Timme, Sascha},
  booktitle={International Congress on Mathematical Software},
  pages={458--465},
  year={2018},
  organization={Springer}
}

References