From 21882a0c5b798c0fc033999afeb2f99fbc91bbe1 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 19 Sep 2023 16:01:07 -0400 Subject: [PATCH 01/30] Initiate --- Project.toml | 6 ++ ext/HomotopyContinuationExtension.jl | 11 +++ .../homotopy_continuation_extension.jl | 69 +++++++++++++++++ src/Catalyst.jl | 6 ++ test/extensions/homotopy_continuation.jl | 74 +++++++++++++++++++ test/runtests.jl | 3 + 6 files changed, 169 insertions(+) create mode 100644 ext/HomotopyContinuationExtension.jl create mode 100644 ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl create mode 100644 test/extensions/homotopy_continuation.jl diff --git a/Project.toml b/Project.toml index e1198f1ef7..701ac12ffa 100644 --- a/Project.toml +++ b/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] +HomotopyContinuationExtension = "HomotopyContinuation" + [compat] DataStructures = "0.18" DiffEqBase = "6.83.0" diff --git a/ext/HomotopyContinuationExtension.jl b/ext/HomotopyContinuationExtension.jl new file mode 100644 index 0000000000..b06fb83220 --- /dev/null +++ b/ext/HomotopyContinuationExtension.jl @@ -0,0 +1,11 @@ +module HomotopyContinuationExtension + +# Fetch packages. +using Catalyst +import ModelingToolkit as MT +import HomotopyContinuation as HC + +# Creates and exports hc_steady_states function. +include("HomotopyContinuationExtension/homotopy_continuation_extension.jl") + +end \ No newline at end of file diff --git a/ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl new file mode 100644 index 0000000000..e1e57810ec --- /dev/null +++ b/ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -0,0 +1,69 @@ +### 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 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 <0 is removed from the output. +- `neg_thres=-1e-20`: Determine the minimum values for which a species concentration is to be considred non-negative. Species conentrations ``> neg_thres`` but `< 0.0` are set to `0.0`. +- `u0=typeof(ps)()`: Initial conditions for which we want to find the steady states. For systems with conservation laws this are required to compute conserved quantities. +- `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: +- This is a wrapper around the `solve` function provided by HomotopyContinuation.jl, all credit for this functionality to that package's authors. +``` + """ +function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=typeof(ps)(), 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=neg_thres) : sols) +end + +# For a given reaction system, paraemter 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 = map(pair -> pair::Pair{Num,Float64}, [symmap_to_varmap(rs,u0); symmap_to_varmap(rs,ps)]) # Needed in case ps or u0 are empty (making the combination a Vector{Any}, causing an error). + p_vals = MT.varmap_to_vars(pre_varmap, parameters(ns); defaults = MT.defaults(ns)) + p_dict = Dict(parameters(ns) .=> p_vals) + eqs = vcat(equations(ns), conservedequations(rs)) + eqs = map(eq -> substitute(eq.rhs - eq.lhs, p_dict), eqs) + return Catalyst.to_multivariate_poly(eqs) +end + +# HC orders the solution vector according to the lexiographic values of the variable names. This reorders the output acording 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.0) && (sol[idx] = 0.0) + end + return filter(sol -> all(>=(0.0), sol), sols) +end \ No newline at end of file diff --git a/src/Catalyst.jl b/src/Catalyst.jl index a8763c3044..b33b6782d0 100644 --- a/src/Catalyst.jl +++ b/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 diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl new file mode 100644 index 0000000000..a2ab3785ca --- /dev/null +++ b/test/extensions/homotopy_continuation.jl @@ -0,0 +1,74 @@ +### Fetch Packages ### +using Catalyst, OrdinaryDiffEq, Test +import HomotopyContinuation + +### Run Tests ### + +# Tests for network without conservation laws. +# Tests for Symbol parameter input. +# Tests for Symbolics initial condiiton input. +# Tests for different types (Symbol/Symbolics) for parameters and initial conditions. +let + rs = @reaction_network begin + (k1,k2), X1 <--> X2 + (k3,k4), 2X2 + X3 <--> X2_2X3 + end + @unpack k1, k2, k3, k4 = rs + ps = [k1 => 1.0, k2 => 2.0, k3 => 2.0, k4 => 2.0] + u0 = [:X1 => 2.0, :X2 => 2.0, :X3 => 2.0, :X2_2X3 => 2.0] + + sim_ss = solve(ODEProblem(rs, u0, (0.0,1000.0), ps), Tsit5(); abstol=1e-12, reltol=1e-12)[end] + hc_ss = hc_steady_states(rs, ps; u0=u0)[1] + @test sim_ss ≈ hc_ss +end + +# Tests for network with multiple steady state. +# Tests for Symbol parameter input. +# Tests tha passing kwargs to HC.solve does not error. +let + wilhelm_2009_model = @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_ss_1 = hc_steady_states(wilhelm_2009_model, ps, seed=0x000004d1) + @test sort(hc_ss_1, by=sol->sol[1]) ≈ [[0.0, 0.0], [0.5, 2.0], [4.5, 6.0]] + + hc_ss_2 = hc_steady_states(wilhelm_2009_model, ps, seed=0x000004d2) + hc_ss_3 = hc_steady_states(wilhelm_2009_model, ps, seed=0x000004d2) + @test hc_ss_1 != hc_ss_2 + @test hc_ss_2 == hc_ss_3 +end + +# Tests that reordering is correct. +# Tests corectness in presence of default values. +# Tests where some defaul values are overwritten with other values +let + rs_1 = @reaction_network begin + @parameters kX1=1.0 kX2=2.0 kY1=12345.0 + @species X1(t)=0.1 X2(t)=0.2 Y1(t)=12345.0 + (kX1,kX2), X1 <--> X2 + (kY1,kY2), Y1 <--> Y2 + (kZ1,kZ2), Z1 <--> Z2 + end + ps = [:kY1 => 1.0, :kY2 => 3.0, :kZ1 => 1.0, :kZ2 => 4.0] + u0_1 = [:Y1 => 1.0, :Y2 => 3.0, :Z1 => 10.0, :Z2 =>40.0] + + ss_1 = sort(hc_steady_states(rs_1, ps; u0=u0_1), by=sol->sol[1]) + @test ss_1 ≈ [[0.2, 0.1, 3.0, 1.0, 40.0, 10.0]] + + rs_2 = @reaction_network begin + @parameters kX1=1.0 kX2=2.0 kY1=12345.0 + @species C2(t)=0.1 C1(t)=0.2 B2(t)=12345.0 + (kX1,kX2), C2 <--> C1 + (kY1,kY2), B2 <--> B1 + (kZ1,kZ2), A2 <--> A1 + end + u0_2 = [:B2 => 1.0, :B1 => 3.0, :A2 => 10.0, :A1 =>40.0] + + ss_2 = sort(hc_steady_states(rs_2, ps; u0=u0_2), by=sol->sol[1]) + @test ss_1 ≈ ss_2 +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 5a772722f2..cc7a1135c1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -46,4 +46,7 @@ using SafeTestsets if !Sys.isapple() @time @safetestset "Graphs" begin include("visualization/graphs.jl") end end + + ### Tests extensions. ### + @time @safetestset "Homotopy Continuation Extension" begin include("extensions/homotopy_continuation.jl") end end # @time From d598876899f5967958978212e0640b24b8b84d77 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 19 Sep 2023 16:21:02 -0400 Subject: [PATCH 02/30] Error check when no u0 given when conservation laws exist --- .../homotopy_continuation_extension.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl index e1e57810ec..f6010538ec 100644 --- a/ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -29,6 +29,7 @@ gives [0.5000000000000002, 2.0000000000000004] [0.0, 0.0] [4.499999999999999, 5.999999999999999] +``` Notes: - This is a wrapper around the `solve` function provided by HomotopyContinuation.jl, all credit for this functionality to that package's authors. @@ -45,6 +46,7 @@ end function steady_state_polynomial(rs::ReactionSystem, ps, u0) ns = convert(NonlinearSystem, rs; remove_conserved = true) pre_varmap = map(pair -> pair::Pair{Num,Float64}, [symmap_to_varmap(rs,u0); symmap_to_varmap(rs,ps)]) # Needed in case ps or u0 are empty (making the combination a Vector{Any}, causing an error). + conservationlaw_errorcheck(rs, pre_varmap) p_vals = MT.varmap_to_vars(pre_varmap, parameters(ns); defaults = MT.defaults(ns)) p_dict = Dict(parameters(ns) .=> p_vals) eqs = vcat(equations(ns), conservedequations(rs)) @@ -52,6 +54,14 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0) return Catalyst.to_multivariate_poly(eqs) end +# If u0s are not given while conservation laws are present, throws an error. +function conservationlaw_errorcheck(rs, pre_varmap) + vars_with_vals = union(first.(pre_varmap), keys(ModelingToolkit.defaults(rs))) + isempty(intersect(species(rs), vars_with_vals)) || return + isempty(conservedequations(rs)) && return + error("The system have conservation laws but no initial conditions were provided. Please provide initial conditions.") +end + # HC orders the solution vector according to the lexiographic values of the variable names. This reorders the output acording 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))) From 67a1e7824633ae0f9117f64ea438cabe1f1b05b7 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 19 Sep 2023 21:16:53 -0400 Subject: [PATCH 03/30] enable rational function and parametric exponents. --- Project.toml | 3 +- ext/HomotopyContinuationExtension.jl | 1 + .../homotopy_continuation_extension.jl | 57 +++++++++++++++++-- src/networkapi.jl | 7 +-- test/extensions/homotopy_continuation.jl | 39 ++++++++++--- test/runtests.jl | 3 + 6 files changed, 90 insertions(+), 20 deletions(-) diff --git a/Project.toml b/Project.toml index 701ac12ffa..285968cbae 100644 --- a/Project.toml +++ b/Project.toml @@ -48,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" @@ -63,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"] diff --git a/ext/HomotopyContinuationExtension.jl b/ext/HomotopyContinuationExtension.jl index b06fb83220..fa65c6e481 100644 --- a/ext/HomotopyContinuationExtension.jl +++ b/ext/HomotopyContinuationExtension.jl @@ -4,6 +4,7 @@ module HomotopyContinuationExtension 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("HomotopyContinuationExtension/homotopy_continuation_extension.jl") diff --git a/ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl index f6010538ec..53eb0258bb 100644 --- a/ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -45,13 +45,15 @@ end # For a given reaction system, paraemter 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 = map(pair -> pair::Pair{Num,Float64}, [symmap_to_varmap(rs,u0); symmap_to_varmap(rs,ps)]) # Needed in case ps or u0 are empty (making the combination a Vector{Any}, causing an error). + pre_varmap = map(pair -> pair::Pair{Num,Float64}, [symmap_to_varmap(rs,u0); symmap_to_varmap(rs,ps)]) conservationlaw_errorcheck(rs, pre_varmap) - p_vals = MT.varmap_to_vars(pre_varmap, parameters(ns); defaults = MT.defaults(ns)) - p_dict = Dict(parameters(ns) .=> p_vals) - eqs = vcat(equations(ns), conservedequations(rs)) - eqs = map(eq -> substitute(eq.rhs - eq.lhs, p_dict), eqs) - return Catalyst.to_multivariate_poly(eqs) + 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. @@ -62,6 +64,49 @@ function conservationlaw_errorcheck(rs, pre_varmap) error("The system have conservation laws but no initial conditions were provided. Please provide initial conditions.") 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) == ^) && isinteger(arguments(expr)[2]) + return arguments(expr)[1] ^ Int64(arguments(expr)[2]) + 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 +(remove_denominators.(arguments(s_expr))...) + end + return s_expr +end + # HC orders the solution vector according to the lexiographic values of the variable names. This reorders the output acording 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))) diff --git a/src/networkapi.jl b/src/networkapi.jl index 6d3613292b..7c8b3e6d1d 100644 --- a/src/networkapi.jl +++ b/src/networkapi.jl @@ -399,6 +399,7 @@ function netstoichmat(rn::ReactionSystem; sparse = false) end # the following function is adapted from SymbolicUtils.jl v.19 +# later on (Spetember 2023) modified by Torkel and Shashi (now assumes input not on polynomial form, which is handled elsewhere, previous version failed in these cases anyway). # Copyright (c) 2020: Shashi Gowda, Yingbo Ma, Mason Protter, Julia Computing. # MIT license """ @@ -407,15 +408,13 @@ end Convert the given system of polynomial equations to multivariate polynomial representation. For example, this can be used in HomotopyContinuation.jl functions. """ -function to_multivariate_poly(polyeqs::AbstractVector{BasicSymbolic{Real}}) +function to_multivariate_poly(polyeqs::AbstractVector{Symbolics.BasicSymbolic{Real}}) @assert length(polyeqs)>=1 "At least one expression must be passed to `multivariate_poly`." pvar2sym, sym2term = SymbolicUtils.get_pvar2sym(), SymbolicUtils.get_sym2term() ps = map(polyeqs) do x if istree(x) && operation(x) == (/) - num, den = arguments(x) - PolyForm(num, pvar2sym, sym2term).p / - PolyForm(den, pvar2sym, sym2term).p + error("We should not be able to get here, please contact the package authors.") else PolyForm(x, pvar2sym, sym2term).p end diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index a2ab3785ca..cbcfb6b109 100644 --- a/test/extensions/homotopy_continuation.jl +++ b/test/extensions/homotopy_continuation.jl @@ -1,5 +1,5 @@ ### Fetch Packages ### -using Catalyst, OrdinaryDiffEq, Test +using Catalyst, Test import HomotopyContinuation ### Run Tests ### @@ -8,6 +8,7 @@ import HomotopyContinuation # Tests for Symbol parameter input. # Tests for Symbolics initial condiiton input. # Tests for different types (Symbol/Symbolics) for parameters and initial conditions. +# Tests that attempts to find steady states of system with conservation laws, while u0 is not provided, gives an error. let rs = @reaction_network begin (k1,k2), X1 <--> X2 @@ -17,9 +18,11 @@ let ps = [k1 => 1.0, k2 => 2.0, k3 => 2.0, k4 => 2.0] u0 = [:X1 => 2.0, :X2 => 2.0, :X3 => 2.0, :X2_2X3 => 2.0] - sim_ss = solve(ODEProblem(rs, u0, (0.0,1000.0), ps), Tsit5(); abstol=1e-12, reltol=1e-12)[end] - hc_ss = hc_steady_states(rs, ps; u0=u0)[1] - @test sim_ss ≈ hc_ss + hc_ss = hc_steady_states(rs, ps; u0=u0, show_progress=false)[1] + f = ODEFunction(convert(ODESystem, rs)) + @test f(hc_ss, last.(ps), 0.0)[1] == 0.0 + + @test_throws Exception hc_steady_states(rs, ps; show_progress=false) end # Tests for network with multiple steady state. @@ -34,11 +37,11 @@ let end ps = [:k3 => 1.0, :k2 => 2.0, :k4 => 1.5, :k1=>8.0] - hc_ss_1 = hc_steady_states(wilhelm_2009_model, ps, seed=0x000004d1) + hc_ss_1 = hc_steady_states(wilhelm_2009_model, ps; seed=0x000004d1, show_progress=false) @test sort(hc_ss_1, by=sol->sol[1]) ≈ [[0.0, 0.0], [0.5, 2.0], [4.5, 6.0]] - hc_ss_2 = hc_steady_states(wilhelm_2009_model, ps, seed=0x000004d2) - hc_ss_3 = hc_steady_states(wilhelm_2009_model, ps, seed=0x000004d2) + hc_ss_2 = hc_steady_states(wilhelm_2009_model, ps; seed=0x000004d2, show_progress=false) + hc_ss_3 = hc_steady_states(wilhelm_2009_model, ps; seed=0x000004d2, show_progress=false) @test hc_ss_1 != hc_ss_2 @test hc_ss_2 == hc_ss_3 end @@ -57,7 +60,7 @@ let ps = [:kY1 => 1.0, :kY2 => 3.0, :kZ1 => 1.0, :kZ2 => 4.0] u0_1 = [:Y1 => 1.0, :Y2 => 3.0, :Z1 => 10.0, :Z2 =>40.0] - ss_1 = sort(hc_steady_states(rs_1, ps; u0=u0_1), by=sol->sol[1]) + ss_1 = sort(hc_steady_states(rs_1, ps; u0=u0_1, show_progress=false), by=sol->sol[1]) @test ss_1 ≈ [[0.2, 0.1, 3.0, 1.0, 40.0, 10.0]] rs_2 = @reaction_network begin @@ -69,6 +72,24 @@ let end u0_2 = [:B2 => 1.0, :B1 => 3.0, :A2 => 10.0, :A1 =>40.0] - ss_2 = sort(hc_steady_states(rs_2, ps; u0=u0_2), by=sol->sol[1]) + ss_2 = sort(hc_steady_states(rs_2, ps; u0=u0_2, show_progress=false), by=sol->sol[1]) @test ss_1 ≈ ss_2 +end + +# Tests that non-scalar reaction rates work. +# Tests that rational polynomial steady state systems work. +# Test filter_negative=false works. +let + rs = @reaction_network begin + 0.1 + hill(X,v,K,n), 0 --> X + d, X --> 0 + end + ps = [:v => 5.0, :K => 2.5, :n => 3, :d => 1.0] + sss = hc_steady_states(rs, ps; filter_negative=false, show_progress=false) + + f = ODEFunction(convert(ODESystem, rs)) + @test length(sss) == 4 + for ss in sss + @test f(sss[1], last.(ps), 0.0)[1] == 0.0 + end end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index cc7a1135c1..88532a8b97 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,9 @@ using SafeTestsets ### Run the tests ### @time begin + + ### Tests extensions. ### + @time @safetestset "Homotopy Continuation Extension" begin include("extensions/homotopy_continuation.jl") end ### Tests the properties of ReactionSystems. ### @time @safetestset "ReactionSystem" begin include("reactionsystem_structure/reactionsystem.jl") end From 6bf65cc157968f818959b35f5826c065c22e7e7f Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 19 Sep 2023 21:18:25 -0400 Subject: [PATCH 04/30] testupdate --- test/runtests.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 88532a8b97..cc7a1135c1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,9 +3,6 @@ using SafeTestsets ### Run the tests ### @time begin - - ### Tests extensions. ### - @time @safetestset "Homotopy Continuation Extension" begin include("extensions/homotopy_continuation.jl") end ### Tests the properties of ReactionSystems. ### @time @safetestset "ReactionSystem" begin include("reactionsystem_structure/reactionsystem.jl") end From ddb8696bfc868e61b4c8a2d8df96646bbb19e388 Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 20 Sep 2023 11:03:36 -0400 Subject: [PATCH 05/30] Update history file. --- HISTORY.md | 1 + 1 file changed, 1 insertion(+) diff --git a/HISTORY.md b/HISTORY.md index b3bc6d0366..521518b82d 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -1,6 +1,7 @@ # Breaking updates and feature summaries across releases ## Catalyst unreleased (master branch) +- Added a HomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds teh steady states of a reactin system using the homotopy continuation method. ## Catalyst 13.4 - Added the ability to create species that represent chemical compounds and know From e8b2e67012aa2c2a8318e07aa50cf6fc4ab3a943 Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 20 Sep 2023 11:52:28 -0400 Subject: [PATCH 06/30] doc update --- .../homotopy_continuation.md | 125 +++++------------- 1 file changed, 31 insertions(+), 94 deletions(-) diff --git a/docs/src/catalyst_applications/homotopy_continuation.md b/docs/src/catalyst_applications/homotopy_continuation.md index 4d04108de3..30898f2291 100644 --- a/docs/src/catalyst_applications/homotopy_continuation.md +++ b/docs/src/catalyst_applications/homotopy_continuation.md @@ -8,19 +8,16 @@ 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. ## 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 @@ -28,116 +25,56 @@ wilhelm_2009_model = @reaction_network begin 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. +To find the steady states of these, an initial condition must also be provided (which is used to compute the system's conserved quantities, in this case `X1+X2`): ```@example hc3 -ns = convert(NonlinearSystem, two_state_model; remove_conserved = true) +ps = [:k1 => 2.0, :k2 => 1.0] +u0 = [:X1 => 1.0, :X2 => 1.0] +hc_steady_states(wilhelm_2009_model, ps; u0=u0) ``` -Again, we next create a dictionary for parameter values that we substitute in to -give our final equation. -```@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. -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 polynomial (such as hill functions with integer hill coefficients). +- 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 +## Citations +If you use this functionality in your reasearch, 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) +[^2]: [Breiding, P., Timme, S. *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) From 6a045812a4906a28c367622f5c6b9229d53bd8d5 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 25 Sep 2023 08:55:45 -0400 Subject: [PATCH 07/30] update references --- docs/src/catalyst_applications/homotopy_continuation.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/src/catalyst_applications/homotopy_continuation.md b/docs/src/catalyst_applications/homotopy_continuation.md index 30898f2291..c02d76ebcd 100644 --- a/docs/src/catalyst_applications/homotopy_continuation.md +++ b/docs/src/catalyst_applications/homotopy_continuation.md @@ -76,5 +76,7 @@ If you use this functionality in your reasearch, please cite the following paper ``` ## 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) -[^2]: [Breiding, P., Timme, S. *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) +[^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) \ No newline at end of file From c3cbd96e0b6e539257c9c8678be2f01bfd2d9058 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 25 Sep 2023 12:11:04 -0400 Subject: [PATCH 08/30] Improve extension name (to avoid conflicts) --- Project.toml | 2 +- .../homotopy_continuation_extension.jl | 0 ext/HomotopyContinuationExtension.jl | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) rename ext/{HomotopyContinuationExtension => CatalystHomotopyContinuationExtension}/homotopy_continuation_extension.jl (100%) diff --git a/Project.toml b/Project.toml index 285968cbae..66e1708d97 100644 --- a/Project.toml +++ b/Project.toml @@ -25,7 +25,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" [extensions] -HomotopyContinuationExtension = "HomotopyContinuation" +CatalystHomotopyContinuationExtension = "HomotopyContinuation" [compat] DataStructures = "0.18" diff --git a/ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl similarity index 100% rename from ext/HomotopyContinuationExtension/homotopy_continuation_extension.jl rename to ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl diff --git a/ext/HomotopyContinuationExtension.jl b/ext/HomotopyContinuationExtension.jl index fa65c6e481..84770ff751 100644 --- a/ext/HomotopyContinuationExtension.jl +++ b/ext/HomotopyContinuationExtension.jl @@ -7,6 +7,6 @@ import HomotopyContinuation as HC import Symbolics: unwrap, wrap, Rewriters, symtype, issym, istree # Creates and exports hc_steady_states function. -include("HomotopyContinuationExtension/homotopy_continuation_extension.jl") +include("CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl") end \ No newline at end of file From a6ddb9693e77c8c0be0748506bc46ca85280e2e9 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 25 Sep 2023 17:11:13 -0400 Subject: [PATCH 09/30] complete extension name update --- HISTORY.md | 2 +- ...ionExtension.jl => CatalystHomotopyContinuationExtension.jl} | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) rename ext/{HomotopyContinuationExtension.jl => CatalystHomotopyContinuationExtension.jl} (86%) diff --git a/HISTORY.md b/HISTORY.md index 521518b82d..34440ef5c8 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -1,7 +1,7 @@ # Breaking updates and feature summaries across releases ## Catalyst unreleased (master branch) -- Added a HomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds teh steady states of a reactin system using the homotopy continuation method. +- Added a CatalystHomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds teh steady states of a reactin system using the homotopy continuation method. ## Catalyst 13.4 - Added the ability to create species that represent chemical compounds and know diff --git a/ext/HomotopyContinuationExtension.jl b/ext/CatalystHomotopyContinuationExtension.jl similarity index 86% rename from ext/HomotopyContinuationExtension.jl rename to ext/CatalystHomotopyContinuationExtension.jl index 84770ff751..4a427c8bd1 100644 --- a/ext/HomotopyContinuationExtension.jl +++ b/ext/CatalystHomotopyContinuationExtension.jl @@ -1,4 +1,4 @@ -module HomotopyContinuationExtension +module CatalystHomotopyContinuationExtension # Fetch packages. using Catalyst From 95bc90d608d3dd1607cbdda7889349c5b1fc4de2 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 25 Sep 2023 19:55:57 -0400 Subject: [PATCH 10/30] 1.6 compat --- src/Catalyst.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index b33b6782d0..f660569954 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -106,5 +106,8 @@ export balance_reaction # HomotopyContinuation function hc_steady_states end export hc_steady_states +if !isdefined(Base, :get_extension) # Workaroudn for 1.6 compatibility + include("../ext/CatalystHomotopyContinuationExtension.jl") + end end # module From 0803dc25e7892d10e74f5c7285e4569f7fb2599c Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 25 Sep 2023 19:56:23 -0400 Subject: [PATCH 11/30] 1.6 compat. --- src/Catalyst.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index f660569954..1806ad2175 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -106,8 +106,8 @@ export balance_reaction # HomotopyContinuation function hc_steady_states end export hc_steady_states -if !isdefined(Base, :get_extension) # Workaroudn for 1.6 compatibility +if !isdefined(Base, :get_extension) # Workaround for 1.6 compatibility. include("../ext/CatalystHomotopyContinuationExtension.jl") - end +end end # module From a7c18c51f5000c87f20d754753d21f6eae326c0e Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 25 Sep 2023 21:01:39 -0400 Subject: [PATCH 12/30] Only test extensions when VERSION >= v"1.9" --- test/runtests.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index cc7a1135c1..7d659feb5c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -48,5 +48,7 @@ using SafeTestsets end ### Tests extensions. ### - @time @safetestset "Homotopy Continuation Extension" begin include("extensions/homotopy_continuation.jl") end + if VERSION >= v"1.9" + @time @safetestset "Homotopy Continuation Extension" begin include("extensions/homotopy_continuation.jl") end + end end # @time From d30be1d9e4ee658525eda4bc007c18809fc0111d Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 29 Sep 2023 17:47:29 -0400 Subject: [PATCH 13/30] small doc improvemtns --- docs/src/catalyst_applications/homotopy_continuation.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/catalyst_applications/homotopy_continuation.md b/docs/src/catalyst_applications/homotopy_continuation.md index c02d76ebcd..7754301c7d 100644 --- a/docs/src/catalyst_applications/homotopy_continuation.md +++ b/docs/src/catalyst_applications/homotopy_continuation.md @@ -9,7 +9,7 @@ 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. -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. +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 publication](@ref homotopy_continuation_citation). ## Basic example For this tutorial, we will use a model from Wilhem (2009)[^1] (which @@ -62,7 +62,7 @@ hc_steady_states(wilhelm_2009_model, ps; u0=u0) - 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. --- -## Citations +## [Citation](@id homotopy_continuation_citation) If you use this functionality in your reasearch, please cite the following paper to support the authors of the HomotopyContinuation package: ``` @inproceedings{HomotopyContinuation.jl, From 2614d099a84566c16898c2d38b3133ac683cffae Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 2 Oct 2023 10:33:13 -0400 Subject: [PATCH 14/30] wip --- .../homotopy_continuation_extension.jl | 2 +- src/miscellaneous_functions.jl | 38 +++++++++++++++++++ 2 files changed, 39 insertions(+), 1 deletion(-) create mode 100644 src/miscellaneous_functions.jl diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 53eb0258bb..5f5691773b 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -13,7 +13,7 @@ Arguments: - `u0=typeof(ps)()`: Initial conditions for which we want to find the steady states. For systems with conservation laws this are required to compute conserved quantities. - `kwargs...`: any additional arguments (like `show_progress= true`) are passed into HomotopyContinuation.jl's `solve` call. -Examples: +Examples ```@repl rs = @reaction_network begin k1, Y --> 2X diff --git a/src/miscellaneous_functions.jl b/src/miscellaneous_functions.jl new file mode 100644 index 0000000000..3675fe5fad --- /dev/null +++ b/src/miscellaneous_functions.jl @@ -0,0 +1,38 @@ +# Contains various functions that does not clearly belong to any other files. Hopefully these can eventually be moved to a better place. + +### Steady State Stability ### + +""" + steady_state_stability(s_state, ps, sys) + +Computes the stability of steady states. For `Vector{Float64}` input, returns a `Bool` corresponding the state's stability. For `Vector{Vector{Float64}}` input, returns a `Vector{Bool}` corresponding the states' stability. + +Arguments: +- `s_state`: The steate(s) for which we want to compute the stability. These are assumed to be system steady states. +- `ps` the parameter values for which we wish to compute stability. +- `sys`: The system for which we wish to comptue stability. Can be a `ReactionSystem` or an `ODEFunction`. + +Example +```@repl +rs = @reaction_network begin + (p,d), 0 <--> X +end +ps = [:p => 1.0, :d => 0.2] +steady_state_stability([5.0], ps, rs) +``` +gives +``` +true +``` + +Notes: +- If a `ReactionSystem` is given, it is coverted to a `ODEFunction` internally. If you plan to compute the stability for a large number of steady states. First covert your ReactionSystem to a ODEFunction using `ODEFunction(convert(ODESystem, rs); jac=true)` and then use this in your input. + + +""" +steady_state_stability(s_states::Vector{Float64}, ps, rs::ReactionSystem) = steady_state_stability(s_states, ODEFunction(convert(ODESystem, rs); jac=true)) +# If a ODESystem is given directyl, no conversion is needed. +function steady_state_stability(s_states::Vector{Float64}, ps, ofun::ODEFunction) + maximum(eigen(ofun.jac(ss, last.(ps), Inf)).values) < 0.0 +end +steady_state_stability(s_states::Vector{Vector{Float64}}, args...) = steady_state_stability.(s_states, args...) From 2ba7b1628a317ca75d87dffb61a087bb63878e04 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:03:12 -0400 Subject: [PATCH 15/30] Update docs/src/catalyst_applications/homotopy_continuation.md Co-authored-by: Sam Isaacson --- docs/src/catalyst_applications/homotopy_continuation.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/catalyst_applications/homotopy_continuation.md b/docs/src/catalyst_applications/homotopy_continuation.md index 7754301c7d..7477a4e772 100644 --- a/docs/src/catalyst_applications/homotopy_continuation.md +++ b/docs/src/catalyst_applications/homotopy_continuation.md @@ -50,7 +50,7 @@ two_state_model = @reaction_network begin (k1,k2), X1 <--> X2 end ``` -To find the steady states of these, an initial condition must also be provided (which is used to compute the system's conserved quantities, in this case `X1+X2`): +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 ps = [:k1 => 2.0, :k2 => 1.0] u0 = [:X1 => 1.0, :X2 => 1.0] From 2795ae2b02211ba03c263abdb67e5065571ce5b5 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:03:24 -0400 Subject: [PATCH 16/30] Update docs/src/catalyst_applications/homotopy_continuation.md Co-authored-by: Sam Isaacson --- docs/src/catalyst_applications/homotopy_continuation.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/src/catalyst_applications/homotopy_continuation.md b/docs/src/catalyst_applications/homotopy_continuation.md index 7477a4e772..c53d479d05 100644 --- a/docs/src/catalyst_applications/homotopy_continuation.md +++ b/docs/src/catalyst_applications/homotopy_continuation.md @@ -54,7 +54,8 @@ Catalyst allows the conservation laws of such systems to be computed using the ` ```@example hc3 ps = [:k1 => 2.0, :k2 => 1.0] u0 = [:X1 => 1.0, :X2 => 1.0] -hc_steady_states(wilhelm_2009_model, ps; u0=u0) +hc_steady_states(wilhelm_2009_model, ps; u0) + ``` ## Final notes From d4f154cc13e3a22b50b321747e290cd72fa3a429 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:03:37 -0400 Subject: [PATCH 17/30] Update docs/src/catalyst_applications/homotopy_continuation.md Co-authored-by: Sam Isaacson --- docs/src/catalyst_applications/homotopy_continuation.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/src/catalyst_applications/homotopy_continuation.md b/docs/src/catalyst_applications/homotopy_continuation.md index c53d479d05..b104108e18 100644 --- a/docs/src/catalyst_applications/homotopy_continuation.md +++ b/docs/src/catalyst_applications/homotopy_continuation.md @@ -59,7 +59,8 @@ hc_steady_states(wilhelm_2009_model, ps; u0) ``` ## Final notes -- `hc_steady_states` supports any systems where all rates are systems of rational polynomial (such as hill functions with integer hill coefficients). +- `hc_steady_states` supports any systems where all rates are systems of rational polynomials (such as Hill functions with integer Hill coefficients). + - 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. --- From 13447e12e167ba27461b424ebdb8ff75d247ce2e Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:03:57 -0400 Subject: [PATCH 18/30] Update ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl Co-authored-by: Sam Isaacson --- .../homotopy_continuation_extension.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 53eb0258bb..5dbc15f98e 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -3,7 +3,7 @@ """ hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=typeof(ps)(), kwargs...) -Uses homotopy continuation to find the steady states of the ODE system corresponding to the provided reaction system. +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. From ad1a7c4898971a2aad80e00a48e1de75a16a8a60 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:04:15 -0400 Subject: [PATCH 19/30] Update ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl Co-authored-by: Sam Isaacson --- .../homotopy_continuation_extension.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 5dbc15f98e..c70449f91b 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -32,7 +32,6 @@ gives ``` Notes: -- This is a wrapper around the `solve` function provided by HomotopyContinuation.jl, all credit for this functionality to that package's authors. ``` """ function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=typeof(ps)(), kwargs...) From f70bca6327d77b47308566c4fd1cbc1de32b44cc Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:05:35 -0400 Subject: [PATCH 20/30] Update ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl Co-authored-by: Sam Isaacson --- .../homotopy_continuation_extension.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index c70449f91b..91211dc997 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -38,7 +38,7 @@ function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, 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=neg_thres) : sols) + return (filter_negative ? filter_negative_f(sols; neg_thres) : sols) end # For a given reaction system, paraemter values, and initial conditions, find the polynomial that HC solves to find steady states. From 242380dd72560d21530d00f064a9c7b66bb6d371 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:06:10 -0400 Subject: [PATCH 21/30] Update ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl Co-authored-by: Sam Isaacson --- .../homotopy_continuation_extension.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 91211dc997..283843a582 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -57,7 +57,8 @@ end # If u0s are not given while conservation laws are present, throws an error. function conservationlaw_errorcheck(rs, pre_varmap) - vars_with_vals = union(first.(pre_varmap), keys(ModelingToolkit.defaults(rs))) + vars_with_vals = Set(p[1] for p in pre_varpmap) + union!(vars_with_vals, keys(ModelingToolkit.defaults(rs)) isempty(intersect(species(rs), vars_with_vals)) || return isempty(conservedequations(rs)) && return error("The system have conservation laws but no initial conditions were provided. Please provide initial conditions.") From dd443b1fe43b74624550cb48ba1c8aaa55c63750 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:06:23 -0400 Subject: [PATCH 22/30] Update ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl Co-authored-by: Sam Isaacson --- .../homotopy_continuation_extension.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 283843a582..173483d7d5 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -59,7 +59,7 @@ end function conservationlaw_errorcheck(rs, pre_varmap) vars_with_vals = Set(p[1] for p in pre_varpmap) union!(vars_with_vals, keys(ModelingToolkit.defaults(rs)) - isempty(intersect(species(rs), vars_with_vals)) || return + issubset(species(rs), vars_with_vals) && return isempty(conservedequations(rs)) && return error("The system have conservation laws but no initial conditions were provided. Please provide initial conditions.") end From dbe2f75a5c01562bc42a794bc452aa974f1ef4fc Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:06:38 -0400 Subject: [PATCH 23/30] Update ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl Co-authored-by: Sam Isaacson --- .../homotopy_continuation_extension.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 173483d7d5..a478e3d074 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -60,8 +60,8 @@ function conservationlaw_errorcheck(rs, pre_varmap) vars_with_vals = Set(p[1] for p in pre_varpmap) union!(vars_with_vals, keys(ModelingToolkit.defaults(rs)) issubset(species(rs), vars_with_vals) && return - isempty(conservedequations(rs)) && return - error("The system have conservation laws but no initial conditions were provided. Please provide initial conditions.") + isempty(conservedequations(rs)) || + error("The system has conservation laws but initial conditions were not provided for all species.") end # Unfolds a function (like mm or hill). From 996ade5b77289ac19a9dad96533fdf10d8c463ad Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:07:00 -0400 Subject: [PATCH 24/30] Update ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl Co-authored-by: Sam Isaacson --- .../homotopy_continuation_extension.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index a478e3d074..87b8cda7d1 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -118,7 +118,7 @@ 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.0) && (sol[idx] = 0.0) + (neg_thres < sol[idx] < 0) && (sol[idx] = 0) end return filter(sol -> all(>=(0.0), sol), sols) end \ No newline at end of file From aa0ab4c51420f56a40b7204fbe555ebed7dbd1f4 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:07:08 -0400 Subject: [PATCH 25/30] Update ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl Co-authored-by: Sam Isaacson --- .../homotopy_continuation_extension.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 87b8cda7d1..6124722b46 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -120,5 +120,5 @@ 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.0), sol), sols) + return filter(sol -> all(>=(0), sol), sols) end \ No newline at end of file From fe979ed2e8315059dca162b31f3868b46385f600 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:08:52 -0400 Subject: [PATCH 26/30] updates --- HISTORY.md | 12 +++++++++++- .../catalyst_applications/homotopy_continuation.md | 2 +- docs/src/index.md | 2 +- src/Catalyst.jl | 3 --- 4 files changed, 13 insertions(+), 6 deletions(-) diff --git a/HISTORY.md b/HISTORY.md index 34440ef5c8..785d985f73 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -1,7 +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 teh steady states of a reactin system using the homotopy continuation method. +- 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 diff --git a/docs/src/catalyst_applications/homotopy_continuation.md b/docs/src/catalyst_applications/homotopy_continuation.md index 7754301c7d..ab8a90307b 100644 --- a/docs/src/catalyst_applications/homotopy_continuation.md +++ b/docs/src/catalyst_applications/homotopy_continuation.md @@ -9,7 +9,7 @@ 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. -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 publication](@ref homotopy_continuation_citation). +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 diff --git a/docs/src/index.md b/docs/src/index.md index d182e7a813..f02599146b 100644 --- a/docs/src/index.md +++ b/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: diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 1806ad2175..b33b6782d0 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -106,8 +106,5 @@ export balance_reaction # HomotopyContinuation function hc_steady_states end export hc_steady_states -if !isdefined(Base, :get_extension) # Workaround for 1.6 compatibility. - include("../ext/CatalystHomotopyContinuationExtension.jl") -end end # module From 783cd9abd4fb26aa8c646df3022b98a5107878d5 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:22:45 -0400 Subject: [PATCH 27/30] updates --- .../homotopy_continuation_extension.jl | 17 ++++++++++------- test/extensions/homotopy_continuation.jl | 3 +++ 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 9e873f2b85..234fa4384f 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -8,9 +8,9 @@ Uses homotopy continuation via HomotopyContinuation.jl to find the steady states 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 <0 is removed from the output. +- `filter_negative=true`: If set to true, solutions with any species concentration neg_thres`` but `< 0.0` are set to `0.0`. -- `u0=typeof(ps)()`: Initial conditions for which we want to find the steady states. For systems with conservation laws this are required to compute conserved quantities. +- `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. - `kwargs...`: any additional arguments (like `show_progress= true`) are passed into HomotopyContinuation.jl's `solve` call. Examples @@ -34,7 +34,7 @@ gives Notes: ``` """ -function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=typeof(ps)(), kwargs...) +function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=nothing, kwargs...) ss_poly = steady_state_polynomial(rs, ps, u0) sols = HC.real_solutions(HC.solve(ss_poly; kwargs...)) reorder_sols!(sols, ss_poly, rs) @@ -44,7 +44,7 @@ end # For a given reaction system, paraemter 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 = map(pair -> pair::Pair{Num,Float64}, [symmap_to_varmap(rs,u0); symmap_to_varmap(rs,ps)]) + pre_varmap = map(pair -> pair, [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) @@ -58,7 +58,6 @@ 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_varpmap) - union!(vars_with_vals, keys(ModelingToolkit.defaults(rs)) issubset(species(rs), vars_with_vals) && return isempty(conservedequations(rs)) || error("The system has conservation laws but initial conditions were not provided for all species.") @@ -89,8 +88,12 @@ end 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) == ^) && isinteger(arguments(expr)[2]) - return arguments(expr)[1] ^ Int64(arguments(expr)[2]) + 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 diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index cbcfb6b109..8174f6ec2e 100644 --- a/test/extensions/homotopy_continuation.jl +++ b/test/extensions/homotopy_continuation.jl @@ -79,6 +79,7 @@ end # Tests that non-scalar reaction rates work. # Tests that rational polynomial steady state systems work. # Test filter_negative=false works. +# tests than non-integer exponents throws an error. let rs = @reaction_network begin 0.1 + hill(X,v,K,n), 0 --> X @@ -92,4 +93,6 @@ let for ss in sss @test f(sss[1], last.(ps), 0.0)[1] == 0.0 end + + @test_throws Exception hc_steady_states(rs, [:v => 5.0, :K => 2.5, :n => 2.7, :d => 1.0]; show_progress=false) end \ No newline at end of file From 2804db160baff8ee7566fd655a787b46bd799415 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 3 Oct 2023 19:48:12 -0400 Subject: [PATCH 28/30] updates --- .../homotopy_continuation_extension.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 234fa4384f..40f74debb5 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -34,7 +34,7 @@ gives Notes: ``` """ -function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=nothing, kwargs...) +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) @@ -57,8 +57,8 @@ 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_varpmap) - issubset(species(rs), vars_with_vals) && return + vars_with_vals = Set(p[1] for p in pre_varmap) + isempty(intersect(species(rs), vars_with_vals)) || return isempty(conservedequations(rs)) || error("The system has conservation laws but initial conditions were not provided for all species.") end From 168044059f44196ac28557a12edd37c1ce7a292a Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 4 Oct 2023 09:37:39 -0400 Subject: [PATCH 29/30] update --- .../homotopy_continuation.md | 4 +- .../homotopy_continuation_extension.jl | 16 ++++---- src/miscellaneous_functions.jl | 38 ------------------- test/extensions/homotopy_continuation.jl | 11 +++--- 4 files changed, 16 insertions(+), 53 deletions(-) delete mode 100644 src/miscellaneous_functions.jl diff --git a/docs/src/catalyst_applications/homotopy_continuation.md b/docs/src/catalyst_applications/homotopy_continuation.md index 0ef5c44b36..bda625a9d8 100644 --- a/docs/src/catalyst_applications/homotopy_continuation.md +++ b/docs/src/catalyst_applications/homotopy_continuation.md @@ -60,12 +60,12 @@ hc_steady_states(wilhelm_2009_model, ps; u0) ## 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. --- ## [Citation](@id homotopy_continuation_citation) -If you use this functionality in your reasearch, please cite the following paper to support the authors of the HomotopyContinuation package: +If you use this functionality in your research, please cite the following paper to support the authors of the HomotopyContinuation package: ``` @inproceedings{HomotopyContinuation.jl, title={{H}omotopy{C}ontinuation.jl: {A} {P}ackage for {H}omotopy {C}ontinuation in {J}ulia}, diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 40f74debb5..be29fa28d1 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -9,8 +9,8 @@ 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`` 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. +- `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 @@ -41,10 +41,10 @@ function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, return (filter_negative ? filter_negative_f(sols; neg_thres) : sols) end -# For a given reaction system, paraemter values, and initial conditions, find the polynomial that HC solves to find steady states. +# 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 = map(pair -> pair, [symmap_to_varmap(rs,u0); symmap_to_varmap(rs,ps)]) + pre_varmap = map(pair -> pair, [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) @@ -58,9 +58,9 @@ 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) - isempty(intersect(species(rs), vars_with_vals)) || return + 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 all species.") + error("The system has conservation laws but initial conditions were not provided for some species.") end # Unfolds a function (like mm or hill). @@ -105,12 +105,12 @@ function remove_denominators(expr) return remove_denominators(arguments(s_expr)[1]) end if operation(s_expr) == + - return +(remove_denominators.(arguments(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 lexiographic values of the variable names. This reorders the output acording to the species index in the reaction system species vector. +# 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] diff --git a/src/miscellaneous_functions.jl b/src/miscellaneous_functions.jl deleted file mode 100644 index 3675fe5fad..0000000000 --- a/src/miscellaneous_functions.jl +++ /dev/null @@ -1,38 +0,0 @@ -# Contains various functions that does not clearly belong to any other files. Hopefully these can eventually be moved to a better place. - -### Steady State Stability ### - -""" - steady_state_stability(s_state, ps, sys) - -Computes the stability of steady states. For `Vector{Float64}` input, returns a `Bool` corresponding the state's stability. For `Vector{Vector{Float64}}` input, returns a `Vector{Bool}` corresponding the states' stability. - -Arguments: -- `s_state`: The steate(s) for which we want to compute the stability. These are assumed to be system steady states. -- `ps` the parameter values for which we wish to compute stability. -- `sys`: The system for which we wish to comptue stability. Can be a `ReactionSystem` or an `ODEFunction`. - -Example -```@repl -rs = @reaction_network begin - (p,d), 0 <--> X -end -ps = [:p => 1.0, :d => 0.2] -steady_state_stability([5.0], ps, rs) -``` -gives -``` -true -``` - -Notes: -- If a `ReactionSystem` is given, it is coverted to a `ODEFunction` internally. If you plan to compute the stability for a large number of steady states. First covert your ReactionSystem to a ODEFunction using `ODEFunction(convert(ODESystem, rs); jac=true)` and then use this in your input. - - -""" -steady_state_stability(s_states::Vector{Float64}, ps, rs::ReactionSystem) = steady_state_stability(s_states, ODEFunction(convert(ODESystem, rs); jac=true)) -# If a ODESystem is given directyl, no conversion is needed. -function steady_state_stability(s_states::Vector{Float64}, ps, ofun::ODEFunction) - maximum(eigen(ofun.jac(ss, last.(ps), Inf)).values) < 0.0 -end -steady_state_stability(s_states::Vector{Vector{Float64}}, args...) = steady_state_stability.(s_states, args...) diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index 8174f6ec2e..144774db80 100644 --- a/test/extensions/homotopy_continuation.jl +++ b/test/extensions/homotopy_continuation.jl @@ -6,7 +6,7 @@ import HomotopyContinuation # Tests for network without conservation laws. # Tests for Symbol parameter input. -# Tests for Symbolics initial condiiton input. +# Tests for Symbolics initial condition input. # Tests for different types (Symbol/Symbolics) for parameters and initial conditions. # Tests that attempts to find steady states of system with conservation laws, while u0 is not provided, gives an error. let @@ -47,8 +47,9 @@ let end # Tests that reordering is correct. -# Tests corectness in presence of default values. -# Tests where some defaul values are overwritten with other values +# Tests correctness in presence of default values. +# Tests where some default values are overwritten with other values. +# Tests where input ps/u0 are tuples with mixed types. let rs_1 = @reaction_network begin @parameters kX1=1.0 kX2=2.0 kY1=12345.0 @@ -57,8 +58,8 @@ let (kY1,kY2), Y1 <--> Y2 (kZ1,kZ2), Z1 <--> Z2 end - ps = [:kY1 => 1.0, :kY2 => 3.0, :kZ1 => 1.0, :kZ2 => 4.0] - u0_1 = [:Y1 => 1.0, :Y2 => 3.0, :Z1 => 10.0, :Z2 =>40.0] + ps = (:kY1 => 1.0, :kY2 => 3, :kZ1 => 1.0, :kZ2 => 4.0) + u0_1 = (:Y1 => 1.0, :Y2 => 3, :Z1 => 10, :Z2 =>40.0) ss_1 = sort(hc_steady_states(rs_1, ps; u0=u0_1, show_progress=false), by=sol->sol[1]) @test ss_1 ≈ [[0.2, 0.1, 3.0, 1.0, 40.0, 10.0]] From a03f3a5ff4de0eb2cfb74c5fc66e180695deb145 Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 4 Oct 2023 09:40:56 -0400 Subject: [PATCH 30/30] update --- .../homotopy_continuation_extension.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index be29fa28d1..d3f5343cc1 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -44,7 +44,7 @@ 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 = map(pair -> pair, [symmap_to_varmap(rs,u0)..., symmap_to_varmap(rs,ps)...]) + 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)