init

Author: TorkelE

docs/src/api.md

@@ -166,6 +166,7 @@ speciesmap
 paramsmap
 reactionparamsmap
 isspecies
+is_autonomous
 Catalyst.isconstant
 Catalyst.isbc
 ```

docs/src/catalyst_functionality/steady_state_stability_computation.md

@@ -0,0 +1,54 @@
+# 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.
+
+## 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:
+```@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` directly to this steady state vector, receiving the stability for each:
+```@example stability_1
+steady_state_stability(steady_states, sa_loop, ps)
+```
+
+## 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 the 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
+```
+
+!!! note
+    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.

src/Catalyst.jl

@@ -42,6 +42,7 @@ import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length,
 import MacroTools, Graphs
 import Graphs: DiGraph, SimpleGraph, SimpleDiGraph, vertices, edges, add_vertices!, nv, ne
 import DataStructures: OrderedDict, OrderedSet
+import LinearAlgebra.eigvals
 import Parameters: @with_kw_noshow
 import Symbolics: occursin, wrap
 
@@ -149,6 +150,10 @@ export @compound, @compounds
 export iscompound, components, coefficients, component_coefficients
 export balance_reaction
 
+# for functions I am unsure where to best place them.
+include("steady_state_stability.jl")
+export steady_state_stability, steady_state_jac
+
 ### Extensions ###
 
 # HomotopyContinuation

src/steady_state_stability.jl

