Merge branch 'deficiency-one'

Author: vyudu

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

@@ -886,3 +886,22 @@ 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) 
+    δ = deficiency(rn)
+    δ_l = linkagedeficiencies(rn)
+
+    complexes, D = reactioncomplexes(rn)
+    lcs = linkageclasses(rn); tslcs = terminallinkageclasses(rn)
+
+    δ_l .<= 1 && sum(δ_l) = δ && length(lcs) == length(tslcs)
+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,81 @@ 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