Merge 173df48 into 4024e67

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"
@@ -12,10 +13,12 @@ LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+MetaGraphs = "626554b9-1ddb-594c-aa3c-2596fe9399a5"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -37,8 +40,8 @@ DataStructures = "0.18"
 DiffEqBase = "6.83.0"
 DocStringExtensions = "0.8, 0.9"
 Graphs = "1.4"
-JumpProcesses = "9.3.2"
 HomotopyContinuation = "2.9"
+JumpProcesses = "9.3.2"
 LaTeXStrings = "1.3.0"
 Latexify = "0.14, 0.15, 0.16"
 MacroTools = "0.5.5"

docs/src/catalyst_applications/chemical_master_equation.md

@@ -0,0 +1,17 @@
+Solving the chemical master equation using finite state projection. 
+
+
+
+```@example
+using Catalyst
+import FiniteStateProjection
+
+
+```
+
+---
+## [Citation]
+
+---
+## References
+

src/Project.toml

@@ -0,0 +1,2 @@
+[deps]
+Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"

src/networkapi.jl

@@ -1652,3 +1652,138 @@ function validate(rs::ReactionSystem, info::String = "")
 
     validated
 end
+
+"""
+    complexbalanced(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). 
+"""
+
+using MetaGraphs, Combinatorics, LinearAlgebra
+
+function complexbalanced(rs::ReactionSystem, rates::Vector) 
+    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)
+
+    # 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]
+            end
+        end
+    end
+
+    L = -D*K
+    S = netstoichmat(rs)
+
+    # Compute ρ using the matrix-tree theorem
+    ρ = zeros(nc)
+    lcs = linkageclasses(rs)
+    rwg = rateweightedgraph(rs, rates)
+
+    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
+    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) 
+    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)
+    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, rm[rxn.rate])
+    end
+
+    rwg
+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
+
+    # 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])
+    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")
+    end
+
+    rxn = reactions(rs)[rxn_idx]
+    sm = speciesmap(rs)
+    rate = rxn.rate
+
+    species = rxn.substrates
+    stoich = rxn.substoich
+    species_idx = [sm[s] for s in species]
+
+    dir * rate * prod(conc[species_idx].^stoich)
+end

test/network_analysis/network_properties.jl

@@ -37,6 +37,9 @@ let
     @test isweaklyreversible(MAPK, subnetworks(MAPK)) == false
     cls = conservationlaws(MAPK)
     @test Catalyst.get_networkproperties(MAPK).rank == 15
+
+    rates = rand(numparams(MAPK))
+    @test Catalyst.complexbalanced(MAPK, rates) == false
     # i=0;
     # for lcs in linkageclasses(MAPK)
     #     i=i+1
@@ -74,6 +77,9 @@ let
     @test isweaklyreversible(rn2, subnetworks(rn2)) == false
     cls = conservationlaws(rn2)
     @test Catalyst.get_networkproperties(rn2).rank == 6
+
+    rates = rand(numparams(rn2))
+    @test Catalyst.complexbalanced(rn2, rates) == false
     # i=0;
     # for lcs in linkageclasses(rn2)
     #     i=i+1
@@ -114,6 +120,9 @@ 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
     # i=0;
     # for lcs in linkageclasses(rn3)
     #     i=i+1
@@ -128,3 +137,14 @@ let
     #     println("-----------")
     # end
 end
+
+let
+    rn4 = @reaction_network begin
+        (k1, k2), C1 <--> C2
+        (k3, k4), C2 <--> C3
+        (k5, k6), C3 <--> C1
+    end
+    rates = rand(numparams(rn4))
+    @test Catalyst.complexbalanced(rn4, rates) == true
+end
+