Merge pull request #822 from vyudu/master

Complex balancing

Project.toml

@@ -3,6 +3,7 @@ uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
 version = "13.5.1"
 
 [deps]
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
@@ -38,6 +39,7 @@ CatalystStructuralIdentifiabilityExtension = "StructuralIdentifiability"
 [compat]
 BifurcationKit = "0.3"
 DataStructures = "0.18"
+Combinatorics = "1.0.2"
 DiffEqBase = "6.83.0"
 DocStringExtensions = "0.8, 0.9"
 DynamicPolynomials = "0.5"

src/Catalyst.jl

@@ -6,6 +6,7 @@ module Catalyst
 using DocStringExtensions
 using SparseArrays, DiffEqBase, Reexport, Setfield
 using LaTeXStrings, Latexify, Requires
+using LinearAlgebra, Combinatorics
 using JumpProcesses: JumpProcesses, JumpProblem, 
                      MassActionJump, ConstantRateJump, VariableRateJump,
                      SpatialMassActionJump

src/network_analysis.jl

@@ -718,3 +718,164 @@ function conservationlaw_errorcheck(rs, pre_varmap)
     isempty(conservedequations(Catalyst.flatten(rs))) ||
         error("The system has conservation laws but initial conditions were not provided for some species.")
 end
+
+"""
+    iscomplexbalanced(rs::ReactionSystem, parametermap)
+
+Constructively compute whether a network will have complex-balanced equilibrium
+solutions, following the method in van der Schaft et al., [2015](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
+"""
+
+function iscomplexbalanced(rs::ReactionSystem, parametermap::Dict) 
+    if length(parametermap) != numparams(rs) 
+        error("Incorrect number of parameters specified.")
+    end
+
+    pmap = symmap_to_varmap(rs, parametermap)
+    pmap = Dict(ModelingToolkit.value(k) => v for (k,v) in pmap)
+
+    sm = speciesmap(rs)
+    cm = reactioncomplexmap(rs)
+    complexes, D = reactioncomplexes(rs)
+    rxns = reactions(rs)
+    nc = length(complexes); nr = numreactions(rs); nm = numspecies(rs)
+
+    if !all(r->ismassaction(r, rs), rxns) 
+        error("The supplied ReactionSystem has reactions that are not ismassaction. Testing for being complex balanced is currently only supported for pure mass action networks.") 
+    end
+
+    rates = [substitute(rate, pmap) for rate in reactionrates(rs)]
+
+    # Construct kinetic matrix, K
+    K = zeros(nr, nc) 
+    for c in 1:nc
+        complex = complexes[c]
+        for (r, dir) in cm[complex]
+            rxn = rxns[r]
+            if dir == -1
+                K[r, c] = rates[r]
+            end
+        end
+    end
+
+    L = -D*K
+    S = netstoichmat(rs)
+
+    # Compute ρ using the matrix-tree theorem
+    g = incidencematgraph(rs)
+    R = ratematrix(rs, rates)
+    ρ = matrixtree(g, R)
+
+    # Determine if 1) ρ is positive and 2) D^T Ln ρ lies in the image of S^T
+    if all(>(0), ρ)
+        img = D'*log.(ρ)
+        if rank(S') == rank(hcat(S', img)) return true else return false end 
+    else
+        return false
+    end
+end
+
+function iscomplexbalanced(rs::ReactionSystem, parametermap::Vector{Pair{Symbol, Float64}}) 
+    pdict = Dict(parametermap)
+    iscomplexbalanced(rs, pdict)
+end
+
+function iscomplexbalanced(rs::ReactionSystem, parametermap::Tuple{Pair{Symbol, Float64}}) 
+    pdict = Dict(parametermap)
+    iscomplexbalanced(rs, pdict)
+end
+
+iscomplexbalanced(rs::ReactionSystem, parametermap) = error("Parameter map must be a dictionary, tuple, or vector of symbol/value pairs.")
+
+
+"""
+    ratematrix(rs::ReactionSystem, parametermap)
+
+    Given a reaction system with n complexes, outputs an n-by-n matrix where R_{ij} is the rate constant of the reaction between complex i and complex j. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
+"""
+
+function ratematrix(rs::ReactionSystem, rates::Vector{Float64}) 
+    complexes, D = reactioncomplexes(rs)
+    n = length(complexes)
+    rxns = reactions(rs)
+    ratematrix = zeros(n, n)
+
+    for r in 1:length(rxns)
+        rxn = rxns[r]
+        s = findfirst(==(-1), @view D[:,r])
+        p = findfirst(==(1), @view D[:,r])
+        ratematrix[s, p] = rates[r]
+    end
+    ratematrix
+end
+
+function ratematrix(rs::ReactionSystem, parametermap::Dict) 
+    if length(parametermap) != numparams(rs) 
+        error("Incorrect number of parameters specified.")
+    end
+
+    pmap = symmap_to_varmap(rs, parametermap)
+    pmap = Dict(ModelingToolkit.value(k) => v for (k,v) in pmap)
+
+    rates = [substitute(rate, pmap) for rate in reactionrates(rs)]
+    ratematrix(rs, rates)
+end
+
+function ratematrix(rs::ReactionSystem, parametermap::Vector{Pair{Symbol, Float64}}) 
+    pdict = Dict(parametermap)
+    ratematrix(rs, pdict)
+end
+
+function ratematrix(rs::ReactionSystem, parametermap::Tuple{Pair{Symbol, Float64}}) 
+    pdict = Dict(parametermap)
+    ratematrix(rs, pdict)
+end
+
+ratematrix(rs::ReactionSystem, parametermap) = error("Parameter map must be a dictionary, tuple, or vector of symbol/value pairs.")
+
+### BELOW: Helper functions for iscomplexbalanced
+
+function matrixtree(g::SimpleDiGraph, distmx::Matrix)
+    n = nv(g)
+    if size(distmx) != (n, n)
+        error("Size of distance matrix is incorrect")
+    end
+
+    π = zeros(n) 
+
+    if !Graphs.is_connected(g)
+        ccs = Graphs.connected_components(g)
+        for cc in ccs
+            sg, vmap = Graphs.induced_subgraph(g, cc)
+            distmx_s = distmx[cc, cc]
+            π_j = matrixtree(sg, distmx_s)
+            π[cc] = π_j
+        end
+        return π
+    end
+
+    # generate all spanning trees
+    ug = SimpleGraph(SimpleDiGraph(g))
+    trees = collect(Combinatorics.combinations(collect(edges(ug)), n-1))
+    trees = SimpleGraph.(trees)
+    trees = filter!(t->isempty(Graphs.cycle_basis(t)), trees)
+
+    # constructed rooted trees for every vertex, compute sum
+    for v in 1:n
+        rootedTrees = [reverse(Graphs.bfs_tree(t, v, dir=:in)) for t in trees]
+        π[v] = sum([treeweight(t, g, distmx) for t in rootedTrees])
+    end
+
+    # sum the contributions
+    return π 
+end
+
+function treeweight(t::SimpleDiGraph, g::SimpleDiGraph, distmx::Matrix) 
+    prod = 1
+    for e in edges(t)
+        s = Graphs.src(e); t = Graphs.dst(e)
+        prod *= distmx[s, t] 
+    end
+    prod 
+end
+

