|
| 1 | +# [Finding Steady States using NonlinearSolve.jl](@id nonlinear_solve) |
| 2 | + |
| 3 | +We have previously described how `ReactionSystem` steady states can be found through [homotopy continuation](@ref homotopy_continuation). Catalyst also supports the creation of `NonlinearProblems` corresponding to `ODEProblems` with all left-hand side derivatives set to 0. These can be solved using the variety of nonlinear equation-solving algorithms implemented by [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl), with the solutions corresponding to system steady states. Generally, this approach is preferred when one only wishes to find *a single solution* (with homotopy continuation being preferred when one wishes to find *all solutions*). |
| 4 | + |
| 5 | +This tutorial describes how to create `NonlinearProblem`s from Catalyst's `ReactionSystemn`s, and how to solve them using NonlinearSolve. More extensive descriptions of available solvers and options can be found in [NonlinearSolve's documentation](https://docs.sciml.ai/NonlinearSolve/stable/). If you use this in your research, please [cite the NonlinearSolve.jl](@ref nonlinear_solve_citation) and [Catalyst.jl](@ref catalyst_citation) publications. |
| 6 | + |
| 7 | +## Basic example |
| 8 | +Let us consider a simple dimerisation network, 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. |
| 9 | +```@example nonlinear_solve1 |
| 10 | +using Catalyst |
| 11 | +dimer_production = @reaction_network begin |
| 12 | + pₘ, 0 --> mRNA |
| 13 | + pₚ, mRNA --> mRNA + P |
| 14 | + (k₁, k₂), 2P <--> P₂ |
| 15 | + d, (mRNA, P, P₂) --> 0 |
| 16 | +end |
| 17 | +``` |
| 18 | +This system corresponds to the following ODE: |
| 19 | +```math |
| 20 | +\begin{aligned} |
| 21 | +\frac{dmRNA}{dt} &= pₘ - d \cdot mRNA \\ |
| 22 | +\frac{dP}{dt} &= pₚ \cdot mRNA - k₁ \cdot P + 2k₂ \cdot P₂ - d \cdot P \\ |
| 23 | +\frac{dP₂}{dt} &= k₁ \cdot P + 2k₂ \cdot P₂ \\ |
| 24 | +\end{aligned} |
| 25 | +``` |
| 26 | +To find its steady states we need to solve: |
| 27 | +```math |
| 28 | +\begin{aligned} |
| 29 | +0 &= pₘ - d \cdot mRNA \\ |
| 30 | +0 &= pₚ \cdot mRNA - k₁ \cdot P + 2k₂ \cdot P₂ - d \cdot P \\ |
| 31 | +0 &= k₁ \cdot P + 2k₂ \cdot P₂ \\ |
| 32 | +\end{aligned} |
| 33 | +``` |
| 34 | + |
| 35 | +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](@ref nonlinear_solve_conservation_laws)). Using these, we can create the `NonlinearProblem` that we wish to solve. |
| 36 | +```@example nonlinear_solve1 |
| 37 | +p = [:pₘ => 0.5, :pₚ => 2.0, :k₁ => 5.0, :k₂ => 1.0, :d => 1.0] |
| 38 | +u_guess = [:mRNA => 1.0, :P => 1.0, :P₂ => 1.0] |
| 39 | +nl_prob = NonlinearProblem(dimer_production, u_guess, p) |
| 40 | +``` |
| 41 | +Finally, we can solve it using the `solve` command, returning the steady state solution: |
| 42 | +```@example nonlinear_solve1 |
| 43 | +using NonlinearSolve |
| 44 | +solve(nl_prob) |
| 45 | +``` |
| 46 | + |
| 47 | +NonlinearSolve provides [a wide range of potential solvers](https://docs.sciml.ai/NonlinearSolve/stable/solvers/NonlinearSystemSolvers/). If we wish to designate one, it can be supplied as a second argument to `solve`. Here, we use the Newton Trust Region method: |
| 48 | +```@example nonlinear_solve1 |
| 49 | +solve(nl_prob, TrustRegion()) |
| 50 | +nothing # hide |
| 51 | +``` |
| 52 | + |
| 53 | +## [Finding steady states through ODE simulations](@id nonlinear_solve_ode_simulation_based) |
| 54 | +The `NonlinearProblem`s 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. NonlinearSolve supports this through the `DynamicSS` method. To use it, an appropriate ODE solver (and any options you wish it to use) must also be supplied (with [a large number being available](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)). Here, we will use the `Tsit5` ODE solver to find the steady states of our dimerisation system. |
| 55 | +```@example nonlinear_solve1 |
| 56 | +using OrdinaryDiffEq, SteadyStateDiffEq |
| 57 | +solve(nprob, DynamicSS(Tsit5())) |
| 58 | +``` |
| 59 | +Here, we had to import [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) (to use `Tsit5`) and [SteadyStateDiffEq.jl](https://github.com/SciML/SteadyStateDiffEq.jl) (to use `DynamicSS`). |
| 60 | + |
| 61 | +Like previously, the found steady state is unaffected by the initial `u` guess. However, when [finding the steady states of systems with conservation laws](@ref nonlinear_solve_conservation_laws) it is important to select a `u` guess corresponding to the initial condition of the ODE simulation which steady state is sought. |
| 62 | + |
| 63 | + |
| 64 | +## [Systems with conservation laws](@id nonlinear_solve_conservation_laws) |
| 65 | +As described in the section on homotopy continuation, when finding the steady states of systems with conservation laws, [additional considerations have to be taken](homotopy_continuation_conservation_laws). E.g. consider the following two-state system: |
| 66 | +```@example nonlinear_solve2 |
| 67 | +using Catalyst, NonlinearSolve # hide |
| 68 | +two_state_model = @reaction_network begin |
| 69 | + (k1,k2), X1 <--> X2 |
| 70 | +end |
| 71 | +``` |
| 72 | +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. This can either be done by performing [ODE simulation-based steady state finding](@ref), using the initial condition as the initial `u` guess. Alternatively, any conservation quantities can be eliminated when the `NonlinearProblem` is created. This feature is supported by Catalyst's [conservation law finding and elimination feature](@ref network_analysis_deficiency). |
| 73 | + |
| 74 | +To eliminate conservation laws we simply provide the `remove_conserved = true` argument to `NonlinearProblem`: |
| 75 | +```@example nonlinear_solve2 |
| 76 | +p = [:k1 => 2.0, :k2 => 3.0] |
| 77 | +u_guess = [:X1 => 3.0, :X2 => 1.0] |
| 78 | +nl_prob = NonlinearProblem(two_state_model, u_guess, p; remove_conserved = true) |
| 79 | +nothing # hide |
| 80 | +``` |
| 81 | +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 in the input, and will hence do in the output as well. We can now find the steady states using `solve` like before: |
| 82 | +```@example nonlinear_solve2 |
| 83 | +sol = solve(nl_prob) |
| 84 | +``` |
| 85 | +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 [inference the solution object using these as input](@ref simulation_structure_interfacing_solutions): |
| 86 | +```@example nonlinear_solve2 |
| 87 | +@unpack X1, X2 = two_state_model |
| 88 | +sol[X1] |
| 89 | +``` |
| 90 | +```@example nonlinear_solve2 |
| 91 | +sol[X2] |
| 92 | +``` |
| 93 | + |
| 94 | +--- |
| 95 | +## [Citations](@id nonlinear_solve_citations) |
| 96 | +If you use this functionality in your research, [in addition to Catalyst](@ref catalyst_citation), please cite the following papers to support the authors of the NonlinearSolve.jl package: |
| 97 | +``` |
| 98 | +A NonlinearSolve. jl-related publication is in preparation, once it is available, its details will be added here. |
| 99 | +``` |
0 commit comments