diff --git a/HISTORY.md b/HISTORY.md index b3bc6d0366..785d985f73 100644 --- a/HISTORY.md +++ b/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 diff --git a/Project.toml b/Project.toml index e1198f1ef7..66e1708d97 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] +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"] diff --git a/docs/src/catalyst_applications/homotopy_continuation.md b/docs/src/catalyst_applications/homotopy_continuation.md index 4d04108de3..bda625a9d8 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. 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 @@ -28,116 +25,60 @@ 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. -```@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) \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index add160c4d5..223e5fa0b8 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/ext/CatalystHomotopyContinuationExtension.jl b/ext/CatalystHomotopyContinuationExtension.jl new file mode 100644 index 0000000000..4a427c8bd1 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl new file mode 100644 index 0000000000..d3f5343cc1 --- /dev/null +++ b/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`` 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 \ 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/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 new file mode 100644 index 0000000000..144774db80 --- /dev/null +++ b/test/extensions/homotopy_continuation.jl @@ -0,0 +1,99 @@ +### Fetch Packages ### +using Catalyst, Test +import HomotopyContinuation + +### Run Tests ### + +# Tests for network without conservation laws. +# Tests for Symbol parameter 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 + 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] + + 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. +# 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, 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, 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 + +# Tests that reordering is correct. +# 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 + @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, :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]] + + 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, 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. +# tests than non-integer exponents throws an error. +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 + + @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 diff --git a/test/runtests.jl b/test/runtests.jl index 5a772722f2..7d659feb5c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -46,4 +46,9 @@ using SafeTestsets if !Sys.isapple() @time @safetestset "Graphs" begin include("visualization/graphs.jl") end end + + ### Tests extensions. ### + if VERSION >= v"1.9" + @time @safetestset "Homotopy Continuation Extension" begin include("extensions/homotopy_continuation.jl") end + end end # @time