Merge branch 'deficiency-one' into detailedbalance

src/network_analysis.jl

@@ -378,7 +378,7 @@ stronglinkageclasses(incidencegraph) = Graphs.strongly_connected_components(inci
 """
     terminallinkageclasses(rn::ReactionSystem)
 
-    Return the terminal strongly connected components of a reaction network's incidence graph (i.e. sub-groups of reaction complexes that are 1) strongly connected and 2) every reaction in the component produces a complex in the component).
+    Return the terminal strongly connected components of a reaction network's incidence graph (i.e. sub-groups of reaction complexes that are 1) strongly connected and 2) every outgoing reaction from a complex in the component produces a complex also in the component).
 """
 
 function terminallinkageclasses(rn::ReactionSystem)
@@ -391,11 +391,15 @@ function terminallinkageclasses(rn::ReactionSystem)
     nps.terminallinkageclasses
 end
 
+# Check whether a given linkage class in a reaction network is terminal, i.e. all outgoing reactions from complexes in the linkage class produce a complex also in hte linkage class
 function isterminal(lc::Vector, rn::ReactionSystem)
     imat = incidencemat(rn)
 
     for r in 1:size(imat, 2)
+        # Find the index of the reactant complex for a given reaction
         s = findfirst(==(-1), @view imat[:, r])
+
+        # If the reactant complex is in the linkage class, check whether the product complex is also in the linkage class. If any of them are not, return false. 
         if s in Set(lc)
             p = findfirst(==(1), @view imat[:, r])
             p in Set(lc) ? continue : return false
@@ -446,11 +450,15 @@ rcs,incidencemat = reactioncomplexes(sir)
 """
 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 computeud 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
 
@@ -1043,31 +1051,51 @@ function satisfiesdeficiencyone(rn::ReactionSystem)
     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
 
-#function deficiencyonealgorithm(rn::ReactionSystem) 
-#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. 
+"""
 
 function robustspecies(rn::ReactionSystem)
     complexes, D = reactioncomplexes(rn)
-    robust_species = Int64[]
+    nps = get_networkproperties(rn)
 
     if deficiency(rn) != 1
-        error("This method currently only checks for robust species in networks with deficiency one.")
+        error("This algorithm currently only checks for robust species in networks with deficiency one.")
     end
 
-    tslcs = terminallinkageclasses(rn)
-    Z = complexstoichmat(rn)
-    nonterminal_complexes = deleteat!([1:length(complexes);], vcat(tslcs...))
-
-    for (c_s, c_p) in collect(Combinatorics.combinations(nonterminal_complexes, 2))
-        supp = findall(!=(0), Z[:, c_s] - Z[:, c_p])
-        length(supp) == 1 && supp[1] ∉ robust_species && push!(robust_species, supp...)
+    # A species is concnetration 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 collect(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 
+            supp = findall(!=(0), Z[:, c_s] - Z[:, c_p])
+            # If the support has length one, then they differ in only one species, and that species is concentration robust. 
+            length(supp) == 1 && supp[1] ∉ robust_species && push!(robust_species, supp...)
+        end
+        nps.checkedrobust = true
+        nps.robustspecies = robust_species
     end
-    robust_species
+
+    nps.robustspecies
 end
 
+"""
+    isconcentrationrobust(rn::ReactionSystem, species::Int)
+
+    Given a reaction network and an index of a species, check if that species is concentration robust. 
+"""
+
 function isconcentrationrobust(rn::ReactionSystem, species::Int)
     robust_species = robustspecies(rn)
     return species in robust_species

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
 

test/network_analysis/network_properties.jl

@@ -330,6 +330,8 @@ end
 
 ### STRONG LINKAGE CLASS TESTS
 
+
+# a) Checks that strong/terminal linkage classes are correctly found. Should identify the (A, B+C) linkage class as non-terminal, since B + C produces D
 let
     rn = @reaction_network begin
         (k1, k2), A <--> B + C
@@ -349,6 +351,7 @@ let
     @test issubset([[3,4,5], [6,7]], tslcs) 
 end
 
+# b) Makes the D + E --> G reaction irreversible. Thus, (D+E) becomes a non-terminal linkage class. Checks whether correctly identifies both (A, B+C) and (D+E) as non-terminal 
 let
     rn = @reaction_network begin
         (k1, k2), A <--> B + C
@@ -368,6 +371,7 @@ let
     @test issubset([[3,4,5], [7]], tslcs) 
 end
 
+# From a), makes the B + C <--> D reaction reversible. Thus, the non-terminal (A, B+C) linkage class gets absorbed into the terminal (A, B+C, D, E, 2F) linkage class, and the terminal linkage classes and strong linkage classes coincide. 
 let
     rn = @reaction_network begin
         (k1, k2), A <--> B + C
@@ -387,6 +391,7 @@ let
     @test issubset([[1,2,3,4,5], [6,7]], tslcs) 
 end
 
+# Simple test for strong and terminal linkage classes
 let
     rn = @reaction_network begin
         (k1, k2), A <--> 2B 
@@ -409,6 +414,8 @@ 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
@@ -418,6 +425,7 @@ let
     end
 
     @test Catalyst.robustspecies(IDHKP_IDH) == [2]
+    @test Catalyst.isconcentrationrobust(IDHKP_IDH, 2) == true
 end
 
 let
@@ -431,8 +439,66 @@ let
     end
 
     @test Catalyst.robustspecies(EnvZ_OmpR) == [6]
+    @test Catalyst.isconcentrationrobust(EnvZ_OmpR, 6) == true
+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
+
+
+
+
+
+
 
 ### Complex and detailed balance tests