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.
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:
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
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
I filed a bug report with https://juliahub.com/ui/Packages/General/Symbolics. It is https://github.com/JuliaSymbolics/Symbolics.jl/issues/1301. So far no activity.
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.