Merge pull request #683 from SciML/hc_extension

Homotopy Continuation extension and remake

HISTORY.md

@@ -1,6 +1,17 @@
 # Breaking updates and feature summaries across releases
 
 ## Catalyst unreleased (master branch)
+- Added a CatalystHomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds the steady states of a reactin system using the homotopy continuation method. This feature is only available for julia versions 1.9+. Example: 
+```julia
+wilhelm_2009_model = @reaction_network begin
+    k1, Y --> 2X
+    k2, 2X --> X + Y
+    k3, X + Y --> Y
+    k4, X --> 0
+end
+ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5]
+hc_steady_states(wilhelm_2009_model, ps)
+```
 
 ## Catalyst 13.4
 - Added the ability to create species that represent chemical compounds and know

Project.toml

@@ -21,6 +21,12 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
+[weakdeps]
+HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
+
+[extensions]
+CatalystHomotopyContinuationExtension = "HomotopyContinuation"
+
 [compat]
 DataStructures = "0.18"
 DiffEqBase = "6.83.0"
@@ -42,6 +48,7 @@ julia = "1.6"
 [extras]
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 Graphviz_jll = "3c863552-8265-54e4-a6dc-903eb78fde85"
+HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
@@ -57,4 +64,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["DomainSets", "Graphviz_jll", "NonlinearSolve", "OrdinaryDiffEq", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful"]
+test = ["DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful"]

docs/src/catalyst_applications/homotopy_continuation.md

@@ -8,136 +8,77 @@ example those which are purely mass action or have only have Hill functions with
 integer Hill exponents). The roots of these can reliably be found through a
 *homotopy continuation* algorithm. This is implemented in Julia through the
 [HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org/) package.
-In this tutorial, we will demonstrate how homotopy continuation can be used to
-find the steady states of mass action chemical reaction networks implemented in
-Catalyst.
+
+Catalyst contains a special homotopy continuation extension that is loaded whenever HomotopyContinuation.jl is. This exports a single function, `hc_steady_states`, that can be used to find the steady states of any `ReactionSystem` structure. If you use this in your research, please [cite the HomotopyContinuation.jl](@ref homotopy_continuation_citation) and [Catalyst.jl]() publications.
 
 ## Basic example
 For this tutorial, we will use a model from Wilhem (2009)[^1] (which
 demonstrates bistability in a small chemical reaction network). We declare the
 model and the parameter set for which we want to find the steady states:
 ```@example hc1
-using Catalyst, ModelingToolkit
+using Catalyst
 import HomotopyContinuation
-const MT = ModelingToolkit
-const HC = HomotopyContinuation
 
 wilhelm_2009_model = @reaction_network begin
     k1, Y --> 2X
     k2, 2X --> X + Y
     k3, X + Y --> Y
     k4, X --> 0
 end
-
-# add default parameters values to model
-setdefaults!(wilhelm_2009_model, [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5])
+ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5]
 nothing   # hide
 ```
-Next, we will need to extract the actual equations from our model. In addition,
-we will substitute in our parameter values to these equations.
-```@example hc1
-ns = convert(NonlinearSystem, wilhelm_2009_model)
-
-# this gets the parameter values ordered consistent with parameters(ns)
-pvals = MT.varmap_to_vars([], MT.parameters(ns); defaults = MT.defaults(ns))
+Here, we only run `import HomotopyContinuation` as we do not require any of its functions, and just need it to be present in the current scope for the extension to be activated.
 
-subs = Dict(MT.parameters(ns) .=> pvals)
-neweqs = map(eq -> substitute(eq.rhs, subs), equations(ns))
-```
-Finally, we use Catalyst's `to_multivariate_poly` function to reinterpret our
-symbolic equations in a polynomial representation that is compatible with
-HomotopyContinuation. We can then apply HomotopyContinuation's `solve` command
-to find the roots, using `real_solutions` to filter our non-physical complex
-steady-states:
+Now we can find the steady states using:
 ```@example hc1
-polyeqs = Catalyst.to_multivariate_poly(neweqs)
-sols = HC.real_solutions(HC.solve(polyeqs))
+hc_steady_states(wilhelm_2009_model, ps)
 ```
-Note that HomotopyContinuation orders variables lexicographically, so this will
-be the ordering present in each steady-state solution vector (i.e. `[X1, X2]` is
-the ordering here).
+The order of the species in the output vectors are the same as in `species(wilhelm_2009_model)`.
+
+It should be noted that the steady state multivariate polynomials produced by reaction systems may have both imaginary and negative roots, which are filtered away by `hc_steady_states`. If you want the negative roots, you can use the `hc_steady_states(wilhelm_2009_model, ps; filter_negative=false)` argument.
 