test/network_analysis/network_properties.jl

@@ -1,7 +1,9 @@
 ### Prepares Tests ###
 
 # Fetch packages.
-using Catalyst, LinearAlgebra, Test
+using Catalyst, LinearAlgebra, Test, StableRNGs
+
+rng = StableRNG(514)
 
 ### Basic Tests ###
 
@@ -40,6 +42,10 @@ let
     @test isweaklyreversible(MAPK, subnetworks(MAPK)) == false
     cls = conservationlaws(MAPK)
     @test Catalyst.get_networkproperties(MAPK).rank == 15
+
+    k = rand(rng, numparams(MAPK))
+    rates = Dict(zip(parameters(MAPK), k))
+    @test Catalyst.iscomplexbalanced(MAPK, rates) == false
     # i=0;
     # for lcs in linkageclasses(MAPK)
     #     i=i+1
@@ -77,6 +83,10 @@ let
     @test isweaklyreversible(rn2, subnetworks(rn2)) == false
     cls = conservationlaws(rn2)
     @test Catalyst.get_networkproperties(rn2).rank == 6
+
+    k = rand(rng, numparams(rn2))
+    rates = Dict(zip(parameters(rn2), k))
+    @test Catalyst.iscomplexbalanced(rn2, rates) == false
     # i=0;
     # for lcs in linkageclasses(rn2)
     #     i=i+1
