Finding Steady States using NonlinearSolve.jl and SteadyStateDiffEq.jl

Catalyst ReactionSystem models can be converted to ODEs (through the reaction rate equation). We have previously described how these ODEs' steady states can be found through homotopy continuation. Generally, homotopy continuation (due to its ability to find all of a system's steady states) is the preferred approach. However, Catalyst supports two additional approaches for finding steady states:

Through solving the nonlinear system produced by setting all ODE differentials to 0^[1].
Through forward ODE simulation from an initial condition until a steady state has been reached.

While these approaches only find a single steady state, they offer two advantages as compared to homotopy continuation:

They are typically much faster.
They can find steady states for models that do not produce multivariate, rational, polynomial systems (which is a requirement for homotopy continuation to work). Examples include models with non-integer hill coefficients.

In practice, model steady states are found through nonlinear system solving by creating a NonlinearProblem, and through forward ODE simulation by creating a SteadyStateProblem. These are then solved through solvers implemented in the NonlinearSolve.jl, package (with the latter approach also requiring the SteadyStateDiffEq.jl package). This tutorial describes how to find steady states through these two approaches. More extensive descriptions of available solvers and options can be found in NonlinearSolve's documentation.

Steady state finding through nonlinear solving

Let us consider a simple dimerisation system, where a protein ($P$) can exist in a monomer and a dimer form. The protein is produced at a constant rate from its mRNA, which is also produced at a constant rate.

using Catalyst
dimer_production = @reaction_network begin
    pₘ, 0 --> mRNA
    pₚ, mRNA --> mRNA + P
    (k₁, k₂), 2P <--> P₂
    d, (mRNA, P, P₂) --> 0
end

\[ \begin{align*} \varnothing &\xrightarrow{\mathtt{p_m}} \mathrm{\mathtt{mRNA}} \\ \mathrm{\mathtt{mRNA}} &\xrightarrow{\mathtt{p_p}} \mathrm{\mathtt{mRNA}} + \mathrm{P} \\ 2 \mathrm{P} &\xrightleftharpoons[\mathtt{k_2}]{\mathtt{k_1}} \mathrm{\mathtt{P_2}} \\ \mathrm{\mathtt{mRNA}} &\xrightarrow{d} \varnothing \\ \mathrm{P} &\xrightarrow{d} \varnothing \\ \mathrm{\mathtt{P_2}} &\xrightarrow{d} \varnothing \end{align*} \]

This system corresponds to the following ODE:

\[\begin{aligned} \frac{dmRNA}{dt} &= pₘ - d \cdot mRNA \\ \frac{dP}{dt} &= pₚ \cdot mRNA - k₁ \cdot P + 2k₂ \cdot P₂ - d \cdot P \\ \frac{dP₂}{dt} &= k₁ \cdot P + 2k₂ \cdot P₂ \\ \end{aligned}\]

To find its steady states we need to solve:

\[\begin{aligned} 0 &= pₘ - d \cdot mRNA \\ 0 &= pₚ \cdot mRNA - k₁ \cdot P + 2k₂ \cdot P₂ - d \cdot P \\ 0 &= k₁ \cdot P + 2k₂ \cdot P₂ \\ \end{aligned}\]

To solve this problem, we must first designate our parameter values, and also make an initial guess of the solution. Generally, for problems with a single solution (like this one), most arbitrary guesses will work fine (the exception typically being systems with conservation laws). Using these, we can create the NonlinearProblem that we wish to solve.

p = [:pₘ => 0.5, :pₚ => 2.0, :k₁ => 5.0, :k₂ => 1.0, :d => 1.0]
u_guess = [:mRNA => 1.0, :P => 1.0, :P₂ => 1.0]
nlprob = NonlinearProblem(dimer_production, u_guess, p)

Finally, we can solve it using the solve command, returning the steady state solution:

using NonlinearSolve
sol = solve(nlprob)

retcode: Success
u: 3-element Vector{Float64}:
 0.5
 0.46332495807108
 0.26833752096446006

Typically, a good default method is automatically selected for any problem. However, NonlinearSolve does provide a wide range of potential solvers. If we wish to designate one, it can be supplied as a second argument to solve. Here, we use the Newton Trust Region method, and then check that the solution is equal to the previous one.

sol_ntr = solve(nlprob, TrustRegion())
sol ≈ sol_ntr

true

Systems with conservation laws

As described in the section on homotopy continuation, when finding the steady states of systems with conservation laws, additional considerations have to be taken. E.g. consider the following two-state system:

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*} \]

It has an infinite number of steady states. To make steady state finding possible, information of the system's conserved quantities (here $C = X1 + X2$) must be provided. Since these can be computed from system initial conditions (u0, i.e. those provided when performing ODE simulations), designating an u0 is often the best way. There are two ways to do this. First, one can perform forward ODE simulation-based steady state finding, using the initial condition as the initial u guess. Alternatively, any conserved quantities can be eliminated when the NonlinearProblem is created. This feature is supported by Catalyst's conservation law finding and elimination feature.

