Merge 8b26cb7 into 134041d

src/Catalyst.jl

@@ -127,7 +127,7 @@ 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

@@ -359,6 +359,51 @@ 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 reaction in the component produces a complex 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
+
+function isterminal(lc::Vector, rn::ReactionSystem) 
+    imat = incidencemat(rn)
+
+    for r in 1:size(imat, 2)
+        s = findfirst(==(-1), @view imat[:, r])
+        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)
 
@@ -886,3 +931,49 @@ function treeweight(t::SimpleDiGraph, g::SimpleDiGraph, distmx::Matrix)
     end
     prod
 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) 
+    complexes, D = reactioncomplexes(rn)
+    δ = deficiency(rn)
+    δ_l = linkagedeficiencies(rn)
+
+    lcs = linkageclasses(rn); tslcs = terminallinkageclasses(rn)
+
+    all(<=(1), δ_l) && sum(δ_l) == δ && length(lcs) == length(tslcs)
+end
+
+#function deficiencyonealgorithm() 
+
+#end
+
+function robustspecies(rn::ReactionSystem) 
+    complexes, D = reactioncomplexes(rn)
+
+    if deficiency(rn) != 1
+        error("This algorithm 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...))
+    robust_species = Int64[]
+
+    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...)
+    end
+    robust_species
+end
+
+function isconcentrationrobust(rn::ReactionSystem, species::Int) 
+    robust_species = robustspecies(rn)
+    return species in robust_species
+end

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

test/network_analysis/network_properties.jl

@@ -326,3 +326,107 @@ let
     @test Catalyst.iscomplexbalanced(rn, rates) == true 
 end
 
+### STRONG LINKAGE CLASS TESTS
+
+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
+
+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
+
+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
+
+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
+
+### CONCENTRATION ROBUSTNESS TESTS
+
+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