From 8b26cb7945c053cf1b457e2329f08f765b22dd9a Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 18 Jun 2024 17:55:52 -0400 Subject: [PATCH 01/22] implement robustness for deficiency one --- src/network_analysis.jl | 46 +++++++++++++++++++++ test/network_analysis/network_properties.jl | 26 ++++++++++++ 2 files changed, 72 insertions(+) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index b6f2954ce0..07b4df6a59 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -931,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 diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index baee75e32a..0f28a14ada 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -404,3 +404,29 @@ let @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 From ad5d765b9b42325fd170888c0796ce4db4f79ad2 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 19 Jun 2024 12:18:22 -0400 Subject: [PATCH 02/22] JUliaFormatter --- src/Catalyst.jl | 3 ++- src/network_analysis.jl | 25 +++++++++++++------------ 2 files changed, 15 insertions(+), 13 deletions(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index cec9cca52f..84b1a31b71 100644 --- a/src/Catalyst.jl +++ b/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 diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 07b4df6a59..f1dbca2186 100644 --- a/src/network_analysis.jl +++ b/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)) @@ -381,17 +381,17 @@ stronglinkageclasses(incidencegraph) = Graphs.strongly_connected_components(inci 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) +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) +function isterminal(lc::Vector, rn::ReactionSystem) imat = incidencemat(rn) for r in 1:size(imat, 2) @@ -940,40 +940,41 @@ end 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) +function satisfiesdeficiencyone(rn::ReactionSystem) complexes, D = reactioncomplexes(rn) δ = deficiency(rn) δ_l = linkagedeficiencies(rn) - lcs = linkageclasses(rn); tslcs = terminallinkageclasses(rn) + lcs = linkageclasses(rn) + tslcs = terminallinkageclasses(rn) all(<=(1), δ_l) && sum(δ_l) == δ && length(lcs) == length(tslcs) end #function deficiencyonealgorithm() - + #end -function robustspecies(rn::ReactionSystem) +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]) + 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) +function isconcentrationrobust(rn::ReactionSystem, species::Int) robust_species = robustspecies(rn) return species in robust_species end From d9635fc075a80e3d1a005748614038d669bfbe7f Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 21 Jun 2024 11:47:23 -0400 Subject: [PATCH 03/22] added deficiency one theorem test --- src/network_analysis.jl | 6 +-- test/network_analysis/network_properties.jl | 56 +++++++++++++++++++++ 2 files changed, 57 insertions(+), 5 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index f1dbca2186..77dc1cd9b0 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -951,15 +951,11 @@ function satisfiesdeficiencyone(rn::ReactionSystem) 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.") + error("This algorithm currently only checks for robust species in networks with deficiency one.") end tslcs = terminallinkageclasses(rn) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 0f28a14ada..ca510f8d12 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -430,3 +430,59 @@ let @test Catalyst.robustspecies(EnvZ_OmpR) == [6] end + +### DEFICIENCY ONE TESTS + + +let # fails because there are two terminal linkage classes in the linkage class + rn = @reaction_network begin + k1, A + B --> 2B + k2, A + B --> 2A + end + + @test Catalyst.satisfiesdeficiencyone(rn) == false +end + + +let # fails because linkage deficiencies do not sum to total deficiency + rn = @reaction_network begin + (k1, k2), A <--> 2A + (k3, k4), A + B <--> C + (k5, k6), C <--> B + end + + @test Catalyst.satisfiesdeficiencyone(rn) == false +end + +let # fails because a linkage class has deficiency two + rn = @reaction_network begin + k1, 3A --> A + 2B + k2, A + 2B --> 3B + k3, 3B --> 2A + B + k4, 2A + B --> 3A + end + + @test Catalyst.satisfiesdeficiencyone(rn) == false +end + +let + rn = @reaction_network begin + (k1, k2), 2A <--> D + (k3, k4), D <--> A + B + (k5, k6), A + B <--> C + (k7, k8), C <--> 2B + (k9, k10), C + D <--> E + F + (k11, k12), E + F <--> H + (k13, k14), H <--> C + E + (k15, k16), C + E <--> D + F + (k17, k18), A + D <--> G + (k19, k20), G <--> B + H + end + + @test Catalyst.satisfiesdeficiencyone(rn) == true +end + + + + + From db8c410871aa5d876b4b61be6e2a266b596ff4bc Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 21 Jun 2024 15:51:54 -0400 Subject: [PATCH 04/22] cov --- src/network_analysis.jl | 23 +++++++++++++-------- src/reactionsystem.jl | 1 + test/network_analysis/network_properties.jl | 2 ++ 3 files changed, 17 insertions(+), 9 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 77dc1cd9b0..8b26fb4313 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -953,21 +953,26 @@ end function robustspecies(rn::ReactionSystem) complexes, D = reactioncomplexes(rn) - + nps = get_networkproperties(rn) + if deficiency(rn) != 1 error("This algorithm currently 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[] + if isempty(nps.robustspecies) + 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...) + 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 + nps.robustspecies = robust_species end - robust_species + + nps.robustspecies end function isconcentrationrobust(rn::ReactionSystem, species::Int) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index c45231f502..1c96d7a860 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -96,6 +96,7 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re 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) + robustspecies::Vector{Int} = Vector{Int}(undef, 0) deficiency::Int = 0 end #! format: on diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index ca510f8d12..e597cc4577 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -416,6 +416,7 @@ let end @test Catalyst.robustspecies(IDHKP_IDH) == [2] + @test Catalyst.isconcentrationrobust(IDHKP_IDH, 2) == true end let @@ -429,6 +430,7 @@ let end @test Catalyst.robustspecies(EnvZ_OmpR) == [6] + @test Catalyst.isconcentrationrobust(EnvZ_OmpR, 6) == true end ### DEFICIENCY ONE TESTS From 414c7c54ca42fe9bdca3a57aaaa89c6f04e011c0 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 21 Jun 2024 17:23:54 -0400 Subject: [PATCH 05/22] juliaformatter --- src/network_analysis.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 8b26fb4313..fbbfcfb4f5 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -954,12 +954,12 @@ end function robustspecies(rn::ReactionSystem) complexes, D = reactioncomplexes(rn) nps = get_networkproperties(rn) - + if deficiency(rn) != 1 error("This algorithm currently only checks for robust species in networks with deficiency one.") end - if isempty(nps.robustspecies) + if isempty(nps.robustspecies) tslcs = terminallinkageclasses(rn) Z = complexstoichmat(rn) nonterminal_complexes = deleteat!([1:length(complexes);], vcat(tslcs...)) From d429a8c41cd690bd4ddd096eecfba50bee8bd36d Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 25 Jun 2024 13:47:44 -0400 Subject: [PATCH 06/22] commented --- src/network_analysis.jl | 35 +++++++++++++++++---- src/reactionsystem.jl | 4 ++- test/network_analysis/network_properties.jl | 13 +++++--- 3 files changed, 40 insertions(+), 12 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index a0e4110bab..53d387557d 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -437,11 +437,15 @@ rcs,incidencemat = reactioncomplexes(sir) """ function deficiency(rn::ReactionSystem) nps = get_networkproperties(rn) - conservationlaws(rn) - r = nps.rank - ig = incidencematgraph(rn) - lc = linkageclasses(rn) - nps.deficiency = Graphs.nv(ig) - length(lc) - r + + # Check if deficiency has been computeud already (initialized to -1) + if nps.deficiency == -1 + conservationlaws(rn) + r = nps.rank + ig = incidencematgraph(rn) + lc = linkageclasses(rn) + nps.deficiency = Graphs.nv(ig) - length(lc) - r + end nps.deficiency end @@ -953,9 +957,17 @@ function satisfiesdeficiencyone(rn::ReactionSystem) lcs = linkageclasses(rn) tslcs = terminallinkageclasses(rn) + # Check the conditions for the deficiency one theorem: 1) the deficiency of each individual linkage class is at most 1; 2) the sum of the linkage deficiencies is the total deficiency, and 3) there is only one terminal linkage class per linkage class. all(<=(1), δ_l) && sum(δ_l) == δ && length(lcs) == length(tslcs) end + +""" + robustspecies(rn::ReactionSystem) + + Return a vector of indices corresponding to species that are concentration robust, i.e. for every positive equilbrium, the concentration of species s will be the same. +""" + function robustspecies(rn::ReactionSystem) complexes, D = reactioncomplexes(rn) nps = get_networkproperties(rn) @@ -964,22 +976,33 @@ function robustspecies(rn::ReactionSystem) error("This algorithm currently only checks for robust species in networks with deficiency one.") end - if isempty(nps.robustspecies) + # A species is concnetration robust in a deficiency one network if there are two non-terminal complexes (i.e. complexes belonging to a linkage class that is not terminal) that differ only in species s (i.e. their difference is some multiple of s. (A + B, A) differ only in B. (A + 2B, B) differ in both A and B, since A + 2B - B = A + B). + if !nps.checkedrobust tslcs = terminallinkageclasses(rn) Z = complexstoichmat(rn) + # Find the complexes that do not belong to a terminal linkage class nonterminal_complexes = deleteat!([1:length(complexes);], vcat(tslcs...)) robust_species = Int64[] for (c_s, c_p) in collect(Combinatorics.combinations(nonterminal_complexes, 2)) + # Check the difference of all the combinations of complexes. The support is the set of indices that are non-zero supp = findall(!=(0), Z[:, c_s] - Z[:, c_p]) + # If the support has length one, then they differ in only one species, and that species is concentration robust. length(supp) == 1 && supp[1] ∉ robust_species && push!(robust_species, supp...) end + nps.checkedrobust = true nps.robustspecies = robust_species end nps.robustspecies end +""" + isconcentrationrobust(rn::ReactionSystem, species::Int) + + Given a reaction network and an index of a species, check if that species is concentration robust. +""" + function isconcentrationrobust(rn::ReactionSystem, species::Int) robust_species = robustspecies(rn) return species in robust_species diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 1c96d7a860..12e0cd7281 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -96,8 +96,10 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re 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) + + checkedrobust::Bool = false robustspecies::Vector{Int} = Vector{Int}(undef, 0) - deficiency::Int = 0 + deficiency::Int = -1 end #! format: on diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index c169393cb0..c308eebbfd 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -412,6 +412,8 @@ end ### CONCENTRATION ROBUSTNESS TESTS +# Check whether concentration-robust species are correctly identified for two well-known reaction networks: the glyoxylate IDHKP-IDH system, and the EnvZ_OmpR signaling pathway. + let IDHKP_IDH = @reaction_network begin (k1, k2), EIp + I <--> EIpI @@ -440,8 +442,8 @@ end ### DEFICIENCY ONE TESTS - -let # fails because there are two terminal linkage classes in the linkage class +# Fails because there are two terminal linkage classes in the linkage class +let rn = @reaction_network begin k1, A + B --> 2B k2, A + B --> 2A @@ -450,8 +452,8 @@ let # fails because there are two terminal linkage classes in the linkage class @test Catalyst.satisfiesdeficiencyone(rn) == false end - -let # fails because linkage deficiencies do not sum to total deficiency +# Fails because linkage deficiencies do not sum to total deficiency +let rn = @reaction_network begin (k1, k2), A <--> 2A (k3, k4), A + B <--> C @@ -461,7 +463,8 @@ let # fails because linkage deficiencies do not sum to total deficiency @test Catalyst.satisfiesdeficiencyone(rn) == false end -let # fails because a linkage class has deficiency two +# Fails because a linkage class has deficiency two +let rn = @reaction_network begin k1, 3A --> A + 2B k2, A + 2B --> 3B From 5eb54c74c9934588732f7f69067119e98b199a77 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 25 Jun 2024 13:51:43 -0400 Subject: [PATCH 07/22] formatter --- src/network_analysis.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 53d387557d..82a5c03e76 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -391,7 +391,6 @@ function terminallinkageclasses(rn::ReactionSystem) nps.terminallinkageclasses end - # Check whether a given linkage class in a reaction network is terminal, i.e. all outgoing reactions from complexes in the linkage class produce a complex also in hte linkage class function isterminal(lc::Vector, rn::ReactionSystem) imat = incidencemat(rn) @@ -961,7 +960,6 @@ function satisfiesdeficiencyone(rn::ReactionSystem) all(<=(1), δ_l) && sum(δ_l) == δ && length(lcs) == length(tslcs) end - """ robustspecies(rn::ReactionSystem) @@ -977,7 +975,7 @@ function robustspecies(rn::ReactionSystem) end # A species is concnetration robust in a deficiency one network if there are two non-terminal complexes (i.e. complexes belonging to a linkage class that is not terminal) that differ only in species s (i.e. their difference is some multiple of s. (A + B, A) differ only in B. (A + 2B, B) differ in both A and B, since A + 2B - B = A + B). - if !nps.checkedrobust + if !nps.checkedrobust tslcs = terminallinkageclasses(rn) Z = complexstoichmat(rn) # Find the complexes that do not belong to a terminal linkage class @@ -1001,7 +999,7 @@ end isconcentrationrobust(rn::ReactionSystem, species::Int) Given a reaction network and an index of a species, check if that species is concentration robust. -""" +""" function isconcentrationrobust(rn::ReactionSystem, species::Int) robust_species = robustspecies(rn) From 9b77f7d172c254a5407f7da8f5d2289fb60386fe Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 11 Jul 2024 11:05:37 -0400 Subject: [PATCH 08/22] comments --- src/network_analysis.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 82a5c03e76..38fbbed0d8 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -437,7 +437,7 @@ rcs,incidencemat = reactioncomplexes(sir) function deficiency(rn::ReactionSystem) nps = get_networkproperties(rn) - # Check if deficiency has been computeud already (initialized to -1) + # Check if deficiency has been computed already (initialized to -1) if nps.deficiency == -1 conservationlaws(rn) r = nps.rank @@ -956,8 +956,11 @@ function satisfiesdeficiencyone(rn::ReactionSystem) lcs = linkageclasses(rn) tslcs = terminallinkageclasses(rn) - # Check the conditions for the deficiency one theorem: 1) the deficiency of each individual linkage class is at most 1; 2) the sum of the linkage deficiencies is the total deficiency, and 3) there is only one terminal linkage class per linkage class. - all(<=(1), δ_l) && sum(δ_l) == δ && length(lcs) == length(tslcs) + # Check the conditions for the deficiency one theorem: + # 1) the deficiency of each individual linkage class is at most 1; + # 2) the sum of the linkage deficiencies is the total deficiency, and + # 3) there is only one terminal linkage class per linkage class. + all(<=(1), δ_l) && (sum(δ_l) == δ) && (length(lcs) == length(tslcs)) end """ @@ -974,10 +977,14 @@ function robustspecies(rn::ReactionSystem) error("This algorithm currently only checks for robust species in networks with deficiency one.") end - # A species is concnetration robust in a deficiency one network if there are two non-terminal complexes (i.e. complexes belonging to a linkage class that is not terminal) that differ only in species s (i.e. their difference is some multiple of s. (A + B, A) differ only in B. (A + 2B, B) differ in both A and B, since A + 2B - B = A + B). + # A species is concnetration robust in a deficiency one network if there are two non-terminal complexes (i.e. complexes + # belonging to a linkage class that is not terminal) that differ only in species s (i.e. their difference is some + # multiple of s. (A + B, A) differ only in B. (A + 2B, B) differ in both A and B, since A + 2B - B = A + B). + if !nps.checkedrobust tslcs = terminallinkageclasses(rn) Z = complexstoichmat(rn) + # Find the complexes that do not belong to a terminal linkage class nonterminal_complexes = deleteat!([1:length(complexes);], vcat(tslcs...)) robust_species = Int64[] @@ -985,6 +992,7 @@ function robustspecies(rn::ReactionSystem) for (c_s, c_p) in collect(Combinatorics.combinations(nonterminal_complexes, 2)) # Check the difference of all the combinations of complexes. The support is the set of indices that are non-zero supp = findall(!=(0), Z[:, c_s] - Z[:, c_p]) + # If the support has length one, then they differ in only one species, and that species is concentration robust. length(supp) == 1 && supp[1] ∉ robust_species && push!(robust_species, supp...) end From 75f2ed55f9fc3af6952fd7a5735dbbdfcdda2eb7 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 11 Jul 2024 11:07:30 -0400 Subject: [PATCH 09/22] Formatter --- src/network_analysis.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 38fbbed0d8..867e868294 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -980,7 +980,7 @@ function robustspecies(rn::ReactionSystem) # A species is concnetration robust in a deficiency one network if there are two non-terminal complexes (i.e. complexes # belonging to a linkage class that is not terminal) that differ only in species s (i.e. their difference is some # multiple of s. (A + B, A) differ only in B. (A + 2B, B) differ in both A and B, since A + 2B - B = A + B). - + if !nps.checkedrobust tslcs = terminallinkageclasses(rn) Z = complexstoichmat(rn) From dfa3c577e7b0b6a8407844b09b96be8b1f5f34db Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 16 Jul 2024 14:40:33 -0400 Subject: [PATCH 10/22] reset fix --- src/reactionsystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 8c342693e4..5d08228e15 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -140,7 +140,7 @@ function reset!(nps::NetworkProperties{I, V}) where {I, V} empty!(nps.linkageclasses) empty!(nps.stronglinkageclasses) empty!(nps.terminallinkageclasses) - nps.deficiency = 0 + nps.deficiency = -1 # this needs to be last due to setproperty! setting it to false nps.isempty = true From 78ad6472e1045cec8c230358b127208dc7b28696 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 17 Jul 2024 14:39:03 -0400 Subject: [PATCH 11/22] formatter --- src/reactionsystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 5d08228e15..324b946637 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -140,7 +140,7 @@ function reset!(nps::NetworkProperties{I, V}) where {I, V} empty!(nps.linkageclasses) empty!(nps.stronglinkageclasses) empty!(nps.terminallinkageclasses) - nps.deficiency = -1 + nps.deficiency = -1 # this needs to be last due to setproperty! setting it to false nps.isempty = true From a89d09135fb993304c0c7a8aedcce7ba99fbf02c Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 18 Jul 2024 11:23:42 -0400 Subject: [PATCH 12/22] deficiencyzero and doc update --- src/network_analysis.jl | 24 +++++++++--------- test/network_analysis/network_properties.jl | 27 +++++++++++++++++++++ 2 files changed, 40 insertions(+), 11 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index ed380e1e80..d62469e034 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -933,10 +933,23 @@ function satisfiesdeficiencyone(rn::ReactionSystem) all(<=(1), δ_l) && (sum(δ_l) == δ) && (length(lcs) == length(tslcs)) end +""" + satisfiesdeficiencyzero(rn::ReactionSystem) + + Check if a reaction network obeys the conditions of the deficiency zero theorem, which ensures that there is only one equilibrium for every positive stoichiometric compatibility class, this equilibrium is asymptotically stable, and this equilibrium is complex balanced. +""" + +function satisfiesdeficiencyzero(rn::ReactionSystem) + δ = deficiency(rn) + δ == 0 && isweaklyreversible(rn) +end + """ robustspecies(rn::ReactionSystem) Return a vector of indices corresponding to species that are concentration robust, i.e. for every positive equilbrium, the concentration of species s will be the same. + + Note: This function currently only works for networks of deficiency one, and is not currently guaranteed to return *all* the concentration-robust species in the network. Any species returned by the function will be robust, but this may not include all of them. Use with caution. Support for higher deficiency networks and necessary conditions for robustness will be coming in the future. """ function robustspecies(rn::ReactionSystem) @@ -972,14 +985,3 @@ function robustspecies(rn::ReactionSystem) nps.robustspecies end - -""" - isconcentrationrobust(rn::ReactionSystem, species::Int) - - Given a reaction network and an index of a species, check if that species is concentration robust. -""" - -function isconcentrationrobust(rn::ReactionSystem, species::Int) - robust_species = robustspecies(rn) - return species in robust_species -end diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index d2282a64e9..94c0b58e86 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -624,7 +624,34 @@ let @test Catalyst.satisfiesdeficiencyone(rn) == true end +### Some tests for deficiency zero networks. +let + rn = @reaction_network begin + (k1, k2), A <--> 2B + (k3, k4), A + C <--> D + k5, D --> B + E + k6, B + E --> A + C + end + # No longer weakly reversible + rn2 = @reaction_network begin + (k1, k2), A <--> 2B + (k3, k4), A + C <--> D + k5, B + E --> D + k6, B + E --> A + C + end + + # No longer weakly reversible + rn3 = @reaction_network begin + k1, A --> 2B + (k3, k4), A + C <--> D + k5, B + E --> D + k6, B + E --> A + C + end + @test satisfiesdeficiencyzero(rn) == true + @test satisfiesdeficiencyzero(rn2) == false + @test satisfiesdeficiencyzero(rn3) == false +end From 83228b3690f834c3df0db9fb1245a31e1095047d Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 18 Jul 2024 11:24:35 -0400 Subject: [PATCH 13/22] Formatter --- src/network_analysis.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index ace7ec1fef..67601e4d92 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -972,7 +972,7 @@ end Check if a reaction network obeys the conditions of the deficiency zero theorem, which ensures that there is only one equilibrium for every positive stoichiometric compatibility class, this equilibrium is asymptotically stable, and this equilibrium is complex balanced. """ -function satisfiesdeficiencyzero(rn::ReactionSystem) +function satisfiesdeficiencyzero(rn::ReactionSystem) δ = deficiency(rn) δ == 0 && isweaklyreversible(rn) end From c5914c793b5527a14860027bc62487f15736e0e1 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 18 Jul 2024 11:26:55 -0400 Subject: [PATCH 14/22] fix test --- test/network_analysis/network_properties.jl | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 033efce953..bbd52115ec 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -713,12 +713,20 @@ let rn3 = @reaction_network begin k1, A --> 2B (k3, k4), A + C <--> D - k5, B + E --> D + k5, D --> B + E k6, B + E --> A + C end + # Weakly reversible but deficiency one + rn4 = @reaction_network begin + (k1, k2), A <--> 2A + (k3, k4), A + B <--> C + (k5, k6), C <--> B + end + @test satisfiesdeficiencyzero(rn) == true @test satisfiesdeficiencyzero(rn2) == false @test satisfiesdeficiencyzero(rn3) == false + @test satisfiesdeficiencyzero(rn4) == false end From 2dfb133ea7b76f9281e1d7a13f70e0379e49f6d9 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 18 Jul 2024 11:30:00 -0400 Subject: [PATCH 15/22] fix test --- test/network_analysis/network_properties.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index bbd52115ec..c909517d54 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -622,7 +622,6 @@ let end @test Catalyst.robustspecies(IDHKP_IDH) == [2] - @test Catalyst.isconcentrationrobust(IDHKP_IDH, 2) == true end let @@ -636,7 +635,6 @@ let end @test Catalyst.robustspecies(EnvZ_OmpR) == [6] - @test Catalyst.isconcentrationrobust(EnvZ_OmpR, 6) == true end ### DEFICIENCY ONE TESTS From b1774d0f50cdc487fd75871729e93542c3e66b57 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 18 Jul 2024 11:54:23 -0400 Subject: [PATCH 16/22] fix test --- test/network_analysis/network_properties.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index c909517d54..decdbaa0ac 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -722,9 +722,9 @@ let (k5, k6), C <--> B end - @test satisfiesdeficiencyzero(rn) == true - @test satisfiesdeficiencyzero(rn2) == false - @test satisfiesdeficiencyzero(rn3) == false - @test satisfiesdeficiencyzero(rn4) == false + @test Catalyst.satisfiesdeficiencyzero(rn) == true + @test Catalyst.satisfiesdeficiencyzero(rn2) == false + @test Catalyst.satisfiesdeficiencyzero(rn3) == false + @test Catalyst.satisfiesdeficiencyzero(rn4) == false end From afbfc94a1a2de8c50f2350ed92e20faa38758a66 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 18 Jul 2024 14:41:03 -0400 Subject: [PATCH 17/22] test fix --- src/network_analysis.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 67601e4d92..7ee76ee2a8 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -974,7 +974,7 @@ end function satisfiesdeficiencyzero(rn::ReactionSystem) δ = deficiency(rn) - δ == 0 && isweaklyreversible(rn) + δ == 0 && isweaklyreversible(rn, subnetworks(rn)) end """ From 29f6dd23121bb94cc6bf101b3d56e307c7fc4917 Mon Sep 17 00:00:00 2001 From: vyudu Date: Thu, 18 Jul 2024 16:15:48 -0400 Subject: [PATCH 18/22] up --- src/network_analysis.jl | 11 ++++++----- src/reactionsystem.jl | 2 ++ 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 7ee76ee2a8..3b06db7bbd 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -963,7 +963,7 @@ function satisfiesdeficiencyone(rn::ReactionSystem) # 1) the deficiency of each individual linkage class is at most 1; # 2) the sum of the linkage deficiencies is the total deficiency, and # 3) there is only one terminal linkage class per linkage class. - all(<=(1), δ_l) && (sum(δ_l) == δ) && (length(lcs) == length(tslcs)) + all(<=(1), δ_l) && (sum(δ_l) == δ) && (length(lcs) == length(tslcs)) end """ @@ -973,8 +973,9 @@ end """ function satisfiesdeficiencyzero(rn::ReactionSystem) + all(r -> ismassaction(r), reactions(rn)) || error("The deficiency one theorem is only valid for reaction networks that are mass action.") δ = deficiency(rn) - δ == 0 && isweaklyreversible(rn, subnetworks(rn)) + δ == 0 && isweaklyreversible(rn, subnetworks(rn)) && all(r -> ismassaction(r)) end """ @@ -993,7 +994,7 @@ function robustspecies(rn::ReactionSystem) error("This algorithm currently only checks for robust species in networks with deficiency one.") end - # A species is concnetration robust in a deficiency one network if there are two non-terminal complexes (i.e. complexes + # A species is concentration robust in a deficiency one network if there are two non-terminal complexes (i.e. complexes # belonging to a linkage class that is not terminal) that differ only in species s (i.e. their difference is some # multiple of s. (A + B, A) differ only in B. (A + 2B, B) differ in both A and B, since A + 2B - B = A + B). @@ -1005,9 +1006,9 @@ function robustspecies(rn::ReactionSystem) nonterminal_complexes = deleteat!([1:length(complexes);], vcat(tslcs...)) robust_species = Int64[] - for (c_s, c_p) in collect(Combinatorics.combinations(nonterminal_complexes, 2)) + for (c_s, c_p) in Combinatorics.combinations(nonterminal_complexes, 2) # Check the difference of all the combinations of complexes. The support is the set of indices that are non-zero - supp = findall(!=(0), Z[:, c_s] - Z[:, c_p]) + supp = findall(!=(0), @views Z[:, c_s] - Z[:, c_p]) # If the support has length one, then they differ in only one species, and that species is concentration robust. length(supp) == 1 && supp[1] ∉ robust_species && push!(robust_species, supp...) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 324b946637..7497544c24 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -141,6 +141,8 @@ function reset!(nps::NetworkProperties{I, V}) where {I, V} empty!(nps.stronglinkageclasses) empty!(nps.terminallinkageclasses) nps.deficiency = -1 + empty!(nps.robustspecies) + nps.checkedrobust = false # this needs to be last due to setproperty! setting it to false nps.isempty = true From 9dce9ed22a992dcb37eba4b9b1de35090a30c730 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 19 Jul 2024 15:46:44 -0400 Subject: [PATCH 19/22] up --- src/network_analysis.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 3b06db7bbd..f3712ddd40 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -952,6 +952,7 @@ end """ function satisfiesdeficiencyone(rn::ReactionSystem) + all(r -> ismassaction(r, rn), reactions(rn)) || error("The deficiency one theorem is only valid for reaction networks that are mass action.") complexes, D = reactioncomplexes(rn) δ = deficiency(rn) δ_l = linkagedeficiencies(rn) @@ -973,9 +974,9 @@ end """ function satisfiesdeficiencyzero(rn::ReactionSystem) - all(r -> ismassaction(r), reactions(rn)) || error("The deficiency one theorem is only valid for reaction networks that are mass action.") + all(r -> ismassaction(r), reactions(rn)) || error("The deficiency zero theorem is only valid for reaction networks that are mass action.") δ = deficiency(rn) - δ == 0 && isweaklyreversible(rn, subnetworks(rn)) && all(r -> ismassaction(r)) + δ == 0 && isweaklyreversible(rn, subnetworks(rn)) && all(r -> ismassaction(r, rn), reactions(rn)) end """ From 75fd1af67249d36141a7f9191268060b48c3d23b Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 19 Jul 2024 15:48:44 -0400 Subject: [PATCH 20/22] up : --- src/network_analysis.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index f3712ddd40..9bd8fe4438 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -974,9 +974,9 @@ end """ function satisfiesdeficiencyzero(rn::ReactionSystem) - all(r -> ismassaction(r), reactions(rn)) || error("The deficiency zero theorem is only valid for reaction networks that are mass action.") + all(r -> ismassaction(r, rn), reactions(rn)) || error("The deficiency zero theorem is only valid for reaction networks that are mass action.") δ = deficiency(rn) - δ == 0 && isweaklyreversible(rn, subnetworks(rn)) && all(r -> ismassaction(r, rn), reactions(rn)) + δ == 0 && isweaklyreversible(rn, subnetworks(rn)) end """ From c1c8db45460a313c6bb89787b9804bdfdcd70eaf Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 19 Jul 2024 15:49:12 -0400 Subject: [PATCH 21/22] formatter --- src/network_analysis.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 9bd8fe4438..37df0cf089 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -952,7 +952,8 @@ end """ function satisfiesdeficiencyone(rn::ReactionSystem) - all(r -> ismassaction(r, rn), reactions(rn)) || error("The deficiency one theorem is only valid for reaction networks that are mass action.") + all(r -> ismassaction(r, rn), reactions(rn)) || + error("The deficiency one theorem is only valid for reaction networks that are mass action.") complexes, D = reactioncomplexes(rn) δ = deficiency(rn) δ_l = linkagedeficiencies(rn) @@ -964,7 +965,7 @@ function satisfiesdeficiencyone(rn::ReactionSystem) # 1) the deficiency of each individual linkage class is at most 1; # 2) the sum of the linkage deficiencies is the total deficiency, and # 3) there is only one terminal linkage class per linkage class. - all(<=(1), δ_l) && (sum(δ_l) == δ) && (length(lcs) == length(tslcs)) + all(<=(1), δ_l) && (sum(δ_l) == δ) && (length(lcs) == length(tslcs)) end """ @@ -974,9 +975,10 @@ end """ function satisfiesdeficiencyzero(rn::ReactionSystem) - all(r -> ismassaction(r, rn), reactions(rn)) || error("The deficiency zero theorem is only valid for reaction networks that are mass action.") + all(r -> ismassaction(r, rn), reactions(rn)) || + error("The deficiency zero theorem is only valid for reaction networks that are mass action.") δ = deficiency(rn) - δ == 0 && isweaklyreversible(rn, subnetworks(rn)) + δ == 0 && isweaklyreversible(rn, subnetworks(rn)) end """ From 19309dcf1661dc6cbc7e6b2a3f0b72fc8d2194d5 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 22 Jul 2024 19:03:14 -0400 Subject: [PATCH 22/22] up --- src/network_analysis.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 37df0cf089..8fb8f2a208 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -1011,10 +1011,15 @@ function robustspecies(rn::ReactionSystem) for (c_s, c_p) in Combinatorics.combinations(nonterminal_complexes, 2) # Check the difference of all the combinations of complexes. The support is the set of indices that are non-zero - supp = findall(!=(0), @views Z[:, c_s] - Z[:, c_p]) + suppcnt = 0 + supp = 0 + for i in 1:size(Z, 1) + (Z[i, c_s] != Z[i, c_p]) && (suppcnt += 1; supp = i) + (suppcnt > 1) && break + end # If the support has length one, then they differ in only one species, and that species is concentration robust. - length(supp) == 1 && supp[1] ∉ robust_species && push!(robust_species, supp...) + (suppcnt == 1) && (supp ∉ robust_species) && push!(robust_species, supp) end nps.checkedrobust = true nps.robustspecies = robust_species