Merge pull request #965 from vyudu/detailedbalance

isaacsas · web-flow · commit f1c3789c2f8a · 2024-11-06T10:08:47.000-05:00
Detailed balance
diff --git a/src/network_analysis.jl b/src/network_analysis.jl
@@ -398,6 +398,19 @@ function isterminal(lc::Vector, rn::ReactionSystem)
     true
 end
 
+function isforestlike(rn::ReactionSystem)
+    subnets = subnetworks(rn)
+    nps = get_networkproperties(rn)
+
+    isempty(nps.incidencemat) && reactioncomplexes(rn)
+    sparseig = issparse(nps.incidencemat)
+    for subnet in subnets
+        nps = get_networkproperties(subnet)
+        isempty(nps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig)
+    end
+    all(Graphs.is_tree ∘ SimpleGraph ∘ incidencematgraph, subnets)
+end
+
 @doc raw"""
     deficiency(rn::ReactionSystem)
 
@@ -724,6 +737,87 @@ function conservationlaw_errorcheck(rs, pre_varmap)
         error("The system has conservation laws but initial conditions were not provided for some species.")
 end
 
+"""
+    isdetailedbalanced(rs::ReactionSystem, parametermap; reltol=1e-9, abstol)
+
+Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit detailed-balanced equilibrium
+solutions, using the Wegscheider conditions, [Feinberg, 1989](https://www.sciencedirect.com/science/article/pii/0009250989851243). A detailed-balanced solution is one for which the rate of every forward reaction exactly equals its reverse reaction. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
+"""
+
+function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict; abstol=0, reltol=1e-9)
+    if length(parametermap) != numparams(rs)
+        error("Incorrect number of parameters specified.")
+    elseif !isreversible(rs)
+        return false
+    elseif !all(r -> ismassaction(r, rs), reactions(rs))
+        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
+
+    isforestlike(rs) && deficiency(rs) == 0 && return true
+
+    pmap = symmap_to_varmap(rs, parametermap)
+    pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap)
+
+    # Construct reaction-complex graph 
+    complexes, D = reactioncomplexes(rs)
+    img = incidencematgraph(rs)
+    undir_img = SimpleGraph(incidencematgraph(rs))
+    K = ratematrix(rs, pmap)
+
+    spanning_forest = Graphs.kruskal_mst(undir_img)
+    outofforest_edges = setdiff(collect(edges(undir_img)), spanning_forest)
+
+    # Independent Cycle Conditions: for any cycle we create by adding in an out-of-forest reaction, the product of forward reaction rates over the cycle must equal the product of reverse reaction rates over the cycle.  
+    for edge in outofforest_edges
+        g = SimpleGraph([spanning_forest..., edge])
+        ic = Graphs.cycle_basis(g)[1]
+        fwd = prod([K[ic[r], ic[r + 1]] for r in 1:(length(ic) - 1)]) * K[ic[end], ic[1]]
+        rev = prod([K[ic[r + 1], ic[r]] for r in 1:(length(ic) - 1)]) * K[ic[1], ic[end]]
+        isapprox(fwd, rev; atol = abstol, rtol = reltol) ? continue : return false
+    end
+
+    # Spanning Forest Conditions: for non-deficiency 0 networks, we get an additional δ equations. Choose an orientation for each reaction pair in the spanning forest (we will take the one given by default from kruskal_mst).  
+
+    if deficiency(rs) > 0
+        rxn_idxs = [edgeindex(D, Graphs.src(e), Graphs.dst(e)) for e in spanning_forest]
+        S_F = netstoichmat(rs)[:, rxn_idxs]
+        sols = positive_nullspace(S_F)
+
+        for i in 1:size(sols, 2)
+            α = sols[:, i]
+            fwd = prod([K[Graphs.src(e), Graphs.dst(e)]^α[i]
+                        for (e, i) in zip(spanning_forest, 1:length(α))])
+            rev = prod([K[Graphs.dst(e), Graphs.src(e)]^α[i]
+                        for (e, i) in zip(spanning_forest, 1:length(α))])
+            isapprox(fwd, rev; atol = abstol, rtol = reltol) ? continue : return false
+        end
+    end
+
+    true
+end
+
+# Helper to find the index of the reaction with a given reactant and product complex.
+function edgeindex(imat, src::T, dst::T) where T <: Int
+    for i in 1:size(imat, 2)
+        (imat[src, i] == -1) && (imat[dst, i] == 1) && return i
+    end
+    error("This edge does not exist in this reaction graph.")
+end
+
+function isdetailedbalanced(rs::ReactionSystem, parametermap::Vector{<:Pair})
+    pdict = Dict(parametermap)
+    isdetailedbalanced(rs, pdict)
+end
+
+function isdetailedbalanced(rs::ReactionSystem, parametermap::Tuple{<:Pair})
+    pdict = Dict(parametermap)
+    isdetailedbalanced(rs, pdict)
+end
+
+function isdetailedbalanced(rs::ReactionSystem, parametermap)
+    error("Parameter map must be a dictionary, tuple, or vector of symbol/value pairs.")
+end
+
 """
     iscomplexbalanced(rs::ReactionSystem, parametermap)
 
diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl
@@ -355,6 +355,7 @@ let
     k = rand(rng, numparams(rn))
     rates = Dict(zip(parameters(rn), k))
     @test Catalyst.iscomplexbalanced(rn, rates) == true 
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
 end
 
 ### STRONG LINKAGE CLASS TESTS
@@ -498,6 +499,107 @@ let
     @test isempty(cyclemat)
 end
 
+### Complex and detailed balance tests
+
+# The following network is conditionally complex balanced - it only 
+
+# Reversible, forest-like deficiency zero network - should be detailed balance for any choice of rate constants. 
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B + C
+        (k3, k4), A <--> D
+        (k5, k6), A + D <--> E
+        (k7, k8), A + D <--> G
+        (k9, k10), G <--> 2F
+        (k11, k12), A + E <--> H
+    end
+
+    k1 = rand(rng, numparams(rn))
+    rates1 = Dict(zip(parameters(rn), k1))
+    k2 = rand(StableRNG(232), numparams(rn))
+    rates2 = Dict(zip(parameters(rn), k2))
+
+    rcs, D = reactioncomplexes(rn)
+    @test Catalyst.isforestlike(rn) == true
+    @test Catalyst.isdetailedbalanced(rn, rates1) == true
+    @test Catalyst.isdetailedbalanced(rn, rates2) == true
+end
+
+# Simple connected reversible network
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B
+        (k3, k4), B <--> C
+        (k5, k6), C <--> A
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    rates1 = [:k1=>1.0, :k2=>1.0, :k3=>1.0, :k4=>1.0, :k5=>1.0, :k6=>1.0]
+    @test Catalyst.isdetailedbalanced(rn, rates1) == true
+    rates2 = [:k1=>2.0, :k2=>1.0, :k3=>1.0, :k4=>1.0, :k5=>1.0, :k6=>1.0]
+    @test Catalyst.isdetailedbalanced(rn, rates2) == false 
+end
+
+# Independent cycle tests: the following reaction entwork has 3 out-of-forest reactions. 
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B + C
+        (k3, k4), A <--> D
+        (k5, k6), B + C <--> D
+        (k7, k8), A + D <--> E
+        (k9, k10), G <--> 2F
+        (k11, k12), A + D <--> G
+        (k13, k14), G <--> E
+        (k15, k16), 2F <--> E
+        (k17, k18), A + E <--> H
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    k = rand(rng, numparams(rn))
+    p = parameters(rn)
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
+
+    # Adjust rate constants to obey the independent cycle conditions. 
+    rates[p[6]] = rates[p[1]]*rates[p[4]]*rates[p[5]] / (rates[p[2]]*rates[p[3]]) 
+    rates[p[14]] = rates[p[13]]*rates[p[11]]*rates[p[8]] / (rates[p[12]]*rates[p[7]])
+    rates[p[16]] = rates[p[8]]*rates[p[15]]*rates[p[9]]*rates[p[11]] / (rates[p[7]]*rates[p[12]]*rates[p[10]])
+    @test Catalyst.isdetailedbalanced(rn, rates) == true
+end
+
+# Deficiency two network: the following reaction network must satisfy both the independent cycle conditions and the spanning forest conditions 
+let
+    rn = @reaction_network begin
+        (k1, k2), 3A <--> A + 2B
+        (k3, k4), A + 2B <--> 3B
+        (k5, k6), 3B <--> 2A + B
+        (k7, k8), 2A + B <--> 3A
+        (k9, k10), 3A <--> 3B
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    @test Catalyst.edgeindex(D, 1, 2) == 1
+    @test Catalyst.edgeindex(D, 4, 3) == 6
+    k = rand(rng, numparams(rn))
+    p = parameters(rn)
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
+
+    # Adjust rate constants to fulfill independent cycle conditions. 
+    rates[p[8]] = rates[p[7]]*rates[p[5]]*rates[p[9]] / (rates[p[6]]*rates[p[10]])
+    rates[p[3]] = rates[p[2]]*rates[p[4]]*rates[p[9]] / (rates[p[1]]*rates[p[10]])
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
+    # Should still fail - doesn't satisfy spanning forest conditions. 
+
+    # Adjust rate constants to fulfill spanning forest conditions. 
+    cons = rates[p[6]] / rates[p[5]]
+    rates[p[1]] = rates[p[2]] * cons
+    rates[p[9]] = rates[p[10]] * cons^(3/2)
+    rates[p[8]] = rates[p[7]]*rates[p[5]]*rates[p[9]] / (rates[p[6]]*rates[p[10]])
+    rates[p[3]] = rates[p[2]]*rates[p[4]]*rates[p[9]] / (rates[p[1]]*rates[p[10]])
+    @test Catalyst.isdetailedbalanced(rn, rates) == true
+end
+
 ### Other Network Properties Tests ###
 
 # Tests outgoing complexes matrices (1). 
@@ -637,6 +739,108 @@ let
     @test Catalyst.robustspecies(EnvZ_OmpR) == [6]
 end
 
+
+### Complex and detailed balance tests
+
+# The following network is conditionally complex balanced - it only 
+
+# Reversible, forest-like deficiency zero network - should be detailed balance for any choice of rate constants. 
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B + C
+        (k3, k4), A <--> D
+        (k5, k6), A + D <--> E
+        (k7, k8), A + D <--> G
+        (k9, k10), G <--> 2F
+        (k11, k12), A + E <--> H
+    end
+
+    k1 = rand(rng, numparams(rn))
+    rates1 = Dict(zip(parameters(rn), k1))
+    k2 = rand(StableRNG(232), numparams(rn))
+    rates2 = Dict(zip(parameters(rn), k2))
+
+    rcs, D = reactioncomplexes(rn)
+    @test Catalyst.isforestlike(rn) == true
+    @test Catalyst.isdetailedbalanced(rn, rates1) == true
+    @test Catalyst.isdetailedbalanced(rn, rates2) == true
+end
+
+# Simple connected reversible network
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B
+        (k3, k4), B <--> C
+        (k5, k6), C <--> A
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    rates1 = [:k1=>1.0, :k2=>1.0, :k3=>1.0, :k4=>1.0, :k5=>1.0, :k6=>1.0]
+    @test Catalyst.isdetailedbalanced(rn, rates1) == true
+    rates2 = [:k1=>2.0, :k2=>1.0, :k3=>1.0, :k4=>1.0, :k5=>1.0, :k6=>1.0]
+    @test Catalyst.isdetailedbalanced(rn, rates2) == false 
+end
+
+# Independent cycle tests: the following reaction entwork has 3 out-of-forest reactions. 
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B + C
+        (k3, k4), A <--> D
+        (k5, k6), B + C <--> D
+        (k7, k8), A + D <--> E
+        (k9, k10), G <--> 2F
+        (k11, k12), A + D <--> G
+        (k13, k14), G <--> E
+        (k15, k16), 2F <--> E
+        (k17, k18), A + E <--> H
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    k = rand(rng, numparams(rn))
+    p = parameters(rn)
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
+
+    # Adjust rate constants to obey the independent cycle conditions. 
+    rates[p[6]] = rates[p[1]]*rates[p[4]]*rates[p[5]] / (rates[p[2]]*rates[p[3]]) 
+    rates[p[14]] = rates[p[13]]*rates[p[11]]*rates[p[8]] / (rates[p[12]]*rates[p[7]])
+    rates[p[16]] = rates[p[8]]*rates[p[15]]*rates[p[9]]*rates[p[11]] / (rates[p[7]]*rates[p[12]]*rates[p[10]])
+    @test Catalyst.isdetailedbalanced(rn, rates) == true
+end
+
+# Deficiency two network: the following reaction network must satisfy both the independent cycle conditions and the spanning forest conditions 
+let
+    rn = @reaction_network begin
+        (k1, k2), 3A <--> A + 2B
+        (k3, k4), A + 2B <--> 3B
+        (k5, k6), 3B <--> 2A + B
+        (k7, k8), 2A + B <--> 3A
+        (k9, k10), 3A <--> 3B
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    @test Catalyst.edgeindex(D, 1, 2) == 1
+    @test Catalyst.edgeindex(D, 4, 3) == 6
+    k = rand(rng, numparams(rn))
+    p = parameters(rn)
+    rates = Dict(zip(parameters(rn), k))
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
+
+    # Adjust rate constants to fulfill independent cycle conditions. 
+    rates[p[8]] = rates[p[7]]*rates[p[5]]*rates[p[9]] / (rates[p[6]]*rates[p[10]])
+    rates[p[3]] = rates[p[2]]*rates[p[4]]*rates[p[9]] / (rates[p[1]]*rates[p[10]])
+    @test Catalyst.isdetailedbalanced(rn, rates) == false
+    # Should still fail - doesn't satisfy spanning forest conditions. 
+
+    # Adjust rate constants to fulfill spanning forest conditions. 
+    cons = rates[p[6]] / rates[p[5]]
+    rates[p[1]] = rates[p[2]] * cons
+    rates[p[9]] = rates[p[10]] * cons^(3/2)
+    rates[p[8]] = rates[p[7]]*rates[p[5]]*rates[p[9]] / (rates[p[6]]*rates[p[10]])
+    rates[p[3]] = rates[p[2]]*rates[p[4]]*rates[p[9]] / (rates[p[1]]*rates[p[10]])
+    @test Catalyst.isdetailedbalanced(rn, rates) == true
+end
+
 ### DEFICIENCY ONE TESTS
 
 # Fails because there are two terminal linkage classes in the linkage class
