Skip to content

Commit

Permalink
Merge pull request #963 from vyudu/linkageclasses
Browse files Browse the repository at this point in the history
Strong and terminal linkage classes
  • Loading branch information
isaacsas authored Jul 12, 2024
2 parents 58f5f15 + 28be9eb commit 8cfa3a1
Show file tree
Hide file tree
Showing 4 changed files with 142 additions and 1 deletion.
3 changes: 2 additions & 1 deletion src/Catalyst.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
50 changes: 50 additions & 0 deletions src/network_analysis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
6 changes: 6 additions & 0 deletions src/reactionsystem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
84 changes: 84 additions & 0 deletions test/network_analysis/network_properties.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit 8cfa3a1

Please sign in to comment.