Refactored complex balance to not use MetaGraphs

Author: vyudu

src/networkapi.jl

@@ -1654,33 +1654,35 @@ function validate(rs::ReactionSystem, info::String = "")
 end
 
 """
-    complexbalanced(rs::ReactionSystem, rates::Vector)
+    iscomplexbalanced(rs::ReactionSystem, rates::Vector)
 
 Constructively compute whether a network will have complex-balanced equilibrium
-solutions, following the method in [this paper](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). 
+solutions, following the method in [this paper](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). Accepts a map of rates [k1 => 1.0, k2 => 2.0,...]k2
 """
 
-using MetaGraphs, Combinatorics, LinearAlgebra
+using Combinatorics, LinearAlgebra
 
-function complexbalanced(rs::ReactionSystem, rates::Vector) 
+function iscomplexbalanced(rs::ReactionSystem, rates::Dict{Any, Float64}) 
     if length(rates) != numparams(rs) 
         error("The number of reaction rates must be equal to the number of parameters")
     end
-    rxm = Dict(zip(reactionparams(rs), rates))
+
     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("Not mass action") end
+
     # 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] = rxm[rxn.rate]
+                K[r, c] = rates[rxn.rate]
             end
         end
     end
@@ -1689,93 +1691,122 @@ function complexbalanced(rs::ReactionSystem, rates::Vector)
     S = netstoichmat(rs)
 
     # Compute ρ using the matrix-tree theorem
-    ρ = zeros(nc)
-    lcs = linkageclasses(rs)
-    rwg = rateweightedgraph(rs, rates)
+    g = incidencematgraph(rs); R = ratematrix(rs, rates)
+    ρ = matrixtree(g, R)
+    @assert isapprox(L*ρ, zeros(nc), atol=1e-12) 
 
-    for lc in lcs
-        sg, vmap = Graphs.induced_subgraph(rwg, lc)
-        ρ_j = matrixtree(sg)
-        ρ[lc] = ρ_j
-    end
-    
     # Determine if 1) ρ is positive and 2) D^T Ln ρ lies in the image of S^T
     if all(x -> x > 0, ρ)
         img = D'*log.(ρ)
-        if rank(S') == rank(hcat(S', img))
-            return true
-        else
-            return false
-        end
+        if rank(S') == rank(hcat(S', img)) return true else return false end 
     else
         return false
     end
 end
 
-"""
-    rateweightedgraph(rs::ReactionSystem, rates::Vector)
-
-Generate an annotated reaction complex graph of a reaction system, where the nodes are annotated with the reaction complex they correspond to and the edges are annotated with the reaction they correspond to and the rate of the reaction. 
-"""
-function rateweightedgraph(rs::ReactionSystem, rates::Vector) 
+# """
+#     rateweightedgraph(rs::ReactionSystem, rates::Vector)
+# 
+# Generate an annotated reaction complex graph of a reaction system, where the nodes are annotated with the reaction complex they correspond to and the edges are annotated with the reaction they correspond to and the rate of the reaction. 
+# """
+# 
+# function rateweightedgraph(rs::ReactionSystem, rates::Dict{Any, Float64}) 
+#     if length(rates) != numparams(rs) 
+#         error("The number of reaction rates must be equal to the number of parameters")
+#     end
+# 
+#     complexes, D = reactioncomplexes(rs)
+#     rxns = reactions(rs)
+# 
+#     g = incidencematgraph(rs)
+#     rwg = MetaDiGraph(g)
+# 
+#     for v in vertices(rwg)
+#         set_prop!(rwg, v, :complex, complexes[v])
+#     end
+# 
+#     for (i, e) in collect(enumerate(edges(rwg)))
+#         rxn = rxns[i]
+#         set_prop!(rwg, Graphs.src(e), Graphs.dst(e), :reaction, rxn)
+#         set_prop!(rwg, Graphs.src(e), Graphs.dst(e), :rate, rates[rxn.rate])
+#     end
+# 
+#     rwg
+# end
+
+function ratematrix(rs::ReactionSystem, rates::Dict{Any, Float64}) 
     if length(rates) != numparams(rs) 
         error("The number of reaction rates must be equal to the number of parameters")
     end
-    rm = Dict(zip(reactionparams(rs), rates))
 
     complexes, D = reactioncomplexes(rs)
+    n = length(complexes)
     rxns = reactions(rs)
+    ratematrix = zeros(n, n)
 
-    g = incidencematgraph(rs)
-    rwg = MetaDiGraph(g)
-
-    for v in vertices(rwg)
-        set_prop!(rwg, v, :complex, complexes[v])
+    for r in 1:length(rxns)
+        rxn = rxns[r]
+        s = findfirst(x->x==-1, D[:,r])
+        p = findfirst(x->x==1, D[:,r])
+        ratematrix[s, p] = rates[rxn.rate]
     end
+    ratematrix
+end
 
-    for (i, e) in collect(enumerate(edges(rwg)))
-        rxn = rxns[i]
-        set_prop!(rwg, Graphs.src(e), Graphs.dst(e), :reaction, rxn)
-        set_prop!(rwg, Graphs.src(e), Graphs.dst(e), :rate, rm[rxn.rate])
+"""
+"""
+function matrixtree(g::SimpleDiGraph, distmx::Matrix)
+    n = nv(g)
+    if size(distmx) != (n, n)
+        error("Size of distance matrix is incorrect")
     end
 
