Merge remote-tracking branch 'vyudu/detailedbalance' into detailedbalance

Author: vyudu

test/network_analysis/network_properties.jl

@@ -499,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).