To eliminate conservation laws we simply provide the remove_conserved = true argument to NonlinearProblem:

p = [:k1 => 2.0, :k2 => 3.0]
u_guess = [:X1 => 3.0, :X2 => 1.0]
nl_prob = NonlinearProblem(two_state_model, u_guess, p; remove_conserved = true)

┌ Warning: Note, when constructing `NonlinearSystem`s with `remove_conserved = true`, possibly via calling `NonlinearProblem`, `remake(::NonlinearProblem)` has some limitations. If in `remake` the value of the conserved constant is explicitly updated, it is not possible to have one variable's `u0` value auto-update to preserve the conservation law. See the [FAQ entry](https://docs.sciml.ai/Catalyst/stable/faqs/#faq_remake_nonlinprob) for more details. This warning can be disabled by passing `conseqs_remake_warn = false`.
└ @ Catalyst ~/work/Catalyst.jl/Catalyst.jl/src/reactionsystem_conversions.jl:607

here it is important that the quantities used in u_guess correspond to the conserved quantities we wish to use. E.g. here the conserved quantity $X1 + X2 = 3.0 + 1.0 = 4$ holds for the initial condition, and will hence also hold in the computed steady state as well. We can now find the steady states using solve like before:

sol = solve(nl_prob)

retcode: Success
u: 2-element Vector{Float64}:
 2.4
 1.6

We note that the output only provides a single value. The reason is that the actual system solved only contains a single equation (the other being eliminated with the conserved quantity). To find the values of $X1$ and $X2$ we can directly query the solution object for these species' values, using the species themselves as inputs:

sol[[:X1, :X2]]

2-element Vector{Float64}:
 2.4
 1.6

Finding steady states through ODE simulations

The NonlinearProblems generated by Catalyst corresponds to ODEs. A common method of solving these is to simulate the ODE from an initial condition until a steady state is reached. Here we do so for the dimerisation system considered in the previous section. First, we declare our model, initial condition, and parameter values.

dimer_production = @reaction_network begin
    pₘ, 0 --> mRNA
    pₚ, mRNA --> mRNA + P
    (k₁, k₂), 2P <--> P₂
    d, (mRNA, P, P₂) --> 0
end
p = [:pₘ => 0.5, :pₚ => 2.0, :k₁ => 5.0, :k₂ => 1.0, :d => 1.0]
u0 = [:mRNA => 0.1, :P => 0.0, :P₂ => 0.0]

Next, we provide these as an input to a SteadyStateProblem

ssprob = SteadyStateProblem(dimer_production, u0, p)

Finally, we can find the steady states using the solver command. Since SteadyStateProblems are solved through forward ODE simulation, we must load the sublibrary of the OrdinaryDiffEq.jl package that corresponds to the selected ODE solver. Any available ODE solver can be used, however, it has to be encapsulated by the DynamicSS() function. E.g. here we use the Rodas5P solver which is loaded from the OrdinaryDiffEqRosenbrock sublibrary:

(which requires loading the SteadyStateDiffEq package).

using SteadyStateDiffEq, OrdinaryDiffEqRosenbrock
solve(ssprob, DynamicSS(Rodas5P()))

retcode: Success
u: 3-element Vector{Float64}:
 0.49999999427790287
 0.4633249171558076
 0.26833743088545464

Note that, unlike for nonlinear system solving, u0 is not just an initial guess of the solution, but the initial conditions from which the steady state simulation is carried out. This means that, for a system with multiple steady states, we can determine the steady states associated with specific initial conditions (which is not possible when the nonlinear solving approach is used). This also permits us to easily handle the presence of conservation laws. The forward ODE simulation approach (unlike homotopy continuation and nonlinear solving) cannot find unstable steady states.

Generally, SteadyStateProblems can be solved using the same options that are available for ODE simulations. E.g. here we designate a specific dt step size:

solve(ssprob, DynamicSS(Rodas5P()); dt = 0.01)

It is possible to use solve SteadyStateProblems using a nonlinear solver, and NonlinearProblems using forward ODE simulation solvers:

using NonlinearSolve
solve(ssprob, TrustRegion())

nlprob = NonlinearProblem(dimer_production, u0, p)
solve(nlprob, DynamicSS(Rodas5P()))

However, especially when the forward ODE simulation approach is used, it is recommended to use the problem type which corresponds to the intended solver.

Citations

If you use this functionality in your research, in addition to Catalyst, please cite the following paper to support the authors of the NonlinearSolve.jl package:

@article{pal2024nonlinearsolve,
  title={NonlinearSolve. jl: High-Performance and Robust Solvers for Systems of Nonlinear Equations in Julia},
  author={Pal, Avik and Holtorf, Flemming and Larsson, Axel and Loman, Torkel and Schaefer, Frank and Qu, Qingyu and Edelman, Alan and Rackauckas, Chris and others},
  journal={arXiv preprint arXiv:2403.16341},
  year={2024}
}

References

1 J. Nocedal, S. J. Wright, Numerical Optimization, Springer (2006).