-    rwg
-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
 
-function matrixtree(g::MetaDiGraph)
     # generate all spanning trees
-    # TODO: implement Winter's algorithm for generating spanning trees 
-    n = nv(g)
     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)
-    
-    π = zeros(n) 
-
-    function treeweight(t::SimpleDiGraph) 
-        prod = 1
-        for e in edges(t)
-            rate = Graphs.has_edge(g, Graphs.src(e), Graphs.dst(e)) ? get_prop(g, e, :rate) : 0
-            prod *= rate
-        end
-        prod 
-    end
+    # trees = spanningtrees(g)
 
     # constructed rooted trees for every edge, compute sum
     for v in 1:n
         rootedTrees = [reverse(Graphs.bfs_tree(t, v, dir=:in)) for t in trees]
-        π[v] = sum([treeweight(t) for t in rootedTrees])
+        π[v] = sum([treeweight(t, g, distmx) for t in rootedTrees])
     end
 
     # sum the contributions
     return π 
 end
 
-function massactionrate(rs::ReactionSystem, rxn_idx::Int, dir::Int, conc::Vector) 
-    if dir != -1 && dir != 1
-        error("Direction must be either +1 in the case of production or -1 in the case of consumption")
+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
+
+# TODO: implement Winter's algorithm for generating spanning trees 
+function spanningtrees(g::SimpleGraph) 
+        
+end
+
+# Checks if a unit consist of exponents with base 1 (and is this unitless).
+unitless_exp(u) = istree(u) && (operation(u) == ^) && (arguments(u)[1] == 1)
 
     rxn = reactions(rs)[rxn_idx]
     sm = speciesmap(rs)

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 ###
 
@@ -41,8 +43,9 @@ let
     cls = conservationlaws(MAPK)
     @test Catalyst.get_networkproperties(MAPK).rank == 15
 
-    rates = rand(numparams(MAPK))
-    @test Catalyst.complexbalanced(MAPK, rates) == false
+    k = rand(rng, numparams(MAPK))
+    rates = Dict(zip(reactionparams(MAPK), k))
+    @test Catalyst.iscomplexbalanced(MAPK, rates) == false
     # i=0;
     # for lcs in linkageclasses(MAPK)
     #     i=i+1
@@ -81,8 +84,9 @@ let
     cls = conservationlaws(rn2)
     @test Catalyst.get_networkproperties(rn2).rank == 6
 
-    rates = rand(numparams(rn2))
-    @test Catalyst.complexbalanced(rn2, rates) == false
+    k = rand(rng, numparams(rn2))
+    rates = Dict(zip(reactionparams(rn2), k))
+    @test Catalyst.iscomplexbalanced(rn2, rates) == false
     # i=0;
     # for lcs in linkageclasses(rn2)
     #     i=i+1
@@ -123,9 +127,10 @@ let
     @test isweaklyreversible(rn3, subnetworks(rn3)) == false
     cls = conservationlaws(rn3)
     @test Catalyst.get_networkproperties(rn3).rank == 10
-    
-    rates = rand(numparams(rn3))
-    @test Catalyst.complexbalanced(rn3, rates) == false
+
+    k = rand(rng, numparams(rn3))
+    rates = Dict(zip(reactionparams(rn3), k))
+    @test Catalyst.iscomplexbalanced(rn3, rates) == false
     # i=0;
     # for lcs in linkageclasses(rn3)
     #     i=i+1
@@ -141,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(reactionparams(rn4), k))
+    @test Catalyst.iscomplexbalanced(rn4, rates) == true
+end
+    
 ### Tests Reversibility ###
 
 # Test function.
@@ -163,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(reactionparams(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == false 
 end
+
 let
     rn = @reaction_network begin
         (k2, k1), A1 <--> A2 + A3
@@ -176,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(reactionparams(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == false 
 end
 let
     rn = @reaction_network begin
@@ -185,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(reactionparams(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == false 
 end
 let
     rn = @reaction_network begin
@@ -195,6 +224,10 @@ let
     rev = false
     weak_rev = false
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(reactionparams(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == false 
 end
 let
     rn = @reaction_network begin
@@ -206,6 +239,11 @@ let
     rev = false
     weak_rev = true
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    # Breaks when a reaction has multiple rates
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(reactionparams(rn), k))
+    # @test Catalyst.iscomplexbalanced(rn, rates) == true 
 end
 let
     rn = @reaction_network begin
@@ -215,6 +253,10 @@ let
     rev = false
     weak_rev = false
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(reactionparams(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == false 
 end
 let
     rn = @reaction_network begin
@@ -224,12 +266,20 @@ let
     rev = true
     weak_rev = true
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(reactionparams(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(reactionparams(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == true 
 end
 let
     rn = @reaction_network begin
@@ -241,6 +291,10 @@ let
     rev = false
     weak_rev = true
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(reactionparams(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == true 
 end
 let
     rn = @reaction_network begin
@@ -252,4 +306,8 @@ let
     rev = false
     weak_rev = false
     testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev)
+
+    k = rand(rng, numparams(rn))
+    rates = Dict(zip(reactionparams(rn), k))
+    @test Catalyst.iscomplexbalanced(rn, rates) == false 
 end