Merge pull request #964 from vyudu/deficiency-one

Deficiency one and concnetration robustness

src/network_analysis.jl

@@ -421,11 +421,15 @@ end
 """
 function deficiency(rn::ReactionSystem)
     nps = get_networkproperties(rn)
-    conservationlaws(rn)
-    r = nps.rank
-    ig = incidencematgraph(rn)
-    lc = linkageclasses(rn)
-    nps.deficiency = Graphs.nv(ig) - length(lc) - r
+
+    # Check if deficiency has been computed already (initialized to -1)
+    if nps.deficiency == -1
+        conservationlaws(rn)
+        r = nps.rank
+        ig = incidencematgraph(rn)
+        lc = linkageclasses(rn)
+        nps.deficiency = Graphs.nv(ig) - length(lc) - r
+    end
     nps.deficiency
 end
 
@@ -938,3 +942,88 @@ end
 function fluxvectors(rs::ReactionSystem)
     cycles(rs)
 end
+
+### Deficiency one
+
+"""
+    satisfiesdeficiencyone(rn::ReactionSystem)
+
+    Check if a reaction network obeys the conditions of the deficiency one theorem, which ensures that there is only one equilibrium for every positive stoichiometric compatibility class.
+"""
+
+function satisfiesdeficiencyone(rn::ReactionSystem)
+    all(r -> ismassaction(r, rn), reactions(rn)) ||
+        error("The deficiency one theorem is only valid for reaction networks that are mass action.")
+    complexes, D = reactioncomplexes(rn)
+    δ = deficiency(rn)
+    δ_l = linkagedeficiencies(rn)
+
+    lcs = linkageclasses(rn)
+    tslcs = terminallinkageclasses(rn)
+
+    # Check the conditions for the deficiency one theorem: 
+    #   1) the deficiency of each individual linkage class is at most 1; 
+    #   2) the sum of the linkage deficiencies is the total deficiency, and 
+    #   3) there is only one terminal linkage class per linkage class. 
+    all(<=(1), δ_l) && (sum(δ_l) == δ) && (length(lcs) == length(tslcs))
+end
+
+"""
+    satisfiesdeficiencyzero(rn::ReactionSystem)
+
+    Check if a reaction network obeys the conditions of the deficiency zero theorem, which ensures that there is only one equilibrium for every positive stoichiometric compatibility class, this equilibrium is asymptotically stable, and this equilibrium is complex balanced.
+"""
+
+function satisfiesdeficiencyzero(rn::ReactionSystem)
+    all(r -> ismassaction(r, rn), reactions(rn)) ||
+        error("The deficiency zero theorem is only valid for reaction networks that are mass action.")
+    δ = deficiency(rn)
+    δ == 0 && isweaklyreversible(rn, subnetworks(rn))
+end
+
+"""
+    robustspecies(rn::ReactionSystem)
+
+    Return a vector of indices corresponding to species that are concentration robust, i.e. for every positive equilbrium, the concentration of species s will be the same. 
+
+    Note: This function currently only works for networks of deficiency one, and is not currently guaranteed to return *all* the concentration-robust species in the network. Any species returned by the function will be robust, but this may not include all of them. Use with caution. Support for higher deficiency networks and necessary conditions for robustness will be coming in the future.  
+"""
+
+function robustspecies(rn::ReactionSystem)
+    complexes, D = reactioncomplexes(rn)
+    nps = get_networkproperties(rn)
+
+    if deficiency(rn) != 1
+        error("This algorithm currently only checks for robust species in networks with deficiency one.")
+    end
+
+    # A species is concentration robust in a deficiency one network if there are two non-terminal complexes (i.e. complexes 
+    # belonging to a linkage class that is not terminal) that differ only in species s (i.e. their difference is some 
+    # multiple of s. (A + B, A) differ only in B. (A + 2B, B) differ in both A and B, since A + 2B - B = A + B). 
+
+    if !nps.checkedrobust
+        tslcs = terminallinkageclasses(rn)
+        Z = complexstoichmat(rn)
+
+        # Find the complexes that do not belong to a terminal linkage class 
+        nonterminal_complexes = deleteat!([1:length(complexes);], vcat(tslcs...))
+        robust_species = Int64[]
+
+        for (c_s, c_p) in Combinatorics.combinations(nonterminal_complexes, 2)
+            # Check the difference of all the combinations of complexes. The support is the set of indices that are non-zero 
+            suppcnt = 0
+            supp = 0
+            for i in 1:size(Z, 1)
+                (Z[i, c_s] != Z[i, c_p]) && (suppcnt += 1; supp = i)
+                (suppcnt > 1) && break
+            end
+
+            # If the support has length one, then they differ in only one species, and that species is concentration robust. 
+            (suppcnt == 1) && (supp ∉ robust_species) && push!(robust_species, supp)
+        end
+        nps.checkedrobust = true
+        nps.robustspecies = robust_species
+    end
+
+    nps.robustspecies
+end

src/reactionsystem.jl

@@ -96,7 +96,10 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     linkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
     stronglinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
     terminallinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
-    deficiency::Int = 0
+
+    checkedrobust::Bool = false 
+    robustspecies::Vector{Int} = Vector{Int}(undef, 0)
+    deficiency::Int = -1 
 end
 #! format: on
 
@@ -137,7 +140,9 @@ function reset!(nps::NetworkProperties{I, V}) where {I, V}
     empty!(nps.linkageclasses)
     empty!(nps.stronglinkageclasses)
     empty!(nps.terminallinkageclasses)
-    nps.deficiency = 0
+    nps.deficiency = -1
+    empty!(nps.robustspecies)
+    nps.checkedrobust = false
 
     # this needs to be last due to setproperty! setting it to false
     nps.isempty = true

test/network_analysis/network_properties.jl

@@ -608,3 +608,123 @@ let
     Catalyst.ratematrix(rn, rates_dict) == rate_mat
     @test_throws Exception Catalyst.iscomplexbalanced(rn, rates_invalid)
 end
+
+### CONCENTRATION ROBUSTNESS TESTS
+
+# Check whether concentration-robust species are correctly identified for two well-known reaction networks: the glyoxylate IDHKP-IDH system, and the EnvZ_OmpR signaling pathway. 
+
+let
+    IDHKP_IDH = @reaction_network begin
+        (k1, k2), EIp + I <--> EIpI
+        k3, EIpI --> EIp + Ip
+        (k4, k5), E + Ip <--> EIp
+        k6, EIp --> E + I
+    end
+
+    @test Catalyst.robustspecies(IDHKP_IDH) == [2]
+end
+
+let
+    EnvZ_OmpR = @reaction_network begin
+        (k1, k2), X <--> XT
+        k3, XT --> Xp
+        (k4, k5), Xp + Y <--> XpY
+        k6, XpY --> X + Yp
+        (k7, k8), XT + Yp <--> XTYp
+        k9, XTYp --> XT + Y
+    end
+
+    @test Catalyst.robustspecies(EnvZ_OmpR) == [6]
+end
+
+### DEFICIENCY ONE TESTS
+
+# Fails because there are two terminal linkage classes in the linkage class
+let 
+    rn = @reaction_network begin
+        k1, A + B --> 2B
+        k2, A + B --> 2A
+    end
+
+    @test Catalyst.satisfiesdeficiencyone(rn) == false
+end
+
+# Fails because linkage deficiencies do not sum to total deficiency
+let 
+    rn = @reaction_network begin
+        (k1, k2), A <--> 2A
+        (k3, k4), A + B <--> C
+        (k5, k6), C <--> B 
+    end
+
+    @test Catalyst.satisfiesdeficiencyone(rn) == false
+end
+
+# Fails because a linkage class has deficiency two
+let 
+    rn = @reaction_network begin
+        k1, 3A --> A + 2B
+        k2, A + 2B --> 3B
+        k3, 3B --> 2A + B
+        k4, 2A + B --> 3A
+    end
+
+    @test Catalyst.satisfiesdeficiencyone(rn) == false
+end
+
+let
+    rn = @reaction_network begin
+        (k1, k2), 2A <--> D
+        (k3, k4), D <--> A + B
+        (k5, k6), A + B <--> C
+        (k7, k8), C <--> 2B
+        (k9, k10), C + D <--> E + F
+        (k11, k12), E + F <--> H
+        (k13, k14), H <--> C + E
+        (k15, k16), C + E <--> D + F
+        (k17, k18), A + D <--> G
+        (k19, k20), G <--> B + H
+    end
+
+    @test Catalyst.satisfiesdeficiencyone(rn) == true
+end
+
+### Some tests for deficiency zero networks. 
+
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> 2B
+        (k3, k4), A + C <--> D
+        k5, D --> B + E
+        k6, B + E --> A + C
+    end
+
+    # No longer weakly reversible
+    rn2 = @reaction_network begin
+        (k1, k2), A <--> 2B
+        (k3, k4), A + C <--> D
+        k5, B + E --> D
+        k6, B + E --> A + C
+    end
+
+    # No longer weakly reversible
+    rn3 = @reaction_network begin
+        k1, A --> 2B
+        (k3, k4), A + C <--> D
+        k5, D --> B + E 
+        k6, B + E --> A + C
+    end
+
+    # Weakly reversible but deficiency one
+    rn4 = @reaction_network begin
+        (k1, k2), A <--> 2A
+        (k3, k4), A + B <--> C
+        (k5, k6), C <--> B
+    end
+
+    @test Catalyst.satisfiesdeficiencyzero(rn) == true
+    @test Catalyst.satisfiesdeficiencyzero(rn2) == false 
+    @test Catalyst.satisfiesdeficiencyzero(rn3) == false 
+    @test Catalyst.satisfiesdeficiencyzero(rn4) == false
+end
+