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, deficiency, subnetworks
+export complexoutgoingmat, incidencematgraph, linkageclasses, stronglinkageclasses,
+       terminallinkageclasses, deficiency, subnetworks
 export linkagedeficiencies, isreversible, isweaklyreversible
 export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants
 

src/network_analysis.jl

@@ -348,6 +348,56 @@ end
 
 linkageclasses(incidencegraph) = Graphs.connected_components(incidencegraph)
 
+"""
+    stronglinkageclasses(rn::ReactionSystem)
+
+    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)
+    nps = get_networkproperties(rn)
+    if isempty(nps.stronglinkageclasses)
+        nps.stronglinkageclasses = stronglinkageclasses(incidencematgraph(rn))
+    end
+    nps.stronglinkageclasses
+end
+
+stronglinkageclasses(incidencegraph) = Graphs.strongly_connected_components(incidencegraph)
+
+"""
+    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 outgoing reaction from a complex in the component produces a complex also in the component).
+"""
+
+function terminallinkageclasses(rn::ReactionSystem)
+    nps = get_networkproperties(rn)
+    if isempty(nps.terminallinkageclasses)
+        slcs = stronglinkageclasses(rn)
+        tslcs = filter(lc -> isterminal(lc, rn), slcs)
+        nps.terminallinkageclasses = tslcs
+    end
+    nps.terminallinkageclasses
+end
+
+# 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
+        end
+    end
+    true
+end
+
 @doc raw"""
     deficiency(rn::ReactionSystem)
 

src/reactionsystem.jl

@@ -78,6 +78,7 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     isempty::Bool = true
     netstoichmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0)
     conservationmat::Matrix{I} = Matrix{I}(undef, 0, 0)
+    cyclemat::Matrix{I} = Matrix{I}(undef, 0, 0)
     col_order::Vector{Int} = Int[]
     rank::Int = 0
     nullity::Int = 0
@@ -93,6 +94,8 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     complexoutgoingmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0)
     incidencegraph::Graphs.SimpleDiGraph{Int} = Graphs.DiGraph()
     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
 end
 #! format: on
@@ -116,6 +119,7 @@ function reset!(nps::NetworkProperties{I, V}) where {I, V}
     nps.isempty && return
     nps.netstoichmat = Matrix{Int}(undef, 0, 0)
     nps.conservationmat = Matrix{I}(undef, 0, 0)
+    nps.cyclemat = Matrix{Int}(undef, 0, 0)
     empty!(nps.col_order)
     nps.rank = 0
     nps.nullity = 0
@@ -131,6 +135,8 @@ function reset!(nps::NetworkProperties{I, V}) where {I, V}
     nps.complexoutgoingmat = Matrix{Int}(undef, 0, 0)
     nps.incidencegraph = Graphs.DiGraph()
     empty!(nps.linkageclasses)
+    empty!(nps.stronglinkageclasses)
+    empty!(nps.terminallinkageclasses)
     nps.deficiency = 0
 
     # this needs to be last due to setproperty! setting it to false

test/network_analysis/network_properties.jl

@@ -325,3 +325,87 @@ let
     rates = Dict(zip(parameters(rn), k))
     @test Catalyst.iscomplexbalanced(rn, rates) == true 
 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
+        k3, B + C --> D
+        k4, D --> E
+        (k5, k6), E <--> 2F
+        k7, 2F --> D
+        (k8, k9), D + E <--> G
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    slcs = stronglinkageclasses(rn)
+    tslcs = terminallinkageclasses(rn)
+    @test length(slcs) == 3
+    @test length(tslcs) == 2
+    @test issubset([[1,2], [3,4,5], [6,7]], slcs)
+    @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
+        k3, B + C --> D
+        k4, D --> E
+        (k5, k6), E <--> 2F
+        k7, 2F --> D
+        (k8, k9), D + E --> G
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    slcs = stronglinkageclasses(rn)
+    tslcs = terminallinkageclasses(rn)
+    @test length(slcs) == 4
+    @test length(tslcs) == 2
+    @test issubset([[1,2], [3,4,5], [6], [7]], slcs)
+    @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
+        (k3, k4), B + C <--> D
+        k5, D --> E
+        (k6, k7), E <--> 2F
+        k8, 2F --> D
+        (k9, k10), D + E <--> G
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    slcs = stronglinkageclasses(rn)
+    tslcs = terminallinkageclasses(rn)
+    @test length(slcs) == 2
+    @test length(tslcs) == 2
+    @test issubset([[1,2,3,4,5], [6,7]], slcs)
+    @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 
+        k3, A --> C + D
+        (k4, k5), C + D <--> E
+        k6, 2B --> F
+        (k7, k8), F <--> 2G
+        (k9, k10), 2G <--> H
+        k11, H --> F
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    slcs = stronglinkageclasses(rn)
+    tslcs = terminallinkageclasses(rn)
+    @test length(slcs) == 3
+    @test length(tslcs) == 2
+    @test issubset([[1,2], [3,4], [5,6,7]], slcs)
+    @test issubset([[3,4], [5,6,7]], tslcs) 
+end