-While it is not the case for this CRN, we note that solutions with negative
-species concentrations can be valid (unphysical) steady-states for certain
-systems. These will need to be filtered out as well.
 
 ## Systems with conservation laws
-Finally, some systems are under-determined, and have an infinite number of
-possible steady states. These are typically systems containing a conservation
+Some systems are under-determined, and have an infinite number of possible steady states. These are typically systems containing a conservation
 law, e.g.
 ```@example hc3
 using Catalyst
 import HomotopyContinuation
-const MT = ModelingToolkit
-const HC = HomotopyContinuation
 
 two_state_model = @reaction_network begin
     (k1,k2), X1 <--> X2
 end
 ```
-Catalyst allows the conservation laws to be computed using the
-`conservationlaws` function. By using these to reduce the dimensionality of the
-system, as well specifying the initial amount of each species,
-HomotopyContinuation can again be used to find steady states. First, we set the
-default values of the system's initial conditions and parameter values. This
-will allow the system to automatically find the conserved amounts.
-```@example hc3
-setdefaults!(two_state_model, [:X1 => 1.0, :X2 => 1.0, :k1 => 2.0, :k2 => 1.0])
-```
-Next, we create a `NonlinearSystem`, while also removing one species via the
-conservation equation.
-```@example hc3
-ns = convert(NonlinearSystem, two_state_model; remove_conserved = true)
-```
-Again, we next create a dictionary for parameter values that we substitute in to
-give our final equation.
+Catalyst allows the conservation laws of such systems to be computed using the `conservationlaws` function. By using these to reduce the dimensionality of the system, as well specifying the initial amount of each species, HomotopyContinuation can again be used to find steady states. To find the steady states using the Catalyst interface to HomotopyContinuation, an initial condition must be provided (which is used to compute the system's conserved quantities, in this case `X1+X2`):
 ```@example hc3
-pvals = MT.varmap_to_vars([], MT.parameters(ns); defaults = MT.defaults(ns))
-subs = Dict(MT.parameters(ns) .=> pvals)
-neweqs = map(eq -> substitute(eq.rhs, subs), equations(ns))
-```
-Notice, our equations are just for `X1` as `X2` was eliminated via the conservation law.
+ps = [:k1 => 2.0, :k2 => 1.0]
+u0 = [:X1 => 1.0, :X2 => 1.0]
+hc_steady_states(wilhelm_2009_model, ps; u0)
 
-Finally, we convert to polynomial form and solve for the steady-states
-```@example hc3
-polyeqs = Catalyst.to_multivariate_poly(neweqs)
-sols = HC.real_solutions(HC.solve(polyeqs))
 ```
 
-If we also want the corresponding value for `X2`, we can substitute
-into the equation for it from the conservation laws:
-```@example hc3
-# get the X2 symbolic variable
-@unpack X2 = two_state_model
-
-# get its algebraic formula in terms of X1 and parameters
-ceqs = conservedequations(two_state_model)
-X2eqidx = findfirst(eq -> isequal(eq.lhs, X2), ceqs)
-X2eq = ceqs[X2eqidx].rhs
+## Final notes
+- `hc_steady_states` supports any systems where all rates are systems of rational polynomials (such as Hill functions with integer Hill coefficients).
+- When providing initial conditions to compute conservation laws, values are only required for those species that are part of conserved quantities. If this set of species is unknown, it is recommended to provide initial conditions for all species. 
+- Additional arguments provided to `hc_steady_states` are automatically passed to HomotopyContinuation's `solve` command. Use e.g. `show_progress=false` to disable the progress bar.
+---
 
-# for each SS, set X1's value in the subs map and calculate X2
-@unpack X1 = two_state_model
-X2 = map(sols) do x
-    X1val = x[1]
-    subs[MT.value(X1)] = X1val
-    substitute(X2eq, subs)
-end
+## [Citation](@id homotopy_continuation_citation)
+If you use this functionality in your research, please cite the following paper to support the authors of the HomotopyContinuation package:
 ```
-giving that the steady-state for `X2` is about `1.33333`.
-
-As an alternative, we could have coupled `neweqs` with the conservation law
-relations to have HomotopyContinuation find the steady-states simultaneously:
-```@example hc3
-# move all the terms in the conserved equations to one side
-# and substitute in the parameter values
-subs = Dict(MT.parameters(ns) .=> pvals)
-conservedrelations = map(eq -> substitute(eq.rhs - eq.lhs, subs), ceqs)
-neweqs = vcat(neweqs, conservedrelations)
-
-# calculate the steady-states
-polyeqs = Catalyst.to_multivariate_poly(neweqs)
-sols = HC.real_solutions(HC.solve(polyeqs))
+@inproceedings{HomotopyContinuation.jl,
+  title={{H}omotopy{C}ontinuation.jl: {A} {P}ackage for {H}omotopy {C}ontinuation in {J}ulia},
+  author={Breiding, Paul and Timme, Sascha},
+  booktitle={International Congress on Mathematical Software},
+  pages={458--465},
+  year={2018},
+  organization={Springer}
+}
 ```
----
+
 ## References
-[^1]: [Wilhelm, T. *The smallest chemical reaction system with bistability*, BMC Systems Biology (2009).](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-050Wilhelm-3-90)
+[^1]: [Thomas Wilhelm, *The smallest chemical reaction system with bistability*, BMC Systems Biology (2009).](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90)
+[^2]: [Paul Breiding, Sascha Timme, *HomotopyContinuation.jl: A Package for Homotopy Continuation in Julia*, International Congress on Mathematical Software (2018).](https://link.springer.com/chapter/10.1007/978-3-319-96418-8_54)
+[^3:] [Andrew J Sommese, Charles W Wampler *The Numerical Solution of Systems of Polynomials Arising in Engineering and Science*, World Scientific (2005).](https://www.worldscientific.com/worldscibooks/10.1142/5763#t=aboutBook)
+[^4:] [Daniel J. Bates, Paul Breiding, Tianran Chen, Jonathan D. Hauenstein, Anton Leykin, Frank Sottile, *Numerical Nonlinear Algebra*, arXiv (2023).](https://arxiv.org/abs/2302.08585)

docs/src/index.md

@@ -130,7 +130,7 @@ Slack's](https://julialang.slack.com) \#sciml-bridged and \#sciml-sysbio channel
 For bugs or feature requests [open an
 issue](https://github.com/SciML/Catalyst.jl/issues).
 
-## Supporting and Citing Catalyst.jl
+## [Supporting and Citing Catalyst.jl](@id catalyst_citation)
 The software in this ecosystem was developed as part of academic research. If you would like to help support it,
 please star the repository as such metrics may help us secure funding in the future. If you use Catalyst as part
 of your research, teaching, or other activities, we would be grateful if you could cite our work:

ext/CatalystHomotopyContinuationExtension.jl

@@ -0,0 +1,12 @@
+module CatalystHomotopyContinuationExtension
+
+# Fetch packages.
+using Catalyst
+import ModelingToolkit as MT
+import HomotopyContinuation as HC
+import Symbolics: unwrap, wrap, Rewriters, symtype, issym, istree
+
+# Creates and exports hc_steady_states function.
+include("CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl")
+
+end

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -0,0 +1,127 @@
+### Homotopy Continuation Based Steady State Finding ###
+
+"""
+    hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=typeof(ps)(), kwargs...)
+
+Uses homotopy continuation via HomotopyContinuation.jl to find the steady states of the ODE system corresponding to the provided reaction system.
+
+Arguments:
+- `rs::ReactionSystem`: The reaction system for which we want to find the steady states.
+- `ps`: The parameter values for which we want to find the steady states.
+- `filter_negative=true`: If set to true, solutions with any species concentration <neg_thres is removed from the output.
+- `neg_thres=-1e-20`: Determine the minimum values for which a species concentration is to be considered non-negative. Species concentrations ``> neg_thres`` but `< 0.0` are set to `0.0`.
+- `u0=nothing`: Initial conditions for which we want to find the steady states. For systems with conservation laws this are required to compute conserved quantities. Initial conditions are not required for all species, only those involved in conserved quantities (if this set is unknown, it is recommended to provide initial conditions for all species). 
+- `kwargs...`: any additional arguments (like `show_progress= true`) are passed into HomotopyContinuation.jl's `solve` call. 
+
+Examples
+```@repl
+rs = @reaction_network begin
+    k1, Y --> 2X
+    k2, 2X --> X + Y
+    k3, X + Y --> Y
+    k4, X --> 0
+end
+ps = [:k3 => 1.0, :k2 => 2.0, :k4 => 1.5, :k1=>8.0]
+hc_sol = hc_steady_states(rs, ps)
+```
+gives
+```
+[0.5000000000000002, 2.0000000000000004]
+[0.0, 0.0]
+[4.499999999999999, 5.999999999999999]
+```
+
+Notes:
+```
+  """
+function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=[], kwargs...)
+    ss_poly = steady_state_polynomial(rs, ps, u0)
+    sols = HC.real_solutions(HC.solve(ss_poly; kwargs...))
+    reorder_sols!(sols, ss_poly, rs)
+    return (filter_negative ? filter_negative_f(sols; neg_thres) : sols)
+end
+
+# For a given reaction system, parameter values, and initial conditions, find the polynomial that HC solves to find steady states.
+function steady_state_polynomial(rs::ReactionSystem, ps, u0)
+    ns = convert(NonlinearSystem, rs; remove_conserved = true)
+    pre_varmap = [symmap_to_varmap(rs,u0)..., symmap_to_varmap(rs,ps)...]
+    conservationlaw_errorcheck(rs, pre_varmap)
+    p_vals = ModelingToolkit.varmap_to_vars(pre_varmap, parameters(ns); defaults = ModelingToolkit.defaults(ns))
+    p_dict  = Dict(parameters(ns) .=> p_vals)
+    eqs_pars_funcs = vcat(equations(ns), conservedequations(rs))
+    eqs_funcs = map(eq -> substitute(eq.rhs - eq.lhs, p_dict), eqs_pars_funcs)
+    eqs = [deregister([mm, mmr, hill, hillr, hillar], eq) for eq in eqs_funcs]
+    eqs_intexp = make_int_exps.(eqs)
+    return Catalyst.to_multivariate_poly(remove_denominators.(eqs_intexp))
+end
+
+# If u0s are not given while conservation laws are present, throws an error.
+function conservationlaw_errorcheck(rs, pre_varmap)
+    vars_with_vals = Set(p[1] for p in pre_varmap)
+    any(s -> s in vars_with_vals, species(rs)) && return
+    isempty(conservedequations(rs)) || 
+        error("The system has conservation laws but initial conditions were not provided for some species.")
+end
+
+# Unfolds a function (like mm or hill). 
+function deregister(fs::Vector{T}, expr) where T
+    for f in fs
+        expr = deregister(f, expr) 
+    end
+    return expr
+end
+# Provided by Shashi Gowda.
+deregister(f, expr) = wrap(Rewriters.Postwalk(Rewriters.PassThrough(___deregister(f)))(unwrap(expr)))
+function ___deregister(f)
+    (expr) ->
+    if istree(expr) && operation(expr) == f
+        args = arguments(expr)
+        invoke_with = map(args) do a
+            t = typeof(a)
+            issym(a) || istree(a) ? wrap(a) => symtype(a) : a => typeof(a)
+        end 
+        invoke(f, Tuple{last.(invoke_with)...}, first.(invoke_with)...)
+    end
+end
+
+# Parses and expression and return a version where any exponents that are Float64 (but an int, like 2.0) are turned into Int64s.
+make_int_exps(expr) = wrap(Rewriters.Postwalk(Rewriters.PassThrough(___make_int_exps))(unwrap(expr))).val
+function ___make_int_exps(expr)
+    !istree(expr) && return expr
+    if (operation(expr) == ^) 
+        if isinteger(arguments(expr)[2])
+            return arguments(expr)[1] ^ Int64(arguments(expr)[2])
+        else
+            error("An non integer ($(arguments(expr)[2])) was found as a variable exponent. Non-integer exponents are not supported for homotopy continuation based steady state finding.")
+        end
+    end
+end
+
+# If the input is a fraction, removes the denominator.
+function remove_denominators(expr)
+    s_expr = simplify_fractions(expr)
+    !istree(expr) && return expr
+    if operation(s_expr) == /
+        return remove_denominators(arguments(s_expr)[1])
+    end
+    if operation(s_expr) == +
+        return sum(remove_denominators(arg) for arg in arguments(s_expr))
+    end
+    return s_expr
+end
+
+# HC orders the solution vector according to the lexicographic values of the variable names. This reorders the output according to the species index in the reaction system species vector.
+function reorder_sols!(sols, ss_poly, rs::ReactionSystem)
+    var_names_extended = String.(Symbol.(HC.variables(ss_poly)))
+    var_names = [Symbol(s[1:prevind(s, findlast('_', s))]) for s in var_names_extended]
+    sort_pattern = indexin(MT.getname.(species(rs)), var_names)
+    foreach(sol -> permute!(sol, sort_pattern), sols)
+end
+
+# Filters away solutions with negative species concentrations (and for neg_thres < val < 0.0, sets val=0.0).
+function filter_negative_f(sols; neg_thres=-1e-20)
+    for sol in sols, idx in 1:length(sol)
+        (neg_thres < sol[idx] < 0) && (sol[idx] = 0)
+    end
+    return filter(sol -> all(>=(0), sol), sols)
+end

src/Catalyst.jl

@@ -101,4 +101,10 @@ export @compound
 export components, iscompound, coefficients
 export balance_reaction
 
+### Extensions ###
+
+# HomotopyContinuation
+function hc_steady_states end
+export hc_steady_states
+
 end # module