You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardExpand all lines: HISTORY.md
+11Lines changed: 11 additions & 0 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -1,6 +1,17 @@
1
1
# Breaking updates and feature summaries across releases
2
2
3
3
## Catalyst unreleased (master branch)
4
+
- 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:
5
+
```julia
6
+
wilhelm_2009_model =@reaction_networkbegin
7
+
k1, Y -->2X
8
+
k2, 2X --> X + Y
9
+
k3, X + Y --> Y
10
+
k4, X -->0
11
+
end
12
+
ps = [:k1=>8.0, :k2=>2.0, :k3=>1.0, :k4=>1.5]
13
+
hc_steady_states(wilhelm_2009_model, ps)
14
+
```
4
15
5
16
## Catalyst 13.4
6
17
- Added the ability to create species that represent chemical compounds and know
In this tutorial, we will demonstrate how homotopy continuation can be used to
12
-
find the steady states of mass action chemical reaction networks implemented in
13
-
Catalyst.
11
+
12
+
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.
14
13
15
14
## Basic example
16
15
For this tutorial, we will use a model from Wilhem (2009)[^1] (which
17
16
demonstrates bistability in a small chemical reaction network). We declare the
18
17
model and the parameter set for which we want to find the steady states:
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.
Finally, we use Catalyst's `to_multivariate_poly` function to reinterpret our
48
-
symbolic equations in a polynomial representation that is compatible with
49
-
HomotopyContinuation. We can then apply HomotopyContinuation's `solve` command
50
-
to find the roots, using `real_solutions` to filter our non-physical complex
51
-
steady-states:
33
+
Now we can find the steady states using:
52
34
```@example hc1
53
-
polyeqs = Catalyst.to_multivariate_poly(neweqs)
54
-
sols = HC.real_solutions(HC.solve(polyeqs))
35
+
hc_steady_states(wilhelm_2009_model, ps)
55
36
```
56
-
Note that HomotopyContinuation orders variables lexicographically, so this will
57
-
be the ordering present in each steady-state solution vector (i.e. `[X1, X2]` is
58
-
the ordering here).
37
+
The order of the species in the output vectors are the same as in `species(wilhelm_2009_model)`.
38
+
39
+
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.
59
40
60
-
While it is not the case for this CRN, we note that solutions with negative
61
-
species concentrations can be valid (unphysical) steady-states for certain
62
-
systems. These will need to be filtered out as well.
63
41
64
42
## Systems with conservation laws
65
-
Finally, some systems are under-determined, and have an infinite number of
66
-
possible steady states. These are typically systems containing a conservation
43
+
Some systems are under-determined, and have an infinite number of possible steady states. These are typically systems containing a conservation
67
44
law, e.g.
68
45
```@example hc3
69
46
using Catalyst
70
47
import HomotopyContinuation
71
-
const MT = ModelingToolkit
72
-
const HC = HomotopyContinuation
73
48
74
49
two_state_model = @reaction_network begin
75
50
(k1,k2), X1 <--> X2
76
51
end
77
52
```
78
-
Catalyst allows the conservation laws to be computed using the
79
-
`conservationlaws` function. By using these to reduce the dimensionality of the
80
-
system, as well specifying the initial amount of each species,
81
-
HomotopyContinuation can again be used to find steady states. First, we set the
82
-
default values of the system's initial conditions and parameter values. This
83
-
will allow the system to automatically find the conserved amounts.
Again, we next create a dictionary for parameter values that we substitute in to
93
-
give our final equation.
53
+
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`):
-`hc_steady_states` supports any systems where all rates are systems of rational polynomials (such as Hill functions with integer Hill coefficients).
63
+
- 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.
64
+
- 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.
65
+
---
117
66
118
-
# for each SS, set X1's value in the subs map and calculate X2
119
-
@unpack X1 = two_state_model
120
-
X2 = map(sols) do x
121
-
X1val = x[1]
122
-
subs[MT.value(X1)] = X1val
123
-
substitute(X2eq, subs)
124
-
end
67
+
## [Citation](@id homotopy_continuation_citation)
68
+
If you use this functionality in your research, please cite the following paper to support the authors of the HomotopyContinuation package:
125
69
```
126
-
giving that the steady-state for `X2` is about `1.33333`.
127
-
128
-
As an alternative, we could have coupled `neweqs` with the conservation law
129
-
relations to have HomotopyContinuation find the steady-states simultaneously:
130
-
```@example hc3
131
-
# move all the terms in the conserved equations to one side
title={{H}omotopy{C}ontinuation.jl: {A} {P}ackage for {H}omotopy {C}ontinuation in {J}ulia},
72
+
author={Breiding, Paul and Timme, Sascha},
73
+
booktitle={International Congress on Mathematical Software},
74
+
pages={458--465},
75
+
year={2018},
76
+
organization={Springer}
77
+
}
140
78
```
141
-
---
79
+
142
80
## References
143
-
[^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)
81
+
[^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)
82
+
[^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)
83
+
[^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)
84
+
[^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)
Uses homotopy continuation via HomotopyContinuation.jl to find the steady states of the ODE system corresponding to the provided reaction system.
7
+
8
+
Arguments:
9
+
- `rs::ReactionSystem`: The reaction system for which we want to find the steady states.
10
+
- `ps`: The parameter values for which we want to find the steady states.
11
+
- `filter_negative=true`: If set to true, solutions with any species concentration <neg_thres is removed from the output.
12
+
- `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`.
13
+
- `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).
14
+
- `kwargs...`: any additional arguments (like `show_progress= true`) are passed into HomotopyContinuation.jl's `solve` call.
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.")
96
+
end
97
+
end
98
+
end
99
+
100
+
# If the input is a fraction, removes the denominator.
101
+
functionremove_denominators(expr)
102
+
s_expr =simplify_fractions(expr)
103
+
!istree(expr) &&return expr
104
+
ifoperation(s_expr) ==/
105
+
returnremove_denominators(arguments(s_expr)[1])
106
+
end
107
+
ifoperation(s_expr) ==+
108
+
returnsum(remove_denominators(arg) for arg inarguments(s_expr))
109
+
end
110
+
return s_expr
111
+
end
112
+
113
+
# 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.
0 commit comments