From 139377754ce29dcb863cb5e447e73ca134a0a705 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 19 Jun 2024 12:01:42 -0400 Subject: [PATCH 01/29] Mulitple equilibria --- src/network_analysis.jl | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 07b4df6a59..75cb80b916 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -950,21 +950,20 @@ function satisfiesdeficiencyone(rn::ReactionSystem) all(<=(1), δ_l) && sum(δ_l) == δ && length(lcs) == length(tslcs) end -#function deficiencyonealgorithm() - +#function deficiencyonealgorithm(rn::ReactionSystem) #end function robustspecies(rn::ReactionSystem) complexes, D = reactioncomplexes(rn) + robust_species = Int64[] if deficiency(rn) != 1 - error("This algorithm only checks for robust species in networks with deficiency one.") + error("This method 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[] for (c_s, c_p) in collect(Combinatorics.combinations(nonterminal_complexes, 2)) supp = findall(!=(0), Z[:,c_s] - Z[:,c_p]) @@ -977,3 +976,18 @@ function isconcentrationrobust(rn::ReactionSystem, species::Int) robust_species = robustspecies(rn) return species in robust_species end + +function hasuniqueequilibria(rn::ReactionSystem) + nps = get_networkproperties(rn) + complexes, D = reactioncomplexes(rn) + δ = deficiency(rn) + concordant = isconcordant(rn) + + δ == 0 && return true + satisfiesdeficiencyone(rn) && return true + concordant && return true + !concordant && CNA.ispositivelydependent(rn) && return false + δ == 1 && return deficiencyonealgorithm(rn) + + error("The network is discordant and high deficiency, but it is inconclusive whether the network will have multiple equilibria.") +end From 6e74d76da48dfcdd6b2b7566cb904712d8180463 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 21 Jun 2024 11:33:28 -0400 Subject: [PATCH 02/29] up --- src/network_analysis.jl | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 75cb80b916..7f56f8f464 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -991,3 +991,31 @@ function hasuniqueequilibria(rn::ReactionSystem) error("The network is discordant and high deficiency, but it is inconclusive whether the network will have multiple equilibria.") end + + +### CONCORDANCE + +function isconcordant(rn::ReactionSystem) + +end + +function isstronglyconcordant(rn::ReactionSystem) + +end + +function isdegenerate(rn::ReactionSystem) + +end + +function fullyopenextension(rn::ReactionSystem) + +end + +function speciesreactiongraph(rn::ReactionSystem) + s = numspecies(rn); sm = speciesmap(rn) + r = numreactions(rn); + + G = Digraph(s+r) + adj = zeros(s+r, s+r) + +end From fb04e2e6e3857d3e5aa97f5f0839094e99d5c08f Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 21 Jun 2024 15:29:47 -0400 Subject: [PATCH 03/29] up --- src/network_analysis.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 7f56f8f464..5d6c3a0d2a 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -977,6 +977,8 @@ function isconcentrationrobust(rn::ReactionSystem, species::Int) return species in robust_species end +""" +""" function hasuniqueequilibria(rn::ReactionSystem) nps = get_networkproperties(rn) complexes, D = reactioncomplexes(rn) @@ -989,7 +991,7 @@ function hasuniqueequilibria(rn::ReactionSystem) !concordant && CNA.ispositivelydependent(rn) && return false δ == 1 && return deficiencyonealgorithm(rn) - error("The network is discordant and high deficiency, but it is inconclusive whether the network will have multiple equilibria.") + error("The network is discordant and high deficiency, but this function currently cannot conclude whether the network has the potential to have multiple equilibria.") end From 7992a48bcd36ec88852b711dc1114493dcf80e81 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 28 Jun 2024 13:08:55 -0400 Subject: [PATCH 04/29] init --- src/network_analysis.jl | 87 ++++++++++++++++++++--------------------- 1 file changed, 42 insertions(+), 45 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 5d6c3a0d2a..8d9da86419 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -765,7 +765,48 @@ end """ iscomplexbalanced(rs::ReactionSystem, parametermap) -Constructively compute whether a network will have complex-balanced equilibrium +Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit detailed-balanced equilibrium +solutions, using the Wegscheider conditions, [. A detailed-balanced solution is one for which the rate of every forward reaction exactly equals its reverse reaction. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. +""" + +function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) + if length(parametermap) != numparams(rs) + error("Incorrect number of parameters specified.") + elseif !isreversible(rs) + return false + elseif !all(r -> ismassaction(r, rs), reactions(rs)) + error("The supplied ReactionSystem has reactions that are not ismassaction. Testing for being complex balanced is currently only supported for pure mass action networks.") + end + + pmap = symmap_to_varmap(rs, parametermap) + pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap) + + sm = speciesmap(rs) + cm = reactioncomplexmap(rs) + complexes, D = reactioncomplexes(rs) + rxns = reactions(rs) + nc = length(complexes); nr = numreactions(rs); nm = numspecies(rs) +end + +function isdetailedbalanced(rs::ReactionSystem, parametermap::Vector{Pair{Symbol, Float64}}) + pdict = Dict(parametermap) + isdetailedbalanced(rs, pdict) +end + +function isdetailedbalanced(rs::ReactionSystem, parametermap::Tuple{Pair{Symbol, Float64}}) + pdict = Dict(parametermap) + isdetailedbalanced(rs, pdict) +end + +function isdetailedbalanced(rs::ReactionSystem, parametermap) + error("Parameter map must be a dictionary, tuple, or vector of symbol/value pairs.") +end + + +""" + iscomplexbalanced(rs::ReactionSystem, parametermap) + +Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit complex-balanced equilibrium solutions, following the method in van der Schaft et al., [2015](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. """ @@ -977,47 +1018,3 @@ function isconcentrationrobust(rn::ReactionSystem, species::Int) return species in robust_species end -""" -""" -function hasuniqueequilibria(rn::ReactionSystem) - nps = get_networkproperties(rn) - complexes, D = reactioncomplexes(rn) - δ = deficiency(rn) - concordant = isconcordant(rn) - - δ == 0 && return true - satisfiesdeficiencyone(rn) && return true - concordant && return true - !concordant && CNA.ispositivelydependent(rn) && return false - δ == 1 && return deficiencyonealgorithm(rn) - - error("The network is discordant and high deficiency, but this function currently cannot conclude whether the network has the potential to have multiple equilibria.") -end - - -### CONCORDANCE - -function isconcordant(rn::ReactionSystem) - -end - -function isstronglyconcordant(rn::ReactionSystem) - -end - -function isdegenerate(rn::ReactionSystem) - -end - -function fullyopenextension(rn::ReactionSystem) - -end - -function speciesreactiongraph(rn::ReactionSystem) - s = numspecies(rn); sm = speciesmap(rn) - r = numreactions(rn); - - G = Digraph(s+r) - adj = zeros(s+r, s+r) - -end From aeb00ab6e209f1b54cc9ed9c22a8ac2ae430208a Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 28 Jun 2024 18:05:24 -0400 Subject: [PATCH 05/29] detailed balance impelmentation (tests not run yet) --- src/network_analysis.jl | 56 ++++++++++-- test/network_analysis/network_properties.jl | 94 +++++++++++++++++++++ 2 files changed, 145 insertions(+), 5 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 8d9da86419..6ddc00729f 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -404,6 +404,20 @@ function isterminal(lc::Vector, rn::ReactionSystem) true end +function isforestlike(rn::ReactionSystem) + subnets = subnetworks(rn) + + im = get_networkproperties(rn).incidencemat + isempty(im) && + error("Error, please call reactioncomplexes(rn::ReactionSystem) to ensure the incidence matrix has been cached.") + sparseig = issparse(im) + for subnet in subnets + nps = get_networkproperties(subnet) + isempty(nps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig) + end + all(Graphs.is_tree ∘ Graph ∘ incidencematgraph, subnets) +end + @doc raw""" deficiency(rn::ReactionSystem) @@ -766,7 +780,7 @@ end iscomplexbalanced(rs::ReactionSystem, parametermap) Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit detailed-balanced equilibrium -solutions, using the Wegscheider conditions, [. A detailed-balanced solution is one for which the rate of every forward reaction exactly equals its reverse reaction. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. +solutions, using the Wegscheider conditions, [Feinberg, 1989](https://www.sciencedirect.com/science/article/pii/0009250989851243). A detailed-balanced solution is one for which the rate of every forward reaction exactly equals its reverse reaction. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. """ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) @@ -778,14 +792,46 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) error("The supplied ReactionSystem has reactions that are not ismassaction. Testing for being complex balanced is currently only supported for pure mass action networks.") end + isforestlike(rs) && deficiency(rs) == 0 && return true + pmap = symmap_to_varmap(rs, parametermap) pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap) - sm = speciesmap(rs) - cm = reactioncomplexmap(rs) complexes, D = reactioncomplexes(rs) - rxns = reactions(rs) - nc = length(complexes); nr = numreactions(rs); nm = numspecies(rs) + img = incidencematgraph(rn) + undir_img = Graph(incidencematgraph(rn)) + K = ratematrix(rn, pmap) + + spanning_forest = kruskal_mst(uimg) + undir_spanning_forest = edges(Graph(spanning_forest)) + outofforest_rxns = findall(∉(undir_spanning_forest), edges(undir_img)) + + # Independent Cycle Conditions: for any cycle we create by adding in an out-of-forest reaction, the product of forward reaction rates over the cycle must equal the product of reverse reaction rates over the cycle. + for rxn in outofforest_rxns + g = Graph([undir_spanning_forest..., collect(edges(g))[rxn]]) + ic = cycle_basis(g) + fwd = prod([K[ic[r], ic[r+1]] for r in 1:length(ic)-1]) * K[ic[end], ic[1]] + rev = prod([K[ic[r+1], ic[r]] for r in 1:length(ic)-1]) * K[ic[1], ic[end]] + fwd == rev ? continue : return false + end + + # Spanning Forest Conditions: for non-deficiency 0 networks, we get an additional δ equations. Choose an orientation for each reaction pair in the spanning forest (by default, we will go from smaller to larger indices). + + if deficiency(rn) > 0 + rxn_idxs = findall(∈(spanning_forest), edges(img)) + S_F = netstoichmat(rn)[:, fwdrxn_idxs] + sols = nullspace(S_F) + @assert size(sols, 2) == deficiency(rn) + + for i in 1:size(sols, 2) + α = sols[:, i] + fwd = prod([K[src(e), dst(e)]^α[i] for (e,i) in zip(spanning_forest,1:length(α))]) + rev = prod([K[dst(e), src(e)]^α[i] for (e,i) in zip(spanning_forest, 1:length(α))]) + fwd == rev ? continue : return false + end + end + + true end function isdetailedbalanced(rs::ReactionSystem, parametermap::Vector{Pair{Symbol, Float64}}) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 0f28a14ada..6dfd070f65 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -295,6 +295,7 @@ let rates = Dict(zip(parameters(rn), k)) @test Catalyst.iscomplexbalanced(rn, rates) == true end + let rn = @reaction_network begin (k2, k1), A + E <--> AE @@ -324,6 +325,7 @@ let k = rand(rng, numparams(rn)) rates = Dict(zip(parameters(rn), k)) @test Catalyst.iscomplexbalanced(rn, rates) == true + @test Catalyst.isdetailedbalanced(rn, rates) == false end ### STRONG LINKAGE CLASS TESTS @@ -430,3 +432,95 @@ let @test Catalyst.robustspecies(EnvZ_OmpR) == [6] end + + +### Complex and detailed balance tests + +# The following network is conditionally complex balanced - it only + + +# Reversible, forest-like deficiency zero network - should be detailed balance for any choice of rate constants. +let + rn = @reaction_network begin + (k1, k2), A <--> B + C + (k3, k4), A <--> D + (k5, k6), A + D <--> E + (k7, k8), A + D <--> G + (k9, k10), G <--> 2F + (k11, k12), A + E <--> H + end + + k1 = rand(rng, numparams(rn)) + rates1 = Dict(zip(parameters(rn), k1)) + k2 = rand(StableRNG(232), numparams(rn)) + rates2 = Dict(zip(parameters(rn), k2)) + + @test isdetailedbalanced(rn, rates1) == true + @test isdetailedbalanced(rn, rates2) == true +end + + +# Simple connected reversible network +let + rn = @reaction_network begin + (k1, k2), A <--> B + (k3, k4), B <--> C + (k5, k6), C <--> A + end + + rates1 = [k1=>1.0, k2=>1.0, k3=>1.0, k4=>1.0, k5=>1.0, k6=>1.0] + @test isdetailedbalanced(rn, rates1) == true + rates2 = [k1=>2.0, k2=>1.0, k3=>1.0, k4=>1.0, k5=>1.0, k6=>1.0] + @test isdetailedbalanced(rn, rates2) == false +end + +# Independent cycle tests: the following reaction entwork has 3 out-of-forest reactions. +let + rn = @reaction_network begin + (k1, k2), A <--> B + C + (k3, k4), A <--> D + (k5, k6), B + C <--> D + (k7, k8), A + D <--> E + (k9, k10), G <--> 2F + (k11, k12), A + D <--> G + (k13, k14), G <--> E + (k15, k16), 2F <--> E + (k17, k18), A + E <--> H + end + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test isdetailedbalanced(rn, rates) == false + + rates[k6] = rates[k1]*rates[k4]*rates[k5] / (rates[k2]*rates[k3]) + rates[k14] = rates[k13]*rates[k11]*rates[k8] / (rates[k12]*rates[k7]) + rates[k12] = rates[k8]*rates[k15]*rates[k9]*rates[k11] / (rates[k7]*rates[k16]*rates[k10]) + @test isdetailedbalanced(rn, rates) == true +end + +# Deficiency two network: the following reaction network must satisfy both the independent cycle conditions and the spanning forest conditions +let + rn = @reaction_network begin + (k1, k2), 3A <--> A + 2B + (k3, k4), A + 2B <--> 3B + (k5, k6), 3B <--> 2A + B + (k7, k8), 2A + B <--> 3A + (k9, k10), 3A <--> 3B + end + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test isdetailedbalanced(rn, rates) == false + + rates[k8] = rates[k7]*rates[k5]*rates[k9] / (rates[k6]*rates[k10]) + rates[k3] = rates[k2]*rates[k4]*rates[k9] / (rates[k1]*rates[k10]) + @test isdetailedbalanced(rn, rates) == false + + cons = rates[k6] / rates[k5] + rates[k1] = rates[k2] * cons + rates[k9] = rates[k10] * cons^(3/2) + rates[k8] = rates[k7]*rates[k5]*rates[k9] / (rates[k6]*rates[k10]) + rates[k3] = rates[k2]*rates[k4]*rates[k9] / (rates[k1]*rates[k10]) + @test isdeatiledbalanced(rn, rates) == true + @test isdetailedbalanced(rn, rates) == false +end From ba1f25086116d337975b38b57cb14dcfc63b3434 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 2 Jul 2024 20:23:08 -0400 Subject: [PATCH 06/29] tests passing --- src/network_analysis.jl | 45 +++++++++------ test/network_analysis/network_properties.jl | 64 ++++++++++++--------- 2 files changed, 63 insertions(+), 46 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 6ddc00729f..df3a9106ef 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -415,7 +415,7 @@ function isforestlike(rn::ReactionSystem) nps = get_networkproperties(subnet) isempty(nps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig) end - all(Graphs.is_tree ∘ Graph ∘ incidencematgraph, subnets) + all(Graphs.is_tree ∘ SimpleGraph ∘ incidencematgraph, subnets) end @doc raw""" @@ -797,43 +797,50 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) pmap = symmap_to_varmap(rs, parametermap) pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap) + # Construct reaction-complex graph complexes, D = reactioncomplexes(rs) - img = incidencematgraph(rn) - undir_img = Graph(incidencematgraph(rn)) - K = ratematrix(rn, pmap) + img = incidencematgraph(rs) + undir_img = SimpleGraph(incidencematgraph(rs)) + K = ratematrix(rs, pmap) - spanning_forest = kruskal_mst(uimg) - undir_spanning_forest = edges(Graph(spanning_forest)) - outofforest_rxns = findall(∉(undir_spanning_forest), edges(undir_img)) + spanning_forest = Graphs.kruskal_mst(undir_img) + outofforest_edges = setdiff(collect(edges(undir_img)), spanning_forest) # Independent Cycle Conditions: for any cycle we create by adding in an out-of-forest reaction, the product of forward reaction rates over the cycle must equal the product of reverse reaction rates over the cycle. - for rxn in outofforest_rxns - g = Graph([undir_spanning_forest..., collect(edges(g))[rxn]]) - ic = cycle_basis(g) + for edge in outofforest_edges + g = SimpleGraph([spanning_forest..., edge]) + ic = Graphs.cycle_basis(g)[1] fwd = prod([K[ic[r], ic[r+1]] for r in 1:length(ic)-1]) * K[ic[end], ic[1]] rev = prod([K[ic[r+1], ic[r]] for r in 1:length(ic)-1]) * K[ic[1], ic[end]] - fwd == rev ? continue : return false + fwd ≈ rev ? continue : return false end - # Spanning Forest Conditions: for non-deficiency 0 networks, we get an additional δ equations. Choose an orientation for each reaction pair in the spanning forest (by default, we will go from smaller to larger indices). + # Spanning Forest Conditions: for non-deficiency 0 networks, we get an additional δ equations. Choose an orientation for each reaction pair in the spanning forest (we will take the one given by default from kruskal_mst). - if deficiency(rn) > 0 - rxn_idxs = findall(∈(spanning_forest), edges(img)) - S_F = netstoichmat(rn)[:, fwdrxn_idxs] + if deficiency(rs) > 0 + rxn_idxs = [edgeindex(D, Graphs.src(e), Graphs.dst(e)) for e in spanning_forest] + S_F = netstoichmat(rs)[:, rxn_idxs] sols = nullspace(S_F) - @assert size(sols, 2) == deficiency(rn) + @assert size(sols, 2) == deficiency(rs) for i in 1:size(sols, 2) α = sols[:, i] - fwd = prod([K[src(e), dst(e)]^α[i] for (e,i) in zip(spanning_forest,1:length(α))]) - rev = prod([K[dst(e), src(e)]^α[i] for (e,i) in zip(spanning_forest, 1:length(α))]) - fwd == rev ? continue : return false + fwd = prod([K[Graphs.src(e), Graphs.dst(e)]^α[i] for (e,i) in zip(spanning_forest,1:length(α))]) + rev = prod([K[Graphs.dst(e), Graphs.src(e)]^α[i] for (e,i) in zip(spanning_forest, 1:length(α))]) + fwd ≈ rev ? continue : return false end end true end +function edgeindex(imat::Matrix{Int64}, src::Int64, dst::Int64) + for i in 1:size(imat, 2) + (imat[src, i] == -1) && (imat[dst, i] == 1) && return i + end + error("This edge does not exist in this reaction graph.") +end + function isdetailedbalanced(rs::ReactionSystem, parametermap::Vector{Pair{Symbol, Float64}}) pdict = Dict(parametermap) isdetailedbalanced(rs, pdict) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 6dfd070f65..35a24199e2 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -325,7 +325,7 @@ let k = rand(rng, numparams(rn)) rates = Dict(zip(parameters(rn), k)) @test Catalyst.iscomplexbalanced(rn, rates) == true - @test Catalyst.isdetailedbalanced(rn, rates) == false + @test Catalyst.Catalyst.isdetailedbalanced(rn, rates) == false end ### STRONG LINKAGE CLASS TESTS @@ -438,7 +438,6 @@ end # The following network is conditionally complex balanced - it only - # Reversible, forest-like deficiency zero network - should be detailed balance for any choice of rate constants. let rn = @reaction_network begin @@ -455,11 +454,12 @@ let k2 = rand(StableRNG(232), numparams(rn)) rates2 = Dict(zip(parameters(rn), k2)) - @test isdetailedbalanced(rn, rates1) == true - @test isdetailedbalanced(rn, rates2) == true + rcs, D = reactioncomplexes(rn) + @test Catalyst.isforestlike(rn) == true + @test Catalyst.isdetailedbalanced(rn, rates1) == true + @test Catalyst.isdetailedbalanced(rn, rates2) == true end - # Simple connected reversible network let rn = @reaction_network begin @@ -468,10 +468,11 @@ let (k5, k6), C <--> A end - rates1 = [k1=>1.0, k2=>1.0, k3=>1.0, k4=>1.0, k5=>1.0, k6=>1.0] - @test isdetailedbalanced(rn, rates1) == true - rates2 = [k1=>2.0, k2=>1.0, k3=>1.0, k4=>1.0, k5=>1.0, k6=>1.0] - @test isdetailedbalanced(rn, rates2) == false + rcs, D = reactioncomplexes(rn) + rates1 = [:k1=>1.0, :k2=>1.0, :k3=>1.0, :k4=>1.0, :k5=>1.0, :k6=>1.0] + @test Catalyst.isdetailedbalanced(rn, rates1) == true + rates2 = [:k1=>2.0, :k2=>1.0, :k3=>1.0, :k4=>1.0, :k5=>1.0, :k6=>1.0] + @test Catalyst.isdetailedbalanced(rn, rates2) == false end # Independent cycle tests: the following reaction entwork has 3 out-of-forest reactions. @@ -488,14 +489,17 @@ let (k17, k18), A + E <--> H end + rcs, D = reactioncomplexes(rn) k = rand(rng, numparams(rn)) + p = parameters(rn) rates = Dict(zip(parameters(rn), k)) - @test isdetailedbalanced(rn, rates) == false + @test Catalyst.isdetailedbalanced(rn, rates) == false - rates[k6] = rates[k1]*rates[k4]*rates[k5] / (rates[k2]*rates[k3]) - rates[k14] = rates[k13]*rates[k11]*rates[k8] / (rates[k12]*rates[k7]) - rates[k12] = rates[k8]*rates[k15]*rates[k9]*rates[k11] / (rates[k7]*rates[k16]*rates[k10]) - @test isdetailedbalanced(rn, rates) == true + # Adjust rate constants to obey the independent cycle conditions. + rates[p[6]] = rates[p[1]]*rates[p[4]]*rates[p[5]] / (rates[p[2]]*rates[p[3]]) + rates[p[14]] = rates[p[13]]*rates[p[11]]*rates[p[8]] / (rates[p[12]]*rates[p[7]]) + rates[p[16]] = rates[p[8]]*rates[p[15]]*rates[p[9]]*rates[p[11]] / (rates[p[7]]*rates[p[12]]*rates[p[10]]) + @test Catalyst.isdetailedbalanced(rn, rates) == true end # Deficiency two network: the following reaction network must satisfy both the independent cycle conditions and the spanning forest conditions @@ -508,19 +512,25 @@ let (k9, k10), 3A <--> 3B end + rcs, D = reactioncomplexes(rn) + @test Catalyst.edgeindex(D, 1, 2) == 1 + @test Catalyst.edgeindex(D, 4, 3) == 6 k = rand(rng, numparams(rn)) + p = parameters(rn) rates = Dict(zip(parameters(rn), k)) - @test isdetailedbalanced(rn, rates) == false - - rates[k8] = rates[k7]*rates[k5]*rates[k9] / (rates[k6]*rates[k10]) - rates[k3] = rates[k2]*rates[k4]*rates[k9] / (rates[k1]*rates[k10]) - @test isdetailedbalanced(rn, rates) == false - - cons = rates[k6] / rates[k5] - rates[k1] = rates[k2] * cons - rates[k9] = rates[k10] * cons^(3/2) - rates[k8] = rates[k7]*rates[k5]*rates[k9] / (rates[k6]*rates[k10]) - rates[k3] = rates[k2]*rates[k4]*rates[k9] / (rates[k1]*rates[k10]) - @test isdeatiledbalanced(rn, rates) == true - @test isdetailedbalanced(rn, rates) == false + @test Catalyst.isdetailedbalanced(rn, rates) == false + + # Adjust rate constants to fulfill independent cycle conditions. + rates[p[8]] = rates[p[7]]*rates[p[5]]*rates[p[9]] / (rates[p[6]]*rates[p[10]]) + rates[p[3]] = rates[p[2]]*rates[p[4]]*rates[p[9]] / (rates[p[1]]*rates[p[10]]) + @test Catalyst.isdetailedbalanced(rn, rates) == false + # Should still fail - doesn't satisfy spanning forest conditions. + + # Adjust rate constants to fulfill spanning forest conditions. + cons = rates[p[6]] / rates[p[5]] + rates[p[1]] = rates[p[2]] * cons + rates[p[9]] = rates[p[10]] * cons^(3/2) + rates[p[8]] = rates[p[7]]*rates[p[5]]*rates[p[9]] / (rates[p[6]]*rates[p[10]]) + rates[p[3]] = rates[p[2]]*rates[p[4]]*rates[p[9]] / (rates[p[1]]*rates[p[10]]) + @test Catalyst.isdetailedbalanced(rn, rates) == true end From f308aea55164c5a490655b8345421984da7a45db Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 2 Jul 2024 20:23:38 -0400 Subject: [PATCH 07/29] formatter --- src/Catalyst.jl | 3 ++- src/network_analysis.jl | 49 +++++++++++++++++++++-------------------- 2 files changed, 27 insertions(+), 25 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 df3a9106ef..374db93ce4 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) @@ -404,7 +404,7 @@ function isterminal(lc::Vector, rn::ReactionSystem) true end -function isforestlike(rn::ReactionSystem) +function isforestlike(rn::ReactionSystem) subnets = subnetworks(rn) im = get_networkproperties(rn).incidencemat @@ -783,7 +783,7 @@ Constructively compute whether a kinetic system (a reaction network with a set o solutions, using the Wegscheider conditions, [Feinberg, 1989](https://www.sciencedirect.com/science/article/pii/0009250989851243). A detailed-balanced solution is one for which the rate of every forward reaction exactly equals its reverse reaction. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. """ -function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) +function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) if length(parametermap) != numparams(rs) error("Incorrect number of parameters specified.") elseif !isreversible(rs) @@ -796,7 +796,7 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) pmap = symmap_to_varmap(rs, parametermap) pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap) - + # Construct reaction-complex graph complexes, D = reactioncomplexes(rs) img = incidencematgraph(rs) @@ -810,23 +810,25 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) for edge in outofforest_edges g = SimpleGraph([spanning_forest..., edge]) ic = Graphs.cycle_basis(g)[1] - fwd = prod([K[ic[r], ic[r+1]] for r in 1:length(ic)-1]) * K[ic[end], ic[1]] - rev = prod([K[ic[r+1], ic[r]] for r in 1:length(ic)-1]) * K[ic[1], ic[end]] + fwd = prod([K[ic[r], ic[r + 1]] for r in 1:(length(ic) - 1)]) * K[ic[end], ic[1]] + rev = prod([K[ic[r + 1], ic[r]] for r in 1:(length(ic) - 1)]) * K[ic[1], ic[end]] fwd ≈ rev ? continue : return false end - + # Spanning Forest Conditions: for non-deficiency 0 networks, we get an additional δ equations. Choose an orientation for each reaction pair in the spanning forest (we will take the one given by default from kruskal_mst). - + if deficiency(rs) > 0 rxn_idxs = [edgeindex(D, Graphs.src(e), Graphs.dst(e)) for e in spanning_forest] S_F = netstoichmat(rs)[:, rxn_idxs] - sols = nullspace(S_F) + sols = nullspace(S_F) @assert size(sols, 2) == deficiency(rs) for i in 1:size(sols, 2) α = sols[:, i] - fwd = prod([K[Graphs.src(e), Graphs.dst(e)]^α[i] for (e,i) in zip(spanning_forest,1:length(α))]) - rev = prod([K[Graphs.dst(e), Graphs.src(e)]^α[i] for (e,i) in zip(spanning_forest, 1:length(α))]) + fwd = prod([K[Graphs.src(e), Graphs.dst(e)]^α[i] + for (e, i) in zip(spanning_forest, 1:length(α))]) + rev = prod([K[Graphs.dst(e), Graphs.src(e)]^α[i] + for (e, i) in zip(spanning_forest, 1:length(α))]) fwd ≈ rev ? continue : return false end end @@ -834,7 +836,7 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) true end -function edgeindex(imat::Matrix{Int64}, src::Int64, dst::Int64) +function edgeindex(imat::Matrix{Int64}, src::Int64, dst::Int64) for i in 1:size(imat, 2) (imat[src, i] == -1) && (imat[dst, i] == 1) && return i end @@ -855,7 +857,6 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap) error("Parameter map must be a dictionary, tuple, or vector of symbol/value pairs.") end - """ iscomplexbalanced(rs::ReactionSystem, parametermap) @@ -1034,12 +1035,13 @@ 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 @@ -1047,27 +1049,26 @@ end #function deficiencyonealgorithm(rn::ReactionSystem) #end -function robustspecies(rn::ReactionSystem) +function robustspecies(rn::ReactionSystem) complexes, D = reactioncomplexes(rn) robust_species = Int64[] - + if deficiency(rn) != 1 error("This method 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...)) 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 bec6883a13421ed759700681eb400ac038b814ad Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 16 Jul 2024 12:41:47 -0400 Subject: [PATCH 08/29] up --- src/network_analysis.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 1a43813ca3..0e1bddb145 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -410,11 +410,10 @@ end function isforestlike(rn::ReactionSystem) subnets = subnetworks(rn) + nps = get_networkproperties(rn) - im = get_networkproperties(rn).incidencemat - isempty(im) && - error("Error, please call reactioncomplexes(rn::ReactionSystem) to ensure the incidence matrix has been cached.") - sparseig = issparse(im) + isempty(nps.incidencemat) && reactioncomplexes(rn) + sparseig = issparse(incidencemat) for subnet in subnets nps = get_networkproperties(subnet) isempty(nps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig) From d77cb1c4dcd539d8e8bc83df21faf705630d8247 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 16 Jul 2024 14:41:29 -0400 Subject: [PATCH 09/29] 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 867c68e8ed8ea99a7eb1eebe632e27285829c41c Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 17 Jul 2024 14:35:14 -0400 Subject: [PATCH 10/29] merge fix --- test/network_analysis/network_properties.jl | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 37eb5068f4..13a17ccac6 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -432,7 +432,6 @@ let @test issubset([[3,4], [5,6,7]], tslcs) end -<<<<<<< HEAD ### 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. @@ -515,12 +514,6 @@ let @test Catalyst.satisfiesdeficiencyone(rn) == true end - - - - - - ### Complex and detailed balance tests # The following network is conditionally complex balanced - it only @@ -620,7 +613,8 @@ let rates[p[8]] = rates[p[7]]*rates[p[5]]*rates[p[9]] / (rates[p[6]]*rates[p[10]]) rates[p[3]] = rates[p[2]]*rates[p[4]]*rates[p[9]] / (rates[p[1]]*rates[p[10]]) @test Catalyst.isdetailedbalanced(rn, rates) == true -======= +end + ### Other Network Properties Tests ### # Tests outgoing complexes matrices (1). @@ -730,5 +724,4 @@ let Catalyst.ratematrix(rn, rates_tup) == rate_mat Catalyst.ratematrix(rn, rates_dict) == rate_mat @test_throws Exception Catalyst.iscomplexbalanced(rn, rates_invalid) ->>>>>>> origin/master end From 2a2ad06d251466de7bd8cc21bb4045bf26bd66c7 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 17 Jul 2024 14:36:23 -0400 Subject: [PATCH 11/29] merge fix --- src/network_analysis.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 9c309985a0..75d2747df5 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -836,13 +836,9 @@ end """ iscomplexbalanced(rs::ReactionSystem, parametermap) -Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit complex-balanced equilibrium -solutions, following the method in van der Schaft et al., [2015](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. -======= Constructively compute whether a network will have complex-balanced equilibrium solutions, following the method in van der Schaft et al., [2015](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. ->>>>>>> origin/master """ function iscomplexbalanced(rs::ReactionSystem, parametermap::Dict) From 6fbd76882e89fd95f8187fc72a1c825c27fb29a4 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 17 Jul 2024 14:37:27 -0400 Subject: [PATCH 12/29] merge fix --- test/network_analysis/network_properties.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 13a17ccac6..02b1f21504 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -345,7 +345,7 @@ let k = rand(rng, numparams(rn)) rates = Dict(zip(parameters(rn), k)) @test Catalyst.iscomplexbalanced(rn, rates) == true - @test Catalyst.Catalyst.isdetailedbalanced(rn, rates) == false + @test Catalyst.isdetailedbalanced(rn, rates) == false end ### STRONG LINKAGE CLASS TESTS From d580647e993b40a99e02470b2fc9244b0c5667b9 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 17 Jul 2024 14:38:36 -0400 Subject: [PATCH 13/29] fmratter --- 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 fa0df9dc4ac990f5c35eb722260740b88b3ef046 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 17 Jul 2024 14:57:39 -0400 Subject: [PATCH 14/29] 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 75d2747df5..5187a89abe 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -403,7 +403,7 @@ function isforestlike(rn::ReactionSystem) nps = get_networkproperties(rn) isempty(nps.incidencemat) && reactioncomplexes(rn) - sparseig = issparse(incidencemat) + sparseig = issparse(nps.incidencemat) for subnet in subnets nps = get_networkproperties(subnet) isempty(nps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig) From 10bcf07bee354b581917bc14602a52e28495ae76 Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Mon, 4 Nov 2024 00:24:43 +0000 Subject: [PATCH 15/29] CompatHelper: bump compat for NonlinearSolve to 4 for package docs, (keep existing compat) --- docs/Project.toml | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 39cdb81ad0..08c3ed624d 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -35,7 +35,6 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" -#StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] @@ -58,7 +57,7 @@ JumpProcesses = "9.13.2" Latexify = "0.16.5" LinearSolve = "2.30" ModelingToolkit = "9.32" -NonlinearSolve = "3.12" +NonlinearSolve = "3.12, 4" Optim = "1.9" Optimization = "4" OptimizationBBO = "0.4" @@ -74,5 +73,4 @@ SpecialFunctions = "2.4" StaticArrays = "1.9" SteadyStateDiffEq = "2.2" StochasticDiffEq = "6.65" -#StructuralIdentifiability = "0.5.8" Symbolics = "5.30.1, 6" From 93f06efdc3b0515d3bdd7e7d3c71757e93545c0e Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 19 Jun 2024 12:01:42 -0400 Subject: [PATCH 16/29] Mulitple equilibria --- src/network_analysis.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 3e42e064c1..8edb15bab1 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -1027,3 +1027,18 @@ function robustspecies(rn::ReactionSystem) nps.robustspecies end + +function hasuniqueequilibria(rn::ReactionSystem) + nps = get_networkproperties(rn) + complexes, D = reactioncomplexes(rn) + δ = deficiency(rn) + concordant = isconcordant(rn) + + δ == 0 && return true + satisfiesdeficiencyone(rn) && return true + concordant && return true + !concordant && CNA.ispositivelydependent(rn) && return false + δ == 1 && return deficiencyonealgorithm(rn) + + error("The network is discordant and high deficiency, but it is inconclusive whether the network will have multiple equilibria.") +end From db84110207d3677e97a57c3f3745efe865342f37 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 21 Jun 2024 11:33:28 -0400 Subject: [PATCH 17/29] up --- src/network_analysis.jl | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 8edb15bab1..841455ef72 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -1042,3 +1042,31 @@ function hasuniqueequilibria(rn::ReactionSystem) error("The network is discordant and high deficiency, but it is inconclusive whether the network will have multiple equilibria.") end + + +### CONCORDANCE + +function isconcordant(rn::ReactionSystem) + +end + +function isstronglyconcordant(rn::ReactionSystem) + +end + +function isdegenerate(rn::ReactionSystem) + +end + +function fullyopenextension(rn::ReactionSystem) + +end + +function speciesreactiongraph(rn::ReactionSystem) + s = numspecies(rn); sm = speciesmap(rn) + r = numreactions(rn); + + G = Digraph(s+r) + adj = zeros(s+r, s+r) + +end From 9ad2e0b27f50852bcdd0d2eb5e711b124d45529f Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 21 Jun 2024 15:29:47 -0400 Subject: [PATCH 18/29] up --- src/network_analysis.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 841455ef72..eb03e48964 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -1028,6 +1028,8 @@ function robustspecies(rn::ReactionSystem) nps.robustspecies end +""" +""" function hasuniqueequilibria(rn::ReactionSystem) nps = get_networkproperties(rn) complexes, D = reactioncomplexes(rn) @@ -1040,7 +1042,7 @@ function hasuniqueequilibria(rn::ReactionSystem) !concordant && CNA.ispositivelydependent(rn) && return false δ == 1 && return deficiencyonealgorithm(rn) - error("The network is discordant and high deficiency, but it is inconclusive whether the network will have multiple equilibria.") + error("The network is discordant and high deficiency, but this function currently cannot conclude whether the network has the potential to have multiple equilibria.") end From bd7402ff99fcc457d3664f8c7279e400b53e2fb9 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 28 Jun 2024 13:08:55 -0400 Subject: [PATCH 19/29] init --- src/network_analysis.jl | 44 ----------------------------------------- 1 file changed, 44 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index eb03e48964..278217cabc 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -1028,47 +1028,3 @@ function robustspecies(rn::ReactionSystem) nps.robustspecies end -""" -""" -function hasuniqueequilibria(rn::ReactionSystem) - nps = get_networkproperties(rn) - complexes, D = reactioncomplexes(rn) - δ = deficiency(rn) - concordant = isconcordant(rn) - - δ == 0 && return true - satisfiesdeficiencyone(rn) && return true - concordant && return true - !concordant && CNA.ispositivelydependent(rn) && return false - δ == 1 && return deficiencyonealgorithm(rn) - - error("The network is discordant and high deficiency, but this function currently cannot conclude whether the network has the potential to have multiple equilibria.") -end - - -### CONCORDANCE - -function isconcordant(rn::ReactionSystem) - -end - -function isstronglyconcordant(rn::ReactionSystem) - -end - -function isdegenerate(rn::ReactionSystem) - -end - -function fullyopenextension(rn::ReactionSystem) - -end - -function speciesreactiongraph(rn::ReactionSystem) - s = numspecies(rn); sm = speciesmap(rn) - r = numreactions(rn); - - G = Digraph(s+r) - adj = zeros(s+r, s+r) - -end From fe80200c9131ce7087e5f644539ebfbb0040a334 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 28 Jun 2024 18:05:24 -0400 Subject: [PATCH 20/29] detailed balance impelmentation (tests not run yet) --- src/network_analysis.jl | 92 +++++++++++++++++++- test/network_analysis/network_properties.jl | 93 +++++++++++++++++++++ 2 files changed, 182 insertions(+), 3 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 278217cabc..bfd34d2b34 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -398,6 +398,20 @@ function isterminal(lc::Vector, rn::ReactionSystem) true end +function isforestlike(rn::ReactionSystem) + subnets = subnetworks(rn) + + im = get_networkproperties(rn).incidencemat + isempty(im) && + error("Error, please call reactioncomplexes(rn::ReactionSystem) to ensure the incidence matrix has been cached.") + sparseig = issparse(im) + for subnet in subnets + nps = get_networkproperties(subnet) + isempty(nps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig) + end + all(Graphs.is_tree ∘ Graph ∘ incidencematgraph, subnets) +end + @doc raw""" deficiency(rn::ReactionSystem) @@ -727,9 +741,81 @@ end """ iscomplexbalanced(rs::ReactionSystem, parametermap) -Constructively compute whether a network will have complex-balanced equilibrium -solutions, following the method in van der Schaft et al., [2015](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). -Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. +Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit detailed-balanced equilibrium +solutions, using the Wegscheider conditions, [Feinberg, 1989](https://www.sciencedirect.com/science/article/pii/0009250989851243). A detailed-balanced solution is one for which the rate of every forward reaction exactly equals its reverse reaction. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. +""" + +function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) + if length(parametermap) != numparams(rs) + error("Incorrect number of parameters specified.") + elseif !isreversible(rs) + return false + elseif !all(r -> ismassaction(r, rs), reactions(rs)) + error("The supplied ReactionSystem has reactions that are not ismassaction. Testing for being complex balanced is currently only supported for pure mass action networks.") + end + + isforestlike(rs) && deficiency(rs) == 0 && return true + + pmap = symmap_to_varmap(rs, parametermap) + pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap) + + complexes, D = reactioncomplexes(rs) + img = incidencematgraph(rn) + undir_img = Graph(incidencematgraph(rn)) + K = ratematrix(rn, pmap) + + spanning_forest = kruskal_mst(uimg) + undir_spanning_forest = edges(Graph(spanning_forest)) + outofforest_rxns = findall(∉(undir_spanning_forest), edges(undir_img)) + + # Independent Cycle Conditions: for any cycle we create by adding in an out-of-forest reaction, the product of forward reaction rates over the cycle must equal the product of reverse reaction rates over the cycle. + for rxn in outofforest_rxns + g = Graph([undir_spanning_forest..., collect(edges(g))[rxn]]) + ic = cycle_basis(g) + fwd = prod([K[ic[r], ic[r+1]] for r in 1:length(ic)-1]) * K[ic[end], ic[1]] + rev = prod([K[ic[r+1], ic[r]] for r in 1:length(ic)-1]) * K[ic[1], ic[end]] + fwd == rev ? continue : return false + end + + # Spanning Forest Conditions: for non-deficiency 0 networks, we get an additional δ equations. Choose an orientation for each reaction pair in the spanning forest (by default, we will go from smaller to larger indices). + + if deficiency(rn) > 0 + rxn_idxs = findall(∈(spanning_forest), edges(img)) + S_F = netstoichmat(rn)[:, fwdrxn_idxs] + sols = nullspace(S_F) + @assert size(sols, 2) == deficiency(rn) + + for i in 1:size(sols, 2) + α = sols[:, i] + fwd = prod([K[src(e), dst(e)]^α[i] for (e,i) in zip(spanning_forest,1:length(α))]) + rev = prod([K[dst(e), src(e)]^α[i] for (e,i) in zip(spanning_forest, 1:length(α))]) + fwd == rev ? continue : return false + end + end + + true +end + +function isdetailedbalanced(rs::ReactionSystem, parametermap::Vector{Pair{Symbol, Float64}}) + pdict = Dict(parametermap) + isdetailedbalanced(rs, pdict) +end + +function isdetailedbalanced(rs::ReactionSystem, parametermap::Tuple{Pair{Symbol, Float64}}) + pdict = Dict(parametermap) + isdetailedbalanced(rs, pdict) +end + +function isdetailedbalanced(rs::ReactionSystem, parametermap) + error("Parameter map must be a dictionary, tuple, or vector of symbol/value pairs.") +end + + +""" + iscomplexbalanced(rs::ReactionSystem, parametermap) + +Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit complex-balanced equilibrium +solutions, following the method in van der Schaft et al., [2015](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. """ function iscomplexbalanced(rs::ReactionSystem, parametermap::Dict) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index decdbaa0ac..f700bee70b 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -355,6 +355,7 @@ let k = rand(rng, numparams(rn)) rates = Dict(zip(parameters(rn), k)) @test Catalyst.iscomplexbalanced(rn, rates) == true + @test Catalyst.isdetailedbalanced(rn, rates) == false end ### STRONG LINKAGE CLASS TESTS @@ -637,6 +638,98 @@ let @test Catalyst.robustspecies(EnvZ_OmpR) == [6] end + +### Complex and detailed balance tests + +# The following network is conditionally complex balanced - it only + + +# Reversible, forest-like deficiency zero network - should be detailed balance for any choice of rate constants. +let + rn = @reaction_network begin + (k1, k2), A <--> B + C + (k3, k4), A <--> D + (k5, k6), A + D <--> E + (k7, k8), A + D <--> G + (k9, k10), G <--> 2F + (k11, k12), A + E <--> H + end + + k1 = rand(rng, numparams(rn)) + rates1 = Dict(zip(parameters(rn), k1)) + k2 = rand(StableRNG(232), numparams(rn)) + rates2 = Dict(zip(parameters(rn), k2)) + + @test isdetailedbalanced(rn, rates1) == true + @test isdetailedbalanced(rn, rates2) == true +end + + +# Simple connected reversible network +let + rn = @reaction_network begin + (k1, k2), A <--> B + (k3, k4), B <--> C + (k5, k6), C <--> A + end + + rates1 = [k1=>1.0, k2=>1.0, k3=>1.0, k4=>1.0, k5=>1.0, k6=>1.0] + @test isdetailedbalanced(rn, rates1) == true + rates2 = [k1=>2.0, k2=>1.0, k3=>1.0, k4=>1.0, k5=>1.0, k6=>1.0] + @test isdetailedbalanced(rn, rates2) == false +end + +# Independent cycle tests: the following reaction entwork has 3 out-of-forest reactions. +let + rn = @reaction_network begin + (k1, k2), A <--> B + C + (k3, k4), A <--> D + (k5, k6), B + C <--> D + (k7, k8), A + D <--> E + (k9, k10), G <--> 2F + (k11, k12), A + D <--> G + (k13, k14), G <--> E + (k15, k16), 2F <--> E + (k17, k18), A + E <--> H + end + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test isdetailedbalanced(rn, rates) == false + + rates[k6] = rates[k1]*rates[k4]*rates[k5] / (rates[k2]*rates[k3]) + rates[k14] = rates[k13]*rates[k11]*rates[k8] / (rates[k12]*rates[k7]) + rates[k12] = rates[k8]*rates[k15]*rates[k9]*rates[k11] / (rates[k7]*rates[k16]*rates[k10]) + @test isdetailedbalanced(rn, rates) == true +end + +# Deficiency two network: the following reaction network must satisfy both the independent cycle conditions and the spanning forest conditions +let + rn = @reaction_network begin + (k1, k2), 3A <--> A + 2B + (k3, k4), A + 2B <--> 3B + (k5, k6), 3B <--> 2A + B + (k7, k8), 2A + B <--> 3A + (k9, k10), 3A <--> 3B + end + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test isdetailedbalanced(rn, rates) == false + + rates[k8] = rates[k7]*rates[k5]*rates[k9] / (rates[k6]*rates[k10]) + rates[k3] = rates[k2]*rates[k4]*rates[k9] / (rates[k1]*rates[k10]) + @test isdetailedbalanced(rn, rates) == false + + cons = rates[k6] / rates[k5] + rates[k1] = rates[k2] * cons + rates[k9] = rates[k10] * cons^(3/2) + rates[k8] = rates[k7]*rates[k5]*rates[k9] / (rates[k6]*rates[k10]) + rates[k3] = rates[k2]*rates[k4]*rates[k9] / (rates[k1]*rates[k10]) + @test isdeatiledbalanced(rn, rates) == true + @test isdetailedbalanced(rn, rates) == false +end + ### DEFICIENCY ONE TESTS # Fails because there are two terminal linkage classes in the linkage class From 364ee21d28dbe2025b14f8b40d36592583ed5fa9 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 2 Jul 2024 20:23:08 -0400 Subject: [PATCH 21/29] tests passing --- src/network_analysis.jl | 45 +++++++++------ test/network_analysis/network_properties.jl | 64 ++++++++++++--------- 2 files changed, 63 insertions(+), 46 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index bfd34d2b34..9d6d47605e 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -409,7 +409,7 @@ function isforestlike(rn::ReactionSystem) nps = get_networkproperties(subnet) isempty(nps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig) end - all(Graphs.is_tree ∘ Graph ∘ incidencematgraph, subnets) + all(Graphs.is_tree ∘ SimpleGraph ∘ incidencematgraph, subnets) end @doc raw""" @@ -759,43 +759,50 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) pmap = symmap_to_varmap(rs, parametermap) pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap) + # Construct reaction-complex graph complexes, D = reactioncomplexes(rs) - img = incidencematgraph(rn) - undir_img = Graph(incidencematgraph(rn)) - K = ratematrix(rn, pmap) + img = incidencematgraph(rs) + undir_img = SimpleGraph(incidencematgraph(rs)) + K = ratematrix(rs, pmap) - spanning_forest = kruskal_mst(uimg) - undir_spanning_forest = edges(Graph(spanning_forest)) - outofforest_rxns = findall(∉(undir_spanning_forest), edges(undir_img)) + spanning_forest = Graphs.kruskal_mst(undir_img) + outofforest_edges = setdiff(collect(edges(undir_img)), spanning_forest) # Independent Cycle Conditions: for any cycle we create by adding in an out-of-forest reaction, the product of forward reaction rates over the cycle must equal the product of reverse reaction rates over the cycle. - for rxn in outofforest_rxns - g = Graph([undir_spanning_forest..., collect(edges(g))[rxn]]) - ic = cycle_basis(g) + for edge in outofforest_edges + g = SimpleGraph([spanning_forest..., edge]) + ic = Graphs.cycle_basis(g)[1] fwd = prod([K[ic[r], ic[r+1]] for r in 1:length(ic)-1]) * K[ic[end], ic[1]] rev = prod([K[ic[r+1], ic[r]] for r in 1:length(ic)-1]) * K[ic[1], ic[end]] - fwd == rev ? continue : return false + fwd ≈ rev ? continue : return false end - # Spanning Forest Conditions: for non-deficiency 0 networks, we get an additional δ equations. Choose an orientation for each reaction pair in the spanning forest (by default, we will go from smaller to larger indices). + # Spanning Forest Conditions: for non-deficiency 0 networks, we get an additional δ equations. Choose an orientation for each reaction pair in the spanning forest (we will take the one given by default from kruskal_mst). - if deficiency(rn) > 0 - rxn_idxs = findall(∈(spanning_forest), edges(img)) - S_F = netstoichmat(rn)[:, fwdrxn_idxs] + if deficiency(rs) > 0 + rxn_idxs = [edgeindex(D, Graphs.src(e), Graphs.dst(e)) for e in spanning_forest] + S_F = netstoichmat(rs)[:, rxn_idxs] sols = nullspace(S_F) - @assert size(sols, 2) == deficiency(rn) + @assert size(sols, 2) == deficiency(rs) for i in 1:size(sols, 2) α = sols[:, i] - fwd = prod([K[src(e), dst(e)]^α[i] for (e,i) in zip(spanning_forest,1:length(α))]) - rev = prod([K[dst(e), src(e)]^α[i] for (e,i) in zip(spanning_forest, 1:length(α))]) - fwd == rev ? continue : return false + fwd = prod([K[Graphs.src(e), Graphs.dst(e)]^α[i] for (e,i) in zip(spanning_forest,1:length(α))]) + rev = prod([K[Graphs.dst(e), Graphs.src(e)]^α[i] for (e,i) in zip(spanning_forest, 1:length(α))]) + fwd ≈ rev ? continue : return false end end true end +function edgeindex(imat::Matrix{Int64}, src::Int64, dst::Int64) + for i in 1:size(imat, 2) + (imat[src, i] == -1) && (imat[dst, i] == 1) && return i + end + error("This edge does not exist in this reaction graph.") +end + function isdetailedbalanced(rs::ReactionSystem, parametermap::Vector{Pair{Symbol, Float64}}) pdict = Dict(parametermap) isdetailedbalanced(rs, pdict) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index f700bee70b..0cc7a8bd75 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -355,7 +355,7 @@ let k = rand(rng, numparams(rn)) rates = Dict(zip(parameters(rn), k)) @test Catalyst.iscomplexbalanced(rn, rates) == true - @test Catalyst.isdetailedbalanced(rn, rates) == false + @test Catalyst.Catalyst.isdetailedbalanced(rn, rates) == false end ### STRONG LINKAGE CLASS TESTS @@ -643,7 +643,6 @@ end # The following network is conditionally complex balanced - it only - # Reversible, forest-like deficiency zero network - should be detailed balance for any choice of rate constants. let rn = @reaction_network begin @@ -660,11 +659,12 @@ let k2 = rand(StableRNG(232), numparams(rn)) rates2 = Dict(zip(parameters(rn), k2)) - @test isdetailedbalanced(rn, rates1) == true - @test isdetailedbalanced(rn, rates2) == true + rcs, D = reactioncomplexes(rn) + @test Catalyst.isforestlike(rn) == true + @test Catalyst.isdetailedbalanced(rn, rates1) == true + @test Catalyst.isdetailedbalanced(rn, rates2) == true end - # Simple connected reversible network let rn = @reaction_network begin @@ -673,10 +673,11 @@ let (k5, k6), C <--> A end - rates1 = [k1=>1.0, k2=>1.0, k3=>1.0, k4=>1.0, k5=>1.0, k6=>1.0] - @test isdetailedbalanced(rn, rates1) == true - rates2 = [k1=>2.0, k2=>1.0, k3=>1.0, k4=>1.0, k5=>1.0, k6=>1.0] - @test isdetailedbalanced(rn, rates2) == false + rcs, D = reactioncomplexes(rn) + rates1 = [:k1=>1.0, :k2=>1.0, :k3=>1.0, :k4=>1.0, :k5=>1.0, :k6=>1.0] + @test Catalyst.isdetailedbalanced(rn, rates1) == true + rates2 = [:k1=>2.0, :k2=>1.0, :k3=>1.0, :k4=>1.0, :k5=>1.0, :k6=>1.0] + @test Catalyst.isdetailedbalanced(rn, rates2) == false end # Independent cycle tests: the following reaction entwork has 3 out-of-forest reactions. @@ -693,14 +694,17 @@ let (k17, k18), A + E <--> H end + rcs, D = reactioncomplexes(rn) k = rand(rng, numparams(rn)) + p = parameters(rn) rates = Dict(zip(parameters(rn), k)) - @test isdetailedbalanced(rn, rates) == false + @test Catalyst.isdetailedbalanced(rn, rates) == false - rates[k6] = rates[k1]*rates[k4]*rates[k5] / (rates[k2]*rates[k3]) - rates[k14] = rates[k13]*rates[k11]*rates[k8] / (rates[k12]*rates[k7]) - rates[k12] = rates[k8]*rates[k15]*rates[k9]*rates[k11] / (rates[k7]*rates[k16]*rates[k10]) - @test isdetailedbalanced(rn, rates) == true + # Adjust rate constants to obey the independent cycle conditions. + rates[p[6]] = rates[p[1]]*rates[p[4]]*rates[p[5]] / (rates[p[2]]*rates[p[3]]) + rates[p[14]] = rates[p[13]]*rates[p[11]]*rates[p[8]] / (rates[p[12]]*rates[p[7]]) + rates[p[16]] = rates[p[8]]*rates[p[15]]*rates[p[9]]*rates[p[11]] / (rates[p[7]]*rates[p[12]]*rates[p[10]]) + @test Catalyst.isdetailedbalanced(rn, rates) == true end # Deficiency two network: the following reaction network must satisfy both the independent cycle conditions and the spanning forest conditions @@ -713,21 +717,27 @@ let (k9, k10), 3A <--> 3B end + rcs, D = reactioncomplexes(rn) + @test Catalyst.edgeindex(D, 1, 2) == 1 + @test Catalyst.edgeindex(D, 4, 3) == 6 k = rand(rng, numparams(rn)) + p = parameters(rn) rates = Dict(zip(parameters(rn), k)) - @test isdetailedbalanced(rn, rates) == false - - rates[k8] = rates[k7]*rates[k5]*rates[k9] / (rates[k6]*rates[k10]) - rates[k3] = rates[k2]*rates[k4]*rates[k9] / (rates[k1]*rates[k10]) - @test isdetailedbalanced(rn, rates) == false - - cons = rates[k6] / rates[k5] - rates[k1] = rates[k2] * cons - rates[k9] = rates[k10] * cons^(3/2) - rates[k8] = rates[k7]*rates[k5]*rates[k9] / (rates[k6]*rates[k10]) - rates[k3] = rates[k2]*rates[k4]*rates[k9] / (rates[k1]*rates[k10]) - @test isdeatiledbalanced(rn, rates) == true - @test isdetailedbalanced(rn, rates) == false + @test Catalyst.isdetailedbalanced(rn, rates) == false + + # Adjust rate constants to fulfill independent cycle conditions. + rates[p[8]] = rates[p[7]]*rates[p[5]]*rates[p[9]] / (rates[p[6]]*rates[p[10]]) + rates[p[3]] = rates[p[2]]*rates[p[4]]*rates[p[9]] / (rates[p[1]]*rates[p[10]]) + @test Catalyst.isdetailedbalanced(rn, rates) == false + # Should still fail - doesn't satisfy spanning forest conditions. + + # Adjust rate constants to fulfill spanning forest conditions. + cons = rates[p[6]] / rates[p[5]] + rates[p[1]] = rates[p[2]] * cons + rates[p[9]] = rates[p[10]] * cons^(3/2) + rates[p[8]] = rates[p[7]]*rates[p[5]]*rates[p[9]] / (rates[p[6]]*rates[p[10]]) + rates[p[3]] = rates[p[2]]*rates[p[4]]*rates[p[9]] / (rates[p[1]]*rates[p[10]]) + @test Catalyst.isdetailedbalanced(rn, rates) == true end ### DEFICIENCY ONE TESTS From e081ebf00253aef805b93ec0df2f3743161f8485 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 2 Jul 2024 20:23:38 -0400 Subject: [PATCH 22/29] formatter --- src/network_analysis.jl | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 9d6d47605e..dc1f5f70c1 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -398,7 +398,7 @@ function isterminal(lc::Vector, rn::ReactionSystem) true end -function isforestlike(rn::ReactionSystem) +function isforestlike(rn::ReactionSystem) subnets = subnetworks(rn) im = get_networkproperties(rn).incidencemat @@ -745,7 +745,7 @@ Constructively compute whether a kinetic system (a reaction network with a set o solutions, using the Wegscheider conditions, [Feinberg, 1989](https://www.sciencedirect.com/science/article/pii/0009250989851243). A detailed-balanced solution is one for which the rate of every forward reaction exactly equals its reverse reaction. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. """ -function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) +function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) if length(parametermap) != numparams(rs) error("Incorrect number of parameters specified.") elseif !isreversible(rs) @@ -758,7 +758,7 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) pmap = symmap_to_varmap(rs, parametermap) pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap) - + # Construct reaction-complex graph complexes, D = reactioncomplexes(rs) img = incidencematgraph(rs) @@ -772,23 +772,25 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) for edge in outofforest_edges g = SimpleGraph([spanning_forest..., edge]) ic = Graphs.cycle_basis(g)[1] - fwd = prod([K[ic[r], ic[r+1]] for r in 1:length(ic)-1]) * K[ic[end], ic[1]] - rev = prod([K[ic[r+1], ic[r]] for r in 1:length(ic)-1]) * K[ic[1], ic[end]] + fwd = prod([K[ic[r], ic[r + 1]] for r in 1:(length(ic) - 1)]) * K[ic[end], ic[1]] + rev = prod([K[ic[r + 1], ic[r]] for r in 1:(length(ic) - 1)]) * K[ic[1], ic[end]] fwd ≈ rev ? continue : return false end - + # Spanning Forest Conditions: for non-deficiency 0 networks, we get an additional δ equations. Choose an orientation for each reaction pair in the spanning forest (we will take the one given by default from kruskal_mst). - + if deficiency(rs) > 0 rxn_idxs = [edgeindex(D, Graphs.src(e), Graphs.dst(e)) for e in spanning_forest] S_F = netstoichmat(rs)[:, rxn_idxs] - sols = nullspace(S_F) + sols = nullspace(S_F) @assert size(sols, 2) == deficiency(rs) for i in 1:size(sols, 2) α = sols[:, i] - fwd = prod([K[Graphs.src(e), Graphs.dst(e)]^α[i] for (e,i) in zip(spanning_forest,1:length(α))]) - rev = prod([K[Graphs.dst(e), Graphs.src(e)]^α[i] for (e,i) in zip(spanning_forest, 1:length(α))]) + fwd = prod([K[Graphs.src(e), Graphs.dst(e)]^α[i] + for (e, i) in zip(spanning_forest, 1:length(α))]) + rev = prod([K[Graphs.dst(e), Graphs.src(e)]^α[i] + for (e, i) in zip(spanning_forest, 1:length(α))]) fwd ≈ rev ? continue : return false end end @@ -796,7 +798,7 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) true end -function edgeindex(imat::Matrix{Int64}, src::Int64, dst::Int64) +function edgeindex(imat::Matrix{Int64}, src::Int64, dst::Int64) for i in 1:size(imat, 2) (imat[src, i] == -1) && (imat[dst, i] == 1) && return i end @@ -817,7 +819,6 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap) error("Parameter map must be a dictionary, tuple, or vector of symbol/value pairs.") end - """ iscomplexbalanced(rs::ReactionSystem, parametermap) @@ -1120,4 +1121,3 @@ function robustspecies(rn::ReactionSystem) nps.robustspecies end - From 2b35c94deaa8b04d8d58c46b4016ae34569bfc1c Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 16 Jul 2024 12:41:47 -0400 Subject: [PATCH 23/29] up --- src/network_analysis.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index dc1f5f70c1..e4930ecf3b 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -400,11 +400,10 @@ end function isforestlike(rn::ReactionSystem) subnets = subnetworks(rn) + nps = get_networkproperties(rn) - im = get_networkproperties(rn).incidencemat - isempty(im) && - error("Error, please call reactioncomplexes(rn::ReactionSystem) to ensure the incidence matrix has been cached.") - sparseig = issparse(im) + isempty(nps.incidencemat) && reactioncomplexes(rn) + sparseig = issparse(incidencemat) for subnet in subnets nps = get_networkproperties(subnet) isempty(nps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig) From d935dd95432ee10da54650475d46735826589498 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 17 Jul 2024 14:36:23 -0400 Subject: [PATCH 24/29] merge fix --- 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 e4930ecf3b..6427b152b8 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -821,8 +821,9 @@ end """ iscomplexbalanced(rs::ReactionSystem, parametermap) -Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit complex-balanced equilibrium -solutions, following the method in van der Schaft et al., [2015](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. +Constructively compute whether a network will have complex-balanced equilibrium +solutions, following the method in van der Schaft et al., [2015](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). +Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. """ function iscomplexbalanced(rs::ReactionSystem, parametermap::Dict) From 7025afec5a9f4f9e39c006ef5010215b79c78c5b Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 17 Jul 2024 14:37:27 -0400 Subject: [PATCH 25/29] merge fix --- test/network_analysis/network_properties.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 0cc7a8bd75..42591bcd61 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -355,7 +355,7 @@ let k = rand(rng, numparams(rn)) rates = Dict(zip(parameters(rn), k)) @test Catalyst.iscomplexbalanced(rn, rates) == true - @test Catalyst.Catalyst.isdetailedbalanced(rn, rates) == false + @test Catalyst.isdetailedbalanced(rn, rates) == false end ### STRONG LINKAGE CLASS TESTS From f9beabc30d1e91b3cbf569e0a8ce96684fbfe15c Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 17 Jul 2024 14:57:39 -0400 Subject: [PATCH 26/29] 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 6427b152b8..a51badb531 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -403,7 +403,7 @@ function isforestlike(rn::ReactionSystem) nps = get_networkproperties(rn) isempty(nps.incidencemat) && reactioncomplexes(rn) - sparseig = issparse(incidencemat) + sparseig = issparse(nps.incidencemat) for subnet in subnets nps = get_networkproperties(subnet) isempty(nps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig) From a2a9e4d616da580f73201ae0d09287267a7279e1 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 4 Nov 2024 14:05:58 -0500 Subject: [PATCH 27/29] minor fixes --- src/network_analysis.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index a51badb531..f122f2388e 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -738,7 +738,7 @@ function conservationlaw_errorcheck(rs, pre_varmap) end """ - iscomplexbalanced(rs::ReactionSystem, parametermap) + isdetailedbalanced(rs::ReactionSystem, parametermap) Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit detailed-balanced equilibrium solutions, using the Wegscheider conditions, [Feinberg, 1989](https://www.sciencedirect.com/science/article/pii/0009250989851243). A detailed-balanced solution is one for which the rate of every forward reaction exactly equals its reverse reaction. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. @@ -782,7 +782,6 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) rxn_idxs = [edgeindex(D, Graphs.src(e), Graphs.dst(e)) for e in spanning_forest] S_F = netstoichmat(rs)[:, rxn_idxs] sols = nullspace(S_F) - @assert size(sols, 2) == deficiency(rs) for i in 1:size(sols, 2) α = sols[:, i] @@ -797,19 +796,20 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) true end -function edgeindex(imat::Matrix{Int64}, src::Int64, dst::Int64) +# Helper to find the index of the reaction with a given reactnat and product complex. +function edgeindex(imat, src::T, dst::T) where T <: Int for i in 1:size(imat, 2) (imat[src, i] == -1) && (imat[dst, i] == 1) && return i end error("This edge does not exist in this reaction graph.") end -function isdetailedbalanced(rs::ReactionSystem, parametermap::Vector{Pair{Symbol, Float64}}) +function isdetailedbalanced(rs::ReactionSystem, parametermap::Vector{<:Pair}) pdict = Dict(parametermap) isdetailedbalanced(rs, pdict) end -function isdetailedbalanced(rs::ReactionSystem, parametermap::Tuple{Pair{Symbol, Float64}}) +function isdetailedbalanced(rs::ReactionSystem, parametermap::Tuple{<:Pair}) pdict = Dict(parametermap) isdetailedbalanced(rs, pdict) end From 11c13131c126f3b4d05265ce1761d8bf60dbacb0 Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 5 Nov 2024 12:34:44 -0500 Subject: [PATCH 28/29] tolerance and nullspace --- src/network_analysis.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index f122f2388e..d2c3cc735d 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -738,13 +738,13 @@ function conservationlaw_errorcheck(rs, pre_varmap) end """ - isdetailedbalanced(rs::ReactionSystem, parametermap) + isdetailedbalanced(rs::ReactionSystem, parametermap; reltol=1e-9, abstol) Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit detailed-balanced equilibrium solutions, using the Wegscheider conditions, [Feinberg, 1989](https://www.sciencedirect.com/science/article/pii/0009250989851243). A detailed-balanced solution is one for which the rate of every forward reaction exactly equals its reverse reaction. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. """ -function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) +function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict; abstol=0, reltol=1e-9) if length(parametermap) != numparams(rs) error("Incorrect number of parameters specified.") elseif !isreversible(rs) @@ -773,7 +773,7 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) ic = Graphs.cycle_basis(g)[1] fwd = prod([K[ic[r], ic[r + 1]] for r in 1:(length(ic) - 1)]) * K[ic[end], ic[1]] rev = prod([K[ic[r + 1], ic[r]] for r in 1:(length(ic) - 1)]) * K[ic[1], ic[end]] - fwd ≈ rev ? continue : return false + isapprox(fwd, rev; atol = abstol, rtol = reltol) ? continue : return false end # Spanning Forest Conditions: for non-deficiency 0 networks, we get an additional δ equations. Choose an orientation for each reaction pair in the spanning forest (we will take the one given by default from kruskal_mst). @@ -781,7 +781,7 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) if deficiency(rs) > 0 rxn_idxs = [edgeindex(D, Graphs.src(e), Graphs.dst(e)) for e in spanning_forest] S_F = netstoichmat(rs)[:, rxn_idxs] - sols = nullspace(S_F) + sols = positive_nullspace(S_F) for i in 1:size(sols, 2) α = sols[:, i] @@ -789,14 +789,14 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict) for (e, i) in zip(spanning_forest, 1:length(α))]) rev = prod([K[Graphs.dst(e), Graphs.src(e)]^α[i] for (e, i) in zip(spanning_forest, 1:length(α))]) - fwd ≈ rev ? continue : return false + isapprox(fwd, rev; atol = abstol, rtol = reltol) ? continue : return false end end true end -# Helper to find the index of the reaction with a given reactnat and product complex. +# Helper to find the index of the reaction with a given reactant and product complex. function edgeindex(imat, src::T, dst::T) where T <: Int for i in 1:size(imat, 2) (imat[src, i] == -1) && (imat[dst, i] == 1) && return i From af16c6998647123a5d5a60c5c07638dc2e5ba840 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 8 Nov 2024 16:04:10 -0500 Subject: [PATCH 29/29] fix docstrings --- src/network_analysis.jl | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index d2c3cc735d..93bebe1888 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -743,7 +743,6 @@ end Constructively compute whether a kinetic system (a reaction network with a set of rate constants) will admit detailed-balanced equilibrium solutions, using the Wegscheider conditions, [Feinberg, 1989](https://www.sciencedirect.com/science/article/pii/0009250989851243). A detailed-balanced solution is one for which the rate of every forward reaction exactly equals its reverse reaction. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. """ - function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict; abstol=0, reltol=1e-9) if length(parametermap) != numparams(rs) error("Incorrect number of parameters specified.") @@ -825,7 +824,6 @@ Constructively compute whether a network will have complex-balanced equilibrium solutions, following the method in van der Schaft et al., [2015](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. """ - function iscomplexbalanced(rs::ReactionSystem, parametermap::Dict) if length(parametermap) != numparams(rs) error("Incorrect number of parameters specified.") @@ -902,7 +900,6 @@ end constant of the reaction between complex i and complex j. Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. """ - function ratematrix(rs::ReactionSystem, rates::Vector{Float64}) complexes, D = reactioncomplexes(rs) n = length(complexes) @@ -997,7 +994,6 @@ end Returns the matrix of a basis of cycles (or flux vectors), or a basis for reaction fluxes for which the system is at steady state. These correspond to right eigenvectors of the stoichiometric matrix. Equivalent to [`fluxmodebasis`](@ref). """ - function cycles(rs::ReactionSystem) nps = get_networkproperties(rs) nsm = netstoichmat(rs) @@ -1032,7 +1028,6 @@ end See documentation for [`cycles`](@ref). """ - function fluxvectors(rs::ReactionSystem) cycles(rs) end @@ -1044,7 +1039,6 @@ 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) all(r -> ismassaction(r, rn), reactions(rn)) || error("The deficiency one theorem is only valid for reaction networks that are mass action.") @@ -1067,7 +1061,6 @@ 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) all(r -> ismassaction(r, rn), reactions(rn)) || error("The deficiency zero theorem is only valid for reaction networks that are mass action.") @@ -1082,7 +1075,6 @@ end 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) complexes, D = reactioncomplexes(rn) nps = get_networkproperties(rn)