Strong and terminal linkage classes

src/Catalyst.jl

@@ -127,7 +127,8 @@ export @reaction_network, @network_component, @reaction, @species
 include("network_analysis.jl")
 export reactioncomplexmap, reactioncomplexes, incidencemat
 export complexstoichmat
-export complexoutgoingmat, incidencematgraph, linkageclasses, stronglinkageclasses, terminallinkageclasses, deficiency, subnetworks
+export complexoutgoingmat, incidencematgraph, linkageclasses, stronglinkageclasses,
+       terminallinkageclasses, deficiency, subnetworks
 export linkagedeficiencies, isreversible, isweaklyreversible
 export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants
 

src/network_analysis.jl

@@ -365,7 +365,7 @@ linkageclasses(incidencegraph) = Graphs.connected_components(incidencegraph)
     Return the strongly connected components of a reaction network's incidence graph (i.e. sub-groups of reaction complexes such that every complex is reachable from every other one in the sub-group).
 """
 
-function stronglinkageclasses(rn::ReactionSystem) 
+function stronglinkageclasses(rn::ReactionSystem)
     nps = get_networkproperties(rn)
     if isempty(nps.stronglinkageclasses)
         nps.stronglinkageclasses = stronglinkageclasses(incidencematgraph(rn))
@@ -378,24 +378,29 @@ 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) 
+function terminallinkageclasses(rn::ReactionSystem)
     nps = get_networkproperties(rn)
     if isempty(nps.terminallinkageclasses)
         slcs = stronglinkageclasses(rn)
-        tslcs = filter(lc->isterminal(lc, rn), slcs)
+        tslcs = filter(lc -> isterminal(lc, rn), slcs)
         nps.terminallinkageclasses = tslcs
     end
     nps.terminallinkageclasses
 end
 
-function isterminal(lc::Vector, rn::ReactionSystem) 
+# Helper function for terminallinkageclasses. Given a linkage class and a reaction network, say whether the linkage class is terminal, 
+# i.e. all outgoing reactions from complexes in the linkage class produce a complex also in the 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

test/network_analysis/network_properties.jl

@@ -328,6 +328,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
@@ -347,6 +349,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
@@ -366,6 +369,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
@@ -385,6 +389,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