@@ -0,0 +1,100 @@
+### Stability Analysis ###
+
+"""
+    stability(u::Vector{T}, rs::ReactionSystem, p; sparse=false, ss_jac = steady_state_jac(u, rs, p; sparse=sparse), t=Inf, non_autonomous_war=true)
+
+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 reaction system model for which we want to compute stability.
+    - `p`: The parameter set for which we want to compute stability.
+    - `sparse=false`: If we wish to create a sparse Jacobian for the stability computation.
+    - `ss_jac = steady_state_jac(u, rs; sparse=sparse)`: 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.
+    - `t=Inf`: The time point at which stability is computed. 
+    - `non_autonomous_war=true`: If the system is non-autonomous (e.g. a rate depends on t), a warning will be given. Set this to false to remove that. Alternatively, specify a nonInf value for `t`.
+
+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)
+```
+
+Notes:
+    y states (each being a `Vector`) is provided as `u`, stability is computed for each state separately.
+"""
+function steady_state_stability(u::Vector{T}, rs::ReactionSystem, p; 
+                                sparse=false, ss_jac = steady_state_jac(rs; u0=u, sparse=sparse), t=Inf, non_autonomous_war=true) where T
+    # Warning checks.
+    !is_autonomous(rs) && non_autonomous_war && @warn "Attempting to compute stability for a non-autonomous system. Set `non_autonomous_war=false` to disable this warning."
+
+    # Because Jacobian currently requires ps to be a normal vector, can be removed once this get fixed in MTK.
+    if (u isa Vector{<:Pair}) || (u isa Dict) 
+        u_dict = Dict(symmap_to_varmap(rs, u))
+        u = [u_dict[var] for var in states(rs)]        
+    end
+    if (p isa Vector{<:Pair}) || (p isa Dict)
+        p_dict = Dict(symmap_to_varmap(rs, p))
+        p = [p_dict[var] for var in parameters(rs)]
+    end
+
+    # Computes stability
+    jac = ss_jac(u, p, Inf)
+    return maximum(real.(eigvals(jac))) < 0
+end
+# Computes the stability for a vector of steady states.
+function steady_state_stability(us::Vector{Vector{T}}, rs::ReactionSystem, p; sparse=false, ss_jac = steady_state_jac(rs; u0=us[1], sparse=sparse)) where T
+    return [steady_state_stability(u, rs, p) for u in us]
+end
+
+"""
+    steady_state_jac(rs::ReactionSystem; u0=[], sparse=false)
+
+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.
+ 
+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)
+```
+
+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.
+    - `sparse=false`: If we wish to create a sparse Jacobian for the stability computation.
+"""
+function steady_state_jac(rs::ReactionSystem; u0=[], sparse=false) where T
+    # If u0 is a vector of values, must be converted to something MTK understands.
+    if (u0 isa Vector{<:Number})
+        if length(u0) == length(states(rs))
+            u0 = Pair.(states(rs), u0)
+        elseif !isempty(u0)
+            error("You are trying to compute a stability matrix, 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
+
+    # Converts u0 to values MTK understands, and checks that potential conservation laws are accounted for.
+    u0 = symmap_to_varmap(rs, u0)
+    conservationlaw_errorcheck(rs, u0)
+
+    # Creates a ODESystem and jacobian.
+    osys = convert(ODESystem, rs; remove_conserved = true, defaults=Dict(u0))
+    return ODEFunction(osys; jac=true, sparse=sparse).jac
+end

test/miscellaneous_tests/api.jl

@@ -446,3 +446,51 @@ let
     neweqs = getfield.(equations(ns), :rhs)
     @test_throws MethodError Catalyst.to_multivariate_poly(neweqs)
 end
+
+# Tests `is_autonomous` function.
+let 
+    # Using default iv.
+    rn1 = @reaction_network begin
+        (p + X*(p1/(t+p3)),d), 0 <--> X
+        (kB,kD), 2X <--> X
+    end
+    rn2 = @reaction_network begin
+        (hill(X, v/t, K, n),d), 0 <--> X
+        (kB,kD), 2X <--> X
+    end
+    rn3 = @reaction_network begin
+        (p + X*(p1+p2),d), 0 <--> X
+        (kB,kD), 2X <--> X
+    end
+
+    @test !is_autonomous(rn1)
+    @test !is_autonomous(rn2)
+    @test is_autonomous(rn3)
+
+    # Using non-default iv.
+    rn4 = @reaction_network begin
+        @ivs i1 i2
+        (p + X*(p1/(1+i1)),d), 0 <--> X
+        (kB,kD), 2X <--> X
+    end
+    rn5 = @reaction_network begin
+        @ivs i1 i2
+        (p + X*(i2+p2),d), 0 <--> X
+        (kB,kD), 2X <--> X
+    end
+    rn6 = @reaction_network begin
+        @ivs i1 i2
+        (hill(X, v/i1, i2, n),d), 0 <--> X
+        (kB,kD), 2X <--> X
+    end
+    rn7 = @reaction_network begin
+        @ivs i1 i2
+        (p + X*(p1+p2),d), 0 <--> X
+        (kB,kD), 2X <--> X
+    end
+
+    @test !is_autonomous(rn4)
+    @test !is_autonomous(rn5)
+    @test !is_autonomous(rn6)
+    @test is_autonomous(rn7)
+end

test/miscellaneous_tests/stability_computation.jl

@@ -0,0 +1,79 @@
+### Fetch Packages ###
+
+# Fetch packages.
+using Catalyst, OrdinaryDiffEq
+using Random, Test
+
+# Sets rnd number.
+using StableRNGs
+rng = StableRNG(12345)
+
+### Run Tests ###
+
+# Tests that stability is correctly assessed (using simulation) in multi stable system.
+# Tests that `steady_state_jac` function works.
+# Tests with and without sparsity.
+# tests using symbolic input.
+let
+    # System which may have between 1 and 7 fixed points.
+    rn = @reaction_network begin
+        v/20.0 + hillar(X,Y,v,K,n), 0 --> X
+        v/20.0 + hillar(Y,X,v,K,n), 0 --> Y
+        d, (X,Y) --> 0
+    end
+    ss_jac = steady_state_jac(rn)
+    ss_jac_sparse = steady_state_jac(rn; sparse=true)
+
+    # Repeats several times, most cases should be encountered several times.
+    for i = 1:50
+        # Generates random parameter values (which can generate all steady states cases).
+        p = [:v => 1.0 + 3*rand(rng), :K => 0.5 + 2*rand(rng), :n => rand(rng,[1,2,3,4]), :d => 0.5 + rand(rng)]
+
+        sss = hc_steady_states(rn, p)
+        stabs_1 = steady_state_stability(sss, rn, p)
+        stabs_2 = steady_state_stability(sss, rn, p; sparse=true)
+        stabs_3 = steady_state_stability(sss, rn, p; ss_jac=ss_jac)
+        stabs_4 = steady_state_stability(sss, rn, p; ss_jac=ss_jac_sparse)
+
+        for (idx,ss) in enumerate(sss)
+            oprob = ODEProblem(rn, [1.001, 0.999] .* ss, (0.0,1000.0), p)
+            sol_end = solve(oprob, Rosenbrock23())[end]
+            stabs_5 = ss ≈ sol_end
+            @test stabs_1[idx] == stabs_2[idx] == stabs_3[idx] == stabs_4[idx] == stabs_5
+        end
+    end
+end
+
+# Checks stability for system with known stability structure.
+# Tests for system with conservation laws.
+# Tests for various input forms of u0 and ps.
+let
+    rn = complete(@reaction_network begin
+        k1+Z, Y --> 2X
+        k2, 2X --> X + Y
+        k3, X + Y --> Y
+        k4, X --> 0
+        (kD1+X, kD2), 2Z <--> Z2
+    end)
+    @unpack k1, k2, k3, k4, kD1, kD2, X, Y, Z, Z2 = rn
+    u0_1 = [X => 1.0, Y => 1.0, Z => 1.0, Z2 => 1.0]
+    u0_2 = [:X => 1.0, :Y => 1.0, :Z => 1.0, :Z2 => 1.0]
+    u0_3 = [rn.X => 1.0, rn.Y => 1.0, rn.Z => 1.0, rn.Z2 => 1.0]
+    u0_4 = [1.0, 1.0, 1.0, 1.0]
+    ps_1 = [k1 => 8.0, k2 => 2.0, k3 => 1.0, k4 => 1.5, kD1 => 0.5, kD2 => 2.0]
+    ps_2 = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5, :kD1 => 0.5, :kD2 => 2.0]
+    ps_3 = [rn.k1 => 8.0, rn.k2 => 2.0, rn.k3 => 1.0, rn.k4 => 1.5, rn.kD1 => 0.5, rn.kD2 => 4.0]
+    ps_4 = [8.0, 2.0, 1.0, 1.5, 0.5, 4.0]
+    
+    sss = hc_steady_states(rn, ps_1; u0=u0_1)
+    for u0 in [u0_1, u0_2, u0_3, u0_4], ps in [ps_1, ps_2, ps_3, ps_4]
+        stab_1 = steady_state_stability(sss, rn, ps)
+        @test length(stab_1) == 3
+        @test count(stab_1) == 2
+    
+        ss_jac = steady_state_jac(rn; u0=u0)
+        stab_2 = steady_state_stability(sss, rn, ps; ss_jac=ss_jac)
+        @test length(stab_2) == 3
+        @test count(stab_2) == 2
+    end
+end