Merge pull request #739 from SciML/new_stability_computation

[v14 ready] New stability computation

Author: TorkelE

HISTORY.md

@@ -86,6 +86,20 @@ plot(bif_dia; xguide="k1", yguide="X")
 ```
 - Automatically handles elimination of conservation laws for computing bifurcation diagrams.
 - Updated Bifurcation documentation with respect to this new feature.
+- Added function `isautonomous` to check if a `ReactionSystem` is autonomous.
+- Added function `steady_state_stability` to compute stability for steady states. Example:
+```julia
+# Creates model.
+rn = @reaction_network begin
+    (p,d), 0 <--> X
+end
+p = [:p => 1.0, :d => 0.5]
+
+# Finds (the trivial) steady state, and computes stability.
+steady_state = [2.0]
+steady_state_stability(steady_state, rn, p)
+```
+Here, `steady_state_stability` take an optional argument `tol = 10*sqrt(eps())`, which is used to determine whether a eigenvalue real part is reliably less that 0.
 
 ## Catalyst 13.5
 - Added a CatalystHomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds the steady states of a reaction system using the homotopy continuation method. This feature is only available for julia versions 1.9+. Example:

docs/pages.jl

@@ -37,7 +37,7 @@ pages = Any[
     "Steady state analysis" => Any[
         "steady_state_functionality/homotopy_continuation.md",
         "steady_state_functionality/nonlinear_solve.md",
-        # Stability analysis.
+        "steady_state_functionality/steady_state_stability_computation.md",        
         "steady_state_functionality/bifurcation_diagrams.md"
         # Dynamical systems analysis.
     ],

docs/src/api.md

@@ -163,6 +163,7 @@ speciesmap
 paramsmap
 reactionsystemparamsmap
 isspecies
+isautonomous
 Catalyst.isconstant
 Catalyst.isbc
 ```

docs/src/steady_state_functionality/steady_state_stability_computation.md

@@ -0,0 +1,63 @@
+# Steady state stability computation
+After system steady states have been found using [HomotopyContinuation.jl](@ref homotopy_continuation), [NonlinearSolve.jl](@ref nonlinear_solve), or other means, their stability can be computed using Catalyst's `steady_state_stability` function. Systems with conservation laws will automatically have these removed, permitting stability computation on systems with singular Jacobian.
+
+!!! warn 
+    Catalyst currently computes steady state stabilities using the naive approach of checking whether a system's largest eigenvalue real part is negative. While more advanced stability computation methods exist (and would be a welcome addition to Catalyst), there is no direct plans to implement these. Furthermore, Catalyst uses a tolerance `tol = 10*sqrt(eps())` to determine whether a computed eigenvalue is far away enough from 0 to be reliably used. This threshold can be changed through the `tol` keyword argument.
+
+## Basic examples
+Let us consider the following basic example:
+```@example stability_1
+using Catalyst
+rn = @reaction_network begin 
+    (p,d), 0 <--> X
+end
+```
+It has a single (stable) steady state at $X = p/d$. We can confirm stability using the `steady_state_stability` function, to which we provide the steady state, the reaction system, and the parameter values:
+```@example stability_1
+ps = [:p => 2.0, :d => 0.5]
+steady_state = [:X => 4.0]
+steady_state_stability(steady_state, rn, ps)
+```
+
+Next, let us consider the following [self-activation loop](@ref ref):
+```@example stability_1
+sa_loop = @reaction_network begin 
+    (hill(X,v,K,n),d), 0 <--> X
+end
+```
+For certain parameter choices, this system exhibits multi-stability. Here, we can find the steady states [using homotopy continuation](@ref homotopy_continuation):
+```@example stability_1
+import HomotopyContinuation
+ps = [:v => 2.0, :K => 0.5, :n => 3, :d => 1.0]
+steady_states = hc_steady_states(sa_loop, ps)
+```
+Next, we can apply `steady_state_stability` to each steady state yielding a vector of their stabilities:
+```@example stability_1
+[steady_state_stability(sstate, sa_loop, ps) for sstate in steady_states]
+```
+
+Finally, as described above, Catalyst uses an optional argument, `tol`, to determine how strict to make the stability check.  I.e. below we set the tolerance to `1e-6` (a larger value, that is stricter, than the default of `10*sqrt(eps())`)
+```@example stability_1
+[steady_state_stability(sstate, sa_loop, ps; tol = 1e-6) for sstate in steady_states]
+nothing# hide
+```
+
+## Pre-computing the Jacobian to increase performance when computing stability for many steady states
+Catalyst uses the system Jacobian to compute steady state stability, and the Jacobian is computed once for each call to `steady_state_stability`. If you repeatedly compute stability for steady states of the same system, pre-computing the Jacobian and supplying it to the `steady_state_stability` function can improve performance. 
+
+In this example we use the self-activation loop from previously, pre-computes its Jacobian, and uses it to multiple `steady_state_stability` calls:
+```@example stability_1
+ss_jac = steady_state_jac(sa_loop)
+
+ps_1 = [:v => 2.0, :K => 0.5, :n => 3, :d => 1.0]
+steady_states_1 = hc_steady_states(sa_loop, ps)
+stability_1 = steady_state_stability(steady_states_1, sa_loop, ps_1; ss_jac=ss_jac)
+
+ps_2 = [:v => 4.0, :K => 1.5, :n => 2, :d => 1.0]
+steady_states_2 = hc_steady_states(sa_loop, ps)
+stability_2 = steady_state_stability(steady_states_2, sa_loop, ps_2; ss_jac=ss_jac)
+nothing # hide
+```
+
+!!! warn
+    For systems with [conservation laws](@ref homotopy_continuation_conservation_laws), `steady_state_jac` must be supplied a `u0` vector (indicating species concentrations for conservation law computation). This is required to eliminate the conserved quantities, preventing a singular Jacobian. These are supplied using the `u0` optional argument.

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