@@ -117,6 +127,10 @@ let
     @test isweaklyreversible(rn3, subnetworks(rn3)) == false
     cls = conservationlaws(rn3)
     @test Catalyst.get_networkproperties(rn3).rank == 10
+
+    k = rand(rng, numparams(rn3))
+    rates = Dict(zip(parameters(rn3), k))
+    @test Catalyst.iscomplexbalanced(rn3, rates) == false
     # i=0;
     # for lcs in linkageclasses(rn3)
     #     i=i+1
@@ -132,6 +146,18 @@ let
     # end
 end
 
+let
+    rn4 = @reaction_network begin
+        (k1, k2), C1 <--> C2
+        (k3, k4), C2 <--> C3
+        (k5, k6), C3 <--> C1
+    end
+
+    k = rand(rng, numparams(rn4))
+    rates = Dict(zip(parameters(rn4), k))
+    @test Catalyst.iscomplexbalanced(rn4, rates) == true
+end
+
 ### Tests Reversibility ###
 
 # Test function.
@@ -154,7 +180,12 @@ let
     rev = false
     weak_rev = false
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == false 
 end
+
 let
     rn = @reaction_network begin
         (k2, k1), A1 <--> A2 + A3
@@ -167,6 +198,10 @@ let
     rev = false
     weak_rev = false
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == false 
 end
 let
     rn = @reaction_network begin
@@ -176,6 +211,9 @@ let
     rev = false
     weak_rev = false
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == false 
 end
 let
     rn = @reaction_network begin
@@ -186,17 +224,25 @@ let
     rev = false
     weak_rev = false
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == false 
 end
 let
     rn = @reaction_network begin
         (k2, k1), A <--> 2B
-        (k4, k3), A + C --> D
+        (k4, k3), A + C <--> D
         k5, D --> B + E
         k6, B + E --> A + C
     end
     rev = false
     weak_rev = true
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == true 
 end
 let
     rn = @reaction_network begin
@@ -206,6 +252,10 @@ let
     rev = false
     weak_rev = false
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == false 
 end
 let
     rn = @reaction_network begin
@@ -215,12 +265,20 @@ let
     rev = true
     weak_rev = true
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == true  
 end
 let
     rn = @reaction_network begin (k2, k1), A + B <--> 2A end
     rev = true
     weak_rev = true
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == true 
 end
 let
     rn = @reaction_network begin
@@ -232,6 +290,10 @@ let
     rev = false
     weak_rev = true
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == true 
 end
 let
     rn = @reaction_network begin
@@ -243,4 +305,24 @@ let
     rev = false
     weak_rev = false
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
-end
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == false 
+end
+
+let
+    rn = @reaction_network begin
+        k1, 3A + 2B --> 3C 
+        k2, B + 4D --> 2E
+        k3, 2E --> 3C
+        (k4, k5), B + 4D <--> 3A + 2B
+        k6, F --> B + 4D
+        k7, 3C --> F
+    end
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == true 
+end
+