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?

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.

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.

1 Like

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.

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> [eq.rhs for eq in equations(sys)]
```

might do the trick?

`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.

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:

1 Like

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 [inlined] [3] mark_vars(expr::SymbolicUtils.BasicSymbolic{Real}, vars::OrderedCollections.OrderedSet{Vector{SymbolicUtils.BasicSymbolic{Real}}}) @ Symbolics ~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl:208 [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 [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 [16] mark_and_exponentiate(expr::SymbolicUtils.BasicSymbolic{Real}, vars::OrderedCollections.OrderedSet{Vector{SymbolicUtils.BasicSymbolic{Real}}}) @ Symbolics ~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl:144 [17] semipolyform_terms(expr::SymbolicUtils.BasicSymbolic{Real}, vars::OrderedCollections.OrderedSet{Vector{SymbolicUtils.BasicSymbolic{Real}}}) @ Symbolics ~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl:160 [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 [19] semipolynomial_form @ ~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl:257 [inlined] [20] polynomial_coeffs(expr::SymbolicUtils.BasicSymbolic{Real}, vars::Vector{Vector{SymbolicUtils.BasicSymbolic{Real}}}) @ Symbolics ~/.julia/packages/Symbolics/XnDVB/src/semipoly.jl:301 [21] check_poly_inunivar(poly::SymbolicUtils.BasicSymbolic{Real}, var::Vector{SymbolicUtils.BasicSymbolic{Real}}) @ Symbolics ~/.julia/packages/Symbolics/XnDVB/src/solver/solve_helpers.jl:97 [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 [23] symbolic_solve(expr::Vector{Equation}, x::Vector{SymbolicUtils.BasicSymbolic{Real}}) @ Symbolics ~/.julia/packages/Symbolics/XnDVB/src/solver/main.jl:145 [24] top-level scope @ In[23]:1

1 Like

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 ?

It is just a toy model to see how Catalyst and Symbolics work together.

Also, wouldnâ€™t be possible to turn a time-dependent variable into a time-independent one? Obviously for finding steady-states time is irrelevant

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.

I just wanted to point out you are writing extraneous commands in defining your network. All you really need is:

```
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

```
@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).

Sorry for the late reply. I just copied an example found somewhere on the net.

If it was in the current Catalyst docs please do let us know via an issue and we can update it accordingly.

I think it was some older version which came out in duckduckgo.

Do you believe that I should file a bug report or feature request and if yes where? To Symbolics, Catalyst or ModellingToolkit?

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

```
@variables t A(t) B(t)
# some basic calculation with symbolic_solve that fails using A and B
```

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.