@@ -3,6 +3,9 @@
 # Creates a BifurcationProblem, using a ReactionSystem as an input.
 function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...;
                                plot_var=nothing, record_from_solution=BK.record_sol_default, jac=true, u0=[], kwargs...)
+    if !isautonomous(rs) 
+        error("Attempting to create a `BifurcationProblem` for a non-autonomous system (e.g. where some rate depend on $(rs.iv)). This is not possible.")
+    end
 
     # Converts symbols to symbolics.
     (bif_par isa Symbol) && (bif_par = ModelingToolkit.get_var_to_name(rs)[bif_par])

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -36,6 +36,9 @@ Notes:
 ```
   """
 function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=[], kwargs...)
+    if !isautonomous(rs) 
+        error("Attempting to compute steady state for a non-autonomous system (e.g. where some rate depend on $(rs.iv)). This is not possible.")
+    end
     ss_poly = steady_state_polynomial(rs, ps, u0)
     sols = HC.real_solutions(HC.solve(ss_poly; kwargs...))
     reorder_sols!(sols, ss_poly, rs)
@@ -57,14 +60,6 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0)
     return poly_type_convert(ss_poly)
 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
-
 # 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)

src/Catalyst.jl

@@ -18,7 +18,7 @@ using DynamicQuantities#, Unitful # Having Unitful here as well currently gives
 
 @reexport using ModelingToolkit
 using Symbolics
-
+using LinearAlgebra
 using RuntimeGeneratedFunctions
 RuntimeGeneratedFunctions.init(@__MODULE__)
 
@@ -102,6 +102,7 @@ export species, nonspecies, reactions, nonreactions, speciesmap, paramsmap
 export numspecies, numreactions, setdefaults!
 export make_empty_network
 export dependants, dependents, substoichmat, prodstoichmat, netstoichmat
+export isautonomous
 export reactionrates
 export isequivalent
 export set_default_noise_scaling
@@ -149,6 +150,10 @@ export @compound, @compounds
 export iscompound, components, coefficients, component_coefficients
 export balance_reaction, balance_system
 
+# Functionality for computing the stability of system steady states.
+include("steady_state_stability.jl")
+export steady_state_stability, steady_state_jac
+
 ### Extensions ###
 
 # HomotopyContinuation

src/reactionsystem.jl

@@ -1161,6 +1161,36 @@ function dependants(rx, network)
     dependents(rx, network)
 end
 
+"""
+isautonomous(rs::ReactionSystem)
+
+Checks if a system is autonomous (i.e. no rate or equation depend on the independent variable(s)).
+Example:
+```julia
+rs1 = @reaction_system
+    (p,d), 0 <--> X
+end
+isautonomous(rs1) # Returns `true`.
+
+rs2 = @reaction_system
+    (p/t,d), 0 <--> X
+end
+isautonomous(rs2) # Returns `false`.
+```
+"""
+function isautonomous(rs::ReactionSystem)
+    # Get all variables occurring in reactions and equations.
+    vars = Set()
+    for eq in equations(rs)
+        (eq isa Reaction) ? get_variables!(vars, eq.rate) : get_variables!(vars, eq)
+    end
+
+    # Checks for iv and spatial ivs
+    (get_iv(rs) in vars) && return false
+    any(in(vars), get_sivs(rs)) && return false
+    return true
+end
+
 
 ### `ReactionSystem` Remaking ###
 

src/reactionsystem_conversions.jl

@@ -518,8 +518,14 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
                       include_zero_odes = true, remove_conserved = false, checks = false,
                       default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
                       all_differentials_permitted = false, kwargs...)
+    # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, NonlinearSystem)
+    if !isautonomous(rs) 
+        error("Attempting to convert a non-autonomous `ReactionSystem` (e.g. where some rate depend on $(rs.iv)) to a `NonlinearSystem`. This is not possible. if you are intending to compute system steady states, consider creating and solving a `SteadyStateProblem.")
+    end
+
+    # Generates system equations.
     fullrs = Catalyst.flatten(rs)
     remove_conserved && conservationlaws(fullrs)
     ists, ispcs = get_indep_sts(fullrs, remove_conserved)

src/steady_state_stability.jl

@@ -0,0 +1,134 @@
+### Stability Analysis ###
+
+"""
+    stability(u::Vector{T}, rs::ReactionSystem, p; tol = 10*sqrt(eps())
+                ss_jac = steady_state_jac(u, rs, p))
+
+Compute the stability of a steady state (Returned as a `Bool`, with `true` indicating stability).
+
+Arguments:
+- `u`: The steady state for which we want to compute stability.
+- `rs`: The `ReactionSystem` model for which we want to compute stability.
+- `p`: The parameter set for which we want to compute stability.
+- `tol = 10*sqrt(eps())`: The tolerance used for determining whether computed eigenvalues real 
+parts can reliably be considered negative/positive. In practise, a steady state is considered 
+stable if the corresponding Jacobian's maximum eigenvalue real part is < 0. However, if this maximum 
+eigenvalue is in the range `-tol< eig < tol`, and error is throw, as we do not deem that stability
+can be ensured with enough certainty. The choice `tol = 10*sqrt(eps())` has *not* been subject
+to much analysis.
+- `ss_jac = steady_state_jac(u, rs)`: It is possible to pre-compute the 
+Jacobian used for stability computation using `steady_state_jac`. If stability is computed 
+for many states, precomputing the Jacobian may speed up evaluation.
+
+Example:
+```julia
+# Creates model.
+rn = @reaction_network begin
+    (p,d), 0 <--> X
+end
+p = [:p => 1.0, :d => 0.5]
+
+# Finds (the trivial) steady state, and computes its stability.
+steady_state = [2.0]
+steady_state_stability(steady_state, rn, p)
+
+Notes:
+- The input `u` can (currently) be both a vector of values (like `[1.0, 2.0]`) and a vector map 
+(like `[X => 1.0, Y => 2.0]`). The reason is that most methods for computing stability only
+produces vectors of values. However, if possible, it is recommended to work with values on the
+form of maps.
+- Catalyst currently computes steady state stabilities using the naive approach of checking whether
+a system's largest eigenvalue real part is negative. While more advanced stability computation 
+methods exist (and would be a welcome addition to Catalyst), there is no direct plans to implement 
+these. Furthermore, Catalyst uses a tolerance `tol = 10*sqrt(eps())` to determine whether a 
+computed eigenvalue is far away enough from 0 to be reliably used. This selected threshold can be changed through the `tol` argument.
+```
+"""
+function steady_state_stability(u::Vector, rs::ReactionSystem, ps; tol = 10*sqrt(eps(ss_val_type(u))),
+                                ss_jac = steady_state_jac(rs; u0 = u))
+    # Warning checks.
+    if !isautonomous(rs) 
+        error("Attempting to compute stability for a non-autonomous system (e.g. where some rate depend on $(rs.iv)). This is not possible.")
+    end
+
+    # If `u` is a vector of values, we convert it to a map. Also, if there are conservation laws,
+    # corresponding values must be eliminated from `u`.
+    u = steady_state_u_conversion(u, rs)
+    if length(u) > length(unknowns(ss_jac.f.sys))
+        u = filter(u_v -> any(isequal(u_v[1]), unknowns(ss_jac.f.sys)), u)
+    end
+
+    # Generates the Jacobian at the steady state (technically, `ss_jac` is an `ODEProblem` with dummy values for `u0` and `p`).
+    J = zeros(length(u), length(u))
+    ss_jac = remake(ss_jac; u0 = u, p = ps)
+    ss_jac.f.jac(J, ss_jac.u0, ss_jac.p, Inf)
+    
+    # Computes stability (by checking that the real part of all eigenvalues is negative).
+    max_eig = maximum(real(ev) for ev in eigvals(J))
+    if abs(max_eig) < tol
+        error("The system Jacobian's maximum eigenvalue at the steady state is within the tolerance range (abs($max_eig) < $tol). Hence, stability could not be reliably determined. If you still wish to compute the stability, reduce the `tol` argument, but note that floating point error in the eigenvalue calculation could lead to incorrect results.")
+    end
+    return max_eig < 0
+end
+
+# Used to determine the type of the steady states values, which is then used to set the tolerance's
+# type.
+ss_val_type(u::Vector{T}) where {T} = T
+ss_val_type(u::Vector{Pair{S,T}}) where {S,T} = T
+ss_val_type(u::Dict{S,T}) where {S,T} = T
+
+"""
+    steady_state_jac(rs::ReactionSystem; u0 = [])
+
+Creates the Jacobian function which can be used as input to `steady_state_stability`. Useful when 
+a large number of stability computation has to be carried out in a performant manner.
+
+Arguments:
+    - `rs`: The reaction system model for which we want to compute stability.
+    - `u0 = []`: For systems with conservation laws, a `u` is required to compute the conserved quantities.
+ 
+Example:
+```julia
+# Creates model.
+rn = @reaction_network begin
+    (p,d), 0 <--> X
+end
+p = [:p => 1.0, :d => 0.5]
+
+# Creates the steady state jacobian.
+ss_jac = steady_state_jacobian(rn)
+
+# Finds (the trivial) steady state, and computes stability.
+steady_state = [2.0]
+steady_state_stability(steady_state, rn, p; ss_jac)
+
+Notes:
+    - In practise, this function returns an `ODEProblem` (with a computed Jacobian) set up in 
+    such a way that it can be used by the `steady_state_stability` function.
+```
+"""
+function steady_state_jac(rs::ReactionSystem; u0 = [sp => 0.0 for sp in unknowns(rs)], 
+                          combinatoric_ratelaws = get_combinatoric_ratelaws(rs))
+    # If u0 is a vector of values, must be converted to something MTK understands.
+
+    # Converts u0 to values MTK understands, and checks that potential conservation laws are accounted for.
+    u0 = steady_state_u_conversion(u0, rs)
+    conservationlaw_errorcheck(rs, u0)
+
+    # Creates an `ODEProblem` with a Jacobian. Dummy values for `u0` and `ps` must be provided.
+    ps = [p => 0.0 for p in parameters(rs)]
+    return ODEProblem(rs, u0, 0, ps; jac = true, remove_conserved = true, 
+                      combinatoric_ratelaws = combinatoric_ratelaws)
+end
+
+# Converts a `u` vector from a vector of values to a map.
+function steady_state_u_conversion(u, rs::ReactionSystem)
+    if (u isa Vector{<:Number})
+        if length(u) == length(unknowns(rs))
+            u = [sp => v for (sp,v) in zip(unknowns(rs), u)]
+        else
+            error("You are trying to generate a stability Jacobian, providing u0 to compute conservation laws. Your provided u0 vector has length < the number of system states. If you provide a u0 vector, these have to be identical.")
+        end
+    end
+    return symmap_to_varmap(rs, u)
+end