# Finding steady-states symbolicaly
**URL:** https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545
**Category:** Modelling & Simulations
**Tags:** modelingtoolkit, catalyst
**Created:** [September 18, 2024, 11:01am UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545 "2024-09-18T11:01:58Z")
**Posts on this page:** 18
**Page:** 1
### Author:  [@adhalanay](https://discourse.julialang.org/u/adhalanay)
#### Post date: [September 18, 2024, 11:01am UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/1 "2024-09-18T11:01:58Z")
Maybe the question is trivial, but I cannot find an answer. Using Catalyst.jl I introduced a reaction network `rn` and constructed the corresponding system of ODEs with `convert(ODESystem,rn)`. The system has a number of parameters. My question is:
Is there a possibility to find the steady-states of the system as a function of the parameters?
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [September 18, 2024, 12:31pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/2 "2024-09-18T12:31:28Z")
> [@adhalanay](#):
>
> Is there a possibility to find the steady-states of the system as a function of the parameters?
That’s a root-finding problem, i.e. you have \frac{d\vec{x}}{dt} = \vec{f}(\vec{x}, \vec{p}) and you want to find the \vec{x}(\vec{p}) such that \vec{f}(\vec{x}, \vec{p}) = 0 for given parameters \vec{p}.
Whether it is practical to find _all_ of the steady states depends on what \vec{f}(\vec{x}, \vec{p}) looks like as a function of \vec{x}. To find a single root, especially if you have some rough guess for \vec{x}, you can try Newton’s method or similar (e.g. via NLsolve.jl). If \vec{f} is linear or affine, then you just have a linear-algebra problem. If \vec{f} consists of polynomials, you could use [HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org/).
Once you have found any steady state for any choice of \vec{p}, you can “track” \vec{x} as a function of \vec{p} as the parameters vary. This is called a “numerical continuation” problem, and one package that can do this is BifurcationKit.jl.
---
### Author:  [@adhalanay](https://discourse.julialang.org/u/adhalanay)
#### Post date: [September 18, 2024, 6:32pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/3 "2024-09-18T18:32:24Z")
Well, the system has polynomial right-side and I believe `symbolic_solve()` is able to solve this kind of equations. No my problem is how to extract these polynomials in order to feed them to the solver.
---
### Author:  [@johannesnauta](https://discourse.julialang.org/u/johannesnauta)
#### Post date: [September 18, 2024, 7:40pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/4 "2024-09-18T19:40:33Z")
I’m not sure about the details, but normally equations are extracted from the `ODESystem` with `equations(sys)`, and they support `Symbolics` syntax, so something like
```julia
julia> [eq.rhs for eq in equations(sys)]
```
might do the trick?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 18, 2024, 7:51pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/5 "2024-09-18T19:51:51Z")
`nsys = convert(NonlinearSystem, rn)` gives the nonlinear system for the steady state equations. `Symbolics.symbolic_solve(equations(nsys), unknowns(nsys))` should use the new symbolic solver. You may want to `using Nemo` for the Grobner basis techniques. That would then give the symbolic solutions, if possible. Give it a try and post here how that goes, the symbolic solver is pretty new so it doesn’t have a tutorial yet.
If not possible but polynomial, then HomotopyContinuation.jl is how to do it, and there’s a tutorial for that.
> **[Finding Steady States through Homotopy Continuation · Catalyst.jl](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/homotopy_continuation/)**
>
> Documentation for Catalyst.jl.
It’s just `hc_steady_states`.
For the purely numerical form, I would not recommend NLsolve.jl. There’s a lot more setup for NonlinearSolve.jl and it’s required for correctness in many ways (symbolic interface, etc.) The tutorial on doing this can be found here:
> **[Finding Steady States using NonlinearSolve.jl and SteadyStateDiffEq.jl ·...](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/nonlinear_solve/)**
>
> Documentation for Catalyst.jl.
---
### Author:  [@adhalanay](https://discourse.julialang.org/u/adhalanay)
#### Post date: [September 18, 2024, 8:56pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/6 "2024-09-18T20:56:43Z")
Thank you for your answer. Here is my example:
> using Catalyst
> using DifferentialEquations
> using WGLMakie
> using Latexify
> using Symbolics
> @variables t
> @species S I R
> @parameters k1 k2 k3 k4
> rn = @reaction\_network begin
> k1, S+I → 2I
> k2, I–\> R
> k3, R → S
> k4, S → R
> end
> odesys=convert(ODESystem,rn)
> nsys = convert(NonlinearSystem, rn)
> symbolic\_solve(equations(nsys), unknowns(nsys))
The last command fails with
> KeyError: key S not found Stacktrace: [1] getindex @ ./iddict.jl:108 [inlined] [2] linearity\_1 @ [~/.julia/packages/Symbolics/XnDVB/src/linearity.jl:57](http://localhost:8888/lab/tree/~/.julia/packages/Symbolics/XnDVB/src/linearity.jl#line=56) [inlined] [3] mark\_vars(expr::SymbolicUtils.BasicSymbolic{Real}, vars::OrderedCollections.OrderedSet{Vector{SymbolicUtils.BasicSymbolic{Real}}}) @ Symbolics [~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl:208](http://localhost:8888/lab/tree/~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl#line=207) [4] (::Base.Fix2{typeof(Symbolics.mark\_vars), OrderedCollections.OrderedSet{Vector{SymbolicUtils.BasicSymbolic{Real}}}})(y::SymbolicUtils.BasicSymbolic{Real}) @ Base ./operators.jl:1135 [5] iterate @ ./generator.jl:47 [inlined] [6] \_collect(c::Vector{Any}, itr::Base.Generator{Vector{Any}, Base.Fix2{typeof(Symbolics.mark\_vars), OrderedCollections.OrderedSet{Vector{SymbolicUtils.BasicSymbolic{Real}}}}}, ::Base.EltypeUnknown, isz::Base.HasShape{1}) @ Base ./array.jl:854 [7] collect\_similar @ ./array.jl:763 [inlined] [8] map @ ./abstractarray.jl:3285 [inlined] [9] mark\_vars(expr::SymbolicUtils.BasicSymbolic{Real}, vars::OrderedCollections.OrderedSet{Vector{SymbolicUtils.BasicSymbolic{Real}}}) @ Symbolics [~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl:204](http://localhost:8888/lab/tree/~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl#line=203) [10] (::Base.Fix2{typeof(Symbolics.mark\_vars), OrderedCollections.OrderedSet{Vector{SymbolicUtils.BasicSymbolic{Real}}}})(y::SymbolicUtils.BasicSymbolic{Real}) @ Base ./operators.jl:1135 [11] iterate @ ./generator.jl:47 [inlined] [12] \_collect(c::Vector{Any}, itr::Base.Generator{Vector{Any}, Base.Fix2{typeof(Symbolics.mark\_vars), OrderedCollections.OrderedSet{Vector{SymbolicUtils.BasicSymbolic{Real}}}}}, ::Base.EltypeUnknown, isz::Base.HasShape{1}) @ Base ./array.jl:854 [13] collect\_similar @ ./array.jl:763 [inlined] [14] map @ ./abstractarray.jl:3285 [inlined] [15] mark\_vars(expr::SymbolicUtils.BasicSymbolic{Real}, vars::OrderedCollections.OrderedSet{Vector{SymbolicUtils.BasicSymbolic{Real}}}) @ Symbolics [~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl:204](http://localhost:8888/lab/tree/~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl#line=203) [16] mark\_and\_exponentiate(expr::SymbolicUtils.BasicSymbolic{Real}, vars::OrderedCollections.OrderedSet{Vector{SymbolicUtils.BasicSymbolic{Real}}}) @ Symbolics [~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl:144](http://localhost:8888/lab/tree/~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl#line=143) [17] semipolyform\_terms(expr::SymbolicUtils.BasicSymbolic{Real}, vars::OrderedCollections.OrderedSet{Vector{SymbolicUtils.BasicSymbolic{Real}}}) @ Symbolics [~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl:160](http://localhost:8888/lab/tree/~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl#line=159) [18] semipolynomial\_form(expr::SymbolicUtils.BasicSymbolic{Real}, vars::Vector{Vector{SymbolicUtils.BasicSymbolic{Real}}}, degree::Float64; consts::Bool) @ Symbolics [~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl:264](http://localhost:8888/lab/tree/~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl#line=263) [19] semipolynomial\_form @ [~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl:257](http://localhost:8888/lab/tree/~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl#line=256) [inlined] [20] polynomial\_coeffs(expr::SymbolicUtils.BasicSymbolic{Real}, vars::Vector{Vector{SymbolicUtils.BasicSymbolic{Real}}}) @ Symbolics [~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl:301](http://localhost:8888/lab/tree/~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl#line=300) [21] check\_poly\_inunivar(poly::SymbolicUtils.BasicSymbolic{Real}, var::Vector{SymbolicUtils.BasicSymbolic{Real}}) @ Symbolics [~/.julia/packages/Symbolics/XnDVB/src/solver/solve\_helpers.jl:97](http://localhost:8888/lab/tree/~/.julia/packages/Symbolics/XnDVB/src/solver/solve_helpers.jl#line=96) [22] symbolic\_solve(expr::Vector{Equation}, x::Vector{SymbolicUtils.BasicSymbolic{Real}}; dropmultiplicity::Bool, warns::Bool) @ Symbolics [~/.julia/packages/Symbolics/XnDVB/src/solver/main.jl:168](http://localhost:8888/lab/tree/~/.julia/packages/Symbolics/XnDVB/src/solver/main.jl#line=167) [23] symbolic\_solve(expr::Vector{Equation}, x::Vector{SymbolicUtils.BasicSymbolic{Real}}) @ Symbolics [~/.julia/packages/Symbolics/XnDVB/src/solver/main.jl:145](http://localhost:8888/lab/tree/~/.julia/packages/Symbolics/XnDVB/src/solver/main.jl#line=144) [24] top-level scope @ In[23]:1
---
### Author:  [@sumiya11](https://discourse.julialang.org/u/sumiya11)
#### Post date: [September 18, 2024, 9:36pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/7 "2024-09-18T21:36:42Z")
I guess `symbolic_solve` doesn’t handle time-dependent variables like S(t). We will fix this.
Out of curiousity : your system seems to have infinite number of solutions for almost all values of k1,k2,k3,k4. Is this expected ?
---
### Author:  [@adhalanay](https://discourse.julialang.org/u/adhalanay)
#### Post date: [September 19, 2024, 6:33am UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/8 "2024-09-19T06:33:04Z")
It is just a toy model to see how Catalyst and Symbolics work together.
---
### Author:  [@adhalanay](https://discourse.julialang.org/u/adhalanay)
#### Post date: [September 19, 2024, 6:42am UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/9 "2024-09-19T06:42:28Z")
Also, wouldn’t be possible to turn a time-dependent variable into a time-independent one? Obviously for finding steady-states time is irrelevant
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [September 19, 2024, 2:22pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/10 "2024-09-19T14:22:30Z")
Yes, one could certainly drop the arguments at some point in the symbolic workflow, but right now that isn’t done anywhere. In any case though, it seems like it would make sense to have `symbolic_solve` work with unknowns that are functions of independent variables.
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [September 19, 2024, 5:09pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/11 "2024-09-19T17:09:17Z")
I just wanted to point out you are writing extraneous commands in defining your network. All you really need is:
```julia
using Catalyst
rn = @reaction_network begin
k1, S+I --> 2I
k2, I --> R
k3, R --> S
k4, S --> R
end
```
Catalyst automatically makes substrates/products species and each `k` a parameter for you already. If you need to access them later you can use `rn.k1`, `rn.S` or such, or use
```julia
@unpack S, I, k2, k3 = rn
```
Defining them outside of `@reaction_network` has no impact on `@reaction_network`, which won’t actually see the versions you declared externally (it will create its own versions internally within the macro).
---
### Author:  [@adhalanay](https://discourse.julialang.org/u/adhalanay)
#### Post date: [October 7, 2024, 1:27pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/12 "2024-10-07T13:27:59Z")
Sorry for the late reply. I just copied an example found somewhere on the net.
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [October 7, 2024, 1:30pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/13 "2024-10-07T13:30:59Z")
If it was in the current Catalyst docs please do let us know via an issue and we can update it accordingly.
---
### Author:  [@adhalanay](https://discourse.julialang.org/u/adhalanay)
#### Post date: [October 7, 2024, 1:53pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/14 "2024-10-07T13:53:29Z")
I think it was some older version which came out in duckduckgo.
---
### Author:  [@adhalanay](https://discourse.julialang.org/u/adhalanay)
#### Post date: [October 7, 2024, 1:54pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/15 "2024-10-07T13:54:59Z")
Do you believe that I should file a bug report or feature request and if yes where? To Symbolics, Catalyst or ModellingToolkit?
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [October 7, 2024, 3:00pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/16 "2024-10-07T15:00:58Z")
You could open an issue on Symbolics.jl, but it would be good to make a minimal example that only uses Symbolics, i.e. something like
```julia
@variables t A(t) B(t)
# some basic calculation with symbolic_solve that fails using A and B
```
---
### Author:  [@adhalanay](https://discourse.julialang.org/u/adhalanay)
#### Post date: [October 14, 2024, 12:55pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/17 "2024-10-14T12:55:06Z")
I filed a bug report with [https://juliahub.com/ui/Packages/General/Symbolics](https://juliahub.com/ui/Packages/General/Symbolics). It is [https://github.com/JuliaSymbolics/Symbolics.jl/issues/1301](https://github.com/JuliaSymbolics/Symbolics.jl/issues/1301). So far no activity.
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [October 14, 2024, 3:06pm UTC](https://discourse.julialang.org/t/finding-steady-states-symbolicaly/119545/18 "2024-10-14T15:06:05Z")
Sorry, there isn’t much I can do about that. I’m not involved in Symbolics.jl development. I see someone was assigned to the issue just now, so hopefully they can get to it in the near future.