From e3c9d339df5cfc70a39519d6ab85beb7d4aea902 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 22 Nov 2024 16:37:51 -0500 Subject: [PATCH 01/29] init --- src/network_analysis.jl | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 86802fb919..7f14076a50 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -192,6 +192,23 @@ function complexstoichmat(::Type{Matrix{Int}}, rn::ReactionSystem, rcs) Z end +""" + networklaplacian(rn::ReactionSystem; parammap=nothing, sparse=false) + + Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as $$\frac{dx}{dt} = Y A_k Φ(x)$$, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the graph Laplacian. This is an n-by-n matrix, where n is the number of complexes, where A_{ij} = k_{ij} if a reaction exists between the two complexes and 0 otherwise. +""" +function networklaplacian(rn::ReactionSystem; parammap = nothing, sparse = false) + +end + + +""" + monomialvector(rn::ReactionSystem) +""" +function monomialvector(rn::ReactionSystem) + +end + @doc raw""" complexoutgoingmat(network::ReactionSystem; sparse=false) From 17d1cd531afb50dd35b5ce6e40a3a8a5322a0320 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 22 Nov 2024 20:31:28 -0500 Subject: [PATCH 02/29] implement functions --- src/network_analysis.jl | 112 +++++++++++++++++++++++++++++++++++++--- src/reactionsystem.jl | 1 + 2 files changed, 105 insertions(+), 8 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 7f14076a50..7a52c931e0 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -193,20 +193,116 @@ function complexstoichmat(::Type{Matrix{Int}}, rn::ReactionSystem, rcs) end """ - networklaplacian(rn::ReactionSystem; parammap=nothing, sparse=false) + laplacianmat(rn::ReactionSystem, parammap=nothing; sparse=false) - Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as $$\frac{dx}{dt} = Y A_k Φ(x)$$, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the graph Laplacian. This is an n-by-n matrix, where n is the number of complexes, where A_{ij} = k_{ij} if a reaction exists between the two complexes and 0 otherwise. + Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the graph Laplacian, and Φ is the [`massactionvector`](@ref). A_k is an n-by-n matrix, where n is the number of complexes, where A_{ij} = k_{ij} if a reaction exists between the two complexes and 0 otherwise. + Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via parammap. """ -function networklaplacian(rn::ReactionSystem; parammap = nothing, sparse = false) - +function laplacianmat(rn::ReactionSystem, pmap = Dict(); sparse = false) + D = incidencemat(rn; sparse); K = fluxmat(rn, pmap; sparse) + -D*K +end + +""" + fluxmat(rn::ReactionSystem, parammap = nothing; sparse=true) + + Return an r×c matrix K such that, if complex j is the substrate complex of reaction i, then K_{ij} = k, the rate constant for this reaction. Mostly a helper function for the network Laplacian, [`networklaplacianmat`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix. +""" +function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=true) + rates = reactionrates(rn) + + if !isempty(pmap) + length(pmap) != length(parameters(rn)) && error("Incorrect number of parameters were specified.") + pmap = symmap_to_varmap(rs, pmap) + pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap) + rates = [substitute(rate, pmap) for rate in rates] + end + + rcmap = reactioncomplexmap(rn); nc = length(rcmap); nr = length(rates) + mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates) + K = if sparse + fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates) + else + fluxmat(Matrix{mtype}, rcmap, rates) + end +end + +function fluxmat(::Type{SparseMatrixCSC{T, Int}}, rcmap, rates) where T + Is = Int[] + Js = Int[] + Vs = T[] + for (i, (complex, rxs)) in enumerate(rcmap) + for (rx, dir) in rxs + dir == -1 && begin + push!(Is, rx) + push!(Js, i) + push!(Vs, rates[rx]) + end + end + end + Z = sparse(Is, Js, Vs, length(rates), length(rcmap)) +end + +function fluxmat(::Type{Matrix{T}}, rcmap, rates) where T + nr = length(rates); nc = length(rcmap) + K = zeros(T, nr, nc) + for (i, (complex, rxs)) in enumerate(rcmap) + for (rx, dir) in rxs + dir == -1 && (K[rx, i] = rates[rx]) + end + end + K end +function fluxmat(rn::ReactionSystem, pmap::Vector; sparse = false) + pdict = Dict(pmap) + fluxmat(rs, pdict; sparse) +end + +function fluxmat(rn::ReactionSystem, pmap::Tuple; sparse = false) + pdict = Dict(pmap) + fluxmat(rs, pdict; sparse) +end """ - monomialvector(rn::ReactionSystem) + massactionvector(rn::ReactionSystem, scmap = nothing) + + Return the vector whose entries correspond to the "substrate complexes" of each reaction. For example, given the reaction A + B --> C, the corresponding entry of the vector would be A*B, and for the reaction 2X + Y --> 3X, the corresponding entry would be X^2*Y. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the [`networklaplacian`](@ref). This utility returns Φ(x). + Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap. """ -function monomialvector(rn::ReactionSystem) - +function massactionvector(rn::ReactionSystem, scmap::Dict = Dict()) + r = numreactions(rn); rxs = reactions(rn) + sm = speciesmap(rn); specs = species(rn) + + if !all(r -> ismassaction(r, rn), rxs) + error("The supplied ReactionSystem has reactions that are not ismassaction. The monomial vector is only defined for pure mass action networks.") + end + + if !isempty(scmap) + length(scmap) != length(specs) && error("Incorrect number of species concentrations were specified.") + scmap = symmap_to_varmap(rs, scmap) + scmap = Dict(ModelingToolkit.value(k) => v for (k, v) in scmap) + specs = [substitute(s, scmap) for s in specs] + end + + Φ = Vector{eltype(specs)}() + for i in 1:n + subs = rxs[i].substrates + stoich = rxs[i].substoich + push!(Φ, prod([specs[sm[s]]^α for (s, α) in zip(subs, stoich)])) + end + + Φ +end + +function massactionvector(rn::ReactionSystem, scmap::Tuple) + sdict = Dict(scmap) + massactionvector(rs, sdict) +end + +function massactionvector(rn::ReactionSystem, scmap::Vector) + sdict = Dict(scmap) + massactionvector(rs, sdict) end @doc raw""" @@ -804,7 +900,7 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict; abstol=0, re 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.") + error("The supplied ReactionSystem has reactions that are not ismassaction. Testing for being detailed balanced is currently only supported for pure mass action networks.") end isforestlike(rs) && deficiency(rs) == 0 && return true diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 83ed069ea0..6d24d28315 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -92,6 +92,7 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re incidencemat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0) complexstoichmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0) complexoutgoingmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0) + networklaplacian::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0) incidencegraph::Graphs.SimpleDiGraph{Int} = Graphs.DiGraph() linkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0) stronglinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0) From 9e2370584afecdc4981e3cf1809aa3099dece42f Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 22 Nov 2024 22:33:27 -0500 Subject: [PATCH 03/29] tests --- src/network_analysis.jl | 28 +- test/network_analysis/crn_theory.jl | 474 ++++++++++++++ test/network_analysis/network_properties.jl | 682 +++----------------- test/runtests.jl | 1 + 4 files changed, 577 insertions(+), 608 deletions(-) create mode 100644 test/network_analysis/crn_theory.jl diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 7a52c931e0..873922be66 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -200,7 +200,7 @@ end """ function laplacianmat(rn::ReactionSystem, pmap = Dict(); sparse = false) D = incidencemat(rn; sparse); K = fluxmat(rn, pmap; sparse) - -D*K + D*K end """ @@ -213,7 +213,7 @@ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=true) if !isempty(pmap) length(pmap) != length(parameters(rn)) && error("Incorrect number of parameters were specified.") - pmap = symmap_to_varmap(rs, pmap) + pmap = symmap_to_varmap(rn, pmap) pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap) rates = [substitute(rate, pmap) for rate in rates] end @@ -256,18 +256,18 @@ end function fluxmat(rn::ReactionSystem, pmap::Vector; sparse = false) pdict = Dict(pmap) - fluxmat(rs, pdict; sparse) + fluxmat(rn, pdict; sparse) end function fluxmat(rn::ReactionSystem, pmap::Tuple; sparse = false) pdict = Dict(pmap) - fluxmat(rs, pdict; sparse) + fluxmat(rn, pdict; sparse) end """ massactionvector(rn::ReactionSystem, scmap = nothing) - Return the vector whose entries correspond to the "substrate complexes" of each reaction. For example, given the reaction A + B --> C, the corresponding entry of the vector would be A*B, and for the reaction 2X + Y --> 3X, the corresponding entry would be X^2*Y. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the [`networklaplacian`](@ref). This utility returns Φ(x). + Return the vector whose entries correspond to the "mass action products" of each complex. For example, given the complex A + B, the corresponding entry of the vector would be A*B, and for the complex 2X + Y, the corresponding entry would be X^2*Y. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the [`networklaplacian`](@ref). This utility returns Φ(x). Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap. """ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict()) @@ -280,16 +280,18 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict()) if !isempty(scmap) length(scmap) != length(specs) && error("Incorrect number of species concentrations were specified.") - scmap = symmap_to_varmap(rs, scmap) + scmap = symmap_to_varmap(rn, scmap) scmap = Dict(ModelingToolkit.value(k) => v for (k, v) in scmap) specs = [substitute(s, scmap) for s in specs] end - Φ = Vector{eltype(specs)}() - for i in 1:n - subs = rxs[i].substrates - stoich = rxs[i].substoich - push!(Φ, prod([specs[sm[s]]^α for (s, α) in zip(subs, stoich)])) + vtype = eltype(specs) <: Symbolics.BasicSymbolic ? Num : eltype(specs) + Φ = Vector{vtype}() + rcmap = reactioncomplexmap(rn) + for comp in keys(reactioncomplexmap(rn)) + subs = map(ce -> getfield(ce, :speciesid), comp) + stoich = map(ce -> getfield(ce, :speciesstoich), comp) + push!(Φ, prod(vtype[specs[s]^α for (s, α) in zip(subs, stoich)])) end Φ @@ -297,12 +299,12 @@ end function massactionvector(rn::ReactionSystem, scmap::Tuple) sdict = Dict(scmap) - massactionvector(rs, sdict) + massactionvector(rn, sdict) end function massactionvector(rn::ReactionSystem, scmap::Vector) sdict = Dict(scmap) - massactionvector(rs, sdict) + massactionvector(rn, sdict) end @doc raw""" diff --git a/test/network_analysis/crn_theory.jl b/test/network_analysis/crn_theory.jl new file mode 100644 index 0000000000..f2f40040c9 --- /dev/null +++ b/test/network_analysis/crn_theory.jl @@ -0,0 +1,474 @@ +# Tests for properties from chemical reaction network theory: deficiency theorems, complex/detailed balance, etc. + +# Tests that `iscomplexbalanced` works for different rate inputs. +# Tests that non-valid rate input yields and error +let + # Declares network. + rn = @reaction_network begin + k1, 3A + 2B --> 3C + k2, B + 4D --> 2E + k3, 2E --> 3C + (k4, k5), B + 4D <--> 3A + 2B + k6, F --> B + 4D + k7, 3C --> F + end + + # Declares rate alternatives. + k = rand(rng, numparams(rn)) + rates_vec = Pair.(parameters(rn), k) + rates_tup = Tuple(rates_vec) + rates_dict = Dict(rates_vec) + rates_invalid = k + + # Tests that inputs are handled correctly. + @test Catalyst.iscomplexbalanced(rn, rates_vec) == Catalyst.iscomplexbalanced(rn, rates_tup) + @test Catalyst.iscomplexbalanced(rn, rates_tup) == Catalyst.iscomplexbalanced(rn, rates_dict) + @test_throws Exception Catalyst.iscomplexbalanced(rn, k) +end + +# Tests rate matrix computation for various input types. +let + # Declares network and its known rate matrix. + rn = @reaction_network begin + (k2, k1), A1 <--> A2 + A3 + k3, A2 + A3 --> A4 + k4, A4 --> A5 + (k6, k5), A5 <--> 2A6 + k7, 2A6 --> A4 + k8, A4 + A5 --> A7 + end + rate_mat = [ + 0.0 1.0 0.0 0.0 0.0 0.0 0.0; + 2.0 0.0 3.0 0.0 0.0 0.0 0.0; + 0.0 0.0 0.0 4.0 0.0 0.0 0.0; + 0.0 0.0 0.0 0.0 5.0 0.0 0.0; + 0.0 0.0 7.0 6.0 0.0 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 0.0 8.0; + 0.0 0.0 0.0 0.0 0.0 0.0 0.0; + ] + + # Declares rate alternatives. + rate_vals = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0] + rates_vec = Pair.(parameters(rn), rate_vals) + rates_tup = Tuple(rates_vec) + rates_dict = Dict(rates_vec) + rates_invalid = reshape(rate_vals, 1, 8) + + # Tests that all input types generates the correct rate matrix. + Catalyst.ratematrix(rn, rate_vals) == rate_mat + Catalyst.ratematrix(rn, rates_vec) == rate_mat + Catalyst.ratematrix(rn, rates_tup) == rate_mat + Catalyst.ratematrix(rn, rates_dict) == rate_mat + @test_throws Exception Catalyst.iscomplexbalanced(rn, rates_invalid) +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 + 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 + + +### Complex balance and reversibility tests ### + +# Test function. +function testreversibility(rn, B, rev, weak_rev) + @test isreversible(rn) == rev + subrn = subnetworks(rn) + @test isweaklyreversible(rn, subrn) == weak_rev +end + +# Tests reversibility for networks with known reversibility. +let + rn = @reaction_network begin + (k2, k1), A1 <--> A2 + A3 + k3, A2 + A3 --> A4 + k4, A4 --> A5 + (k6, k5), A5 <--> 2A6 + k7, 2A6 --> A4 + k8, A4 + A5 --> A7 + end + rev = false + weak_rev = false + testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test Catalyst.iscomplexbalanced(rn, rates) == false +end + +let + rn4 = @reaction_network begin + (k1, k2), C1 <--> C2 + (k3, k4), C2 <--> C3 + (k5, k6), C3 <--> C1 + end + + k = rand(rng, numparams(rn4)) + rates = Dict(zip(parameters(rn4), k)) + @test Catalyst.iscomplexbalanced(rn4, rates) == true +end + +let + rn = @reaction_network begin + (k2, k1), A1 <--> A2 + A3 + k3, A2 + A3 --> A4 + k4, A4 --> A5 + (k6, k5), A5 <--> 2A6 + k7, A4 --> 2A6 + (k9, k8), A4 + A5 <--> A7 + end + rev = false + weak_rev = false + testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test Catalyst.iscomplexbalanced(rn, rates) == false +end + +let + rn = @reaction_network begin + k1, A --> B + k2, A --> C + end + rev = false + weak_rev = false + testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test Catalyst.iscomplexbalanced(rn, rates) == false +end + +let + rn = @reaction_network begin + k1, A --> B + k2, A --> C + k3, B + C --> 2A + end + rev = false + weak_rev = false + testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test Catalyst.iscomplexbalanced(rn, rates) == false +end + +let + rn = @reaction_network begin + (k2, k1), A <--> 2B + (k4, k3), A + C <--> D + k5, D --> B + E + k6, B + E --> A + C + end + rev = false + weak_rev = true + testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test Catalyst.iscomplexbalanced(rn, rates) == true +end + +let + rn = @reaction_network begin + (k2, k1), A + E <--> AE + k3, AE --> B + E + end + rev = false + weak_rev = false + testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test Catalyst.iscomplexbalanced(rn, rates) == false +end + +let + rn = @reaction_network begin + (k2, k1), A + E <--> AE + (k4, k3), AE <--> B + E + end + rev = true + weak_rev = true + testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test Catalyst.iscomplexbalanced(rn, rates) == true +end + +let + rn = @reaction_network begin (k2, k1), A + B <--> 2A end + rev = true + weak_rev = true + testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test Catalyst.iscomplexbalanced(rn, rates) == true +end + +let + rn = @reaction_network begin + k1, A + B --> 3A + k2, 3A --> 2A + C + k3, 2A + C --> 2B + k4, 2B --> A + B + end + rev = false + weak_rev = true + testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test Catalyst.iscomplexbalanced(rn, rates) == true +end + +let + rn = @reaction_network begin + (k2, k1), A + E <--> AE + (k4, k3), AE <--> B + E + k5, B --> 0 + k6, 0 --> A + end + rev = false + weak_rev = false + testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test Catalyst.iscomplexbalanced(rn, rates) == false +end + +let + rn = @reaction_network begin + k1, 3A + 2B --> 3C + k2, B + 4D --> 2E + k3, 2E --> 3C + (k4, k5), B + 4D <--> 3A + 2B + k6, F --> B + 4D + k7, 3C --> F + end + + k = rand(rng, numparams(rn)) + rates = Dict(zip(parameters(rn), k)) + @test Catalyst.iscomplexbalanced(rn, rates) == true + @test Catalyst.isdetailedbalanced(rn, rates) == false +end + + +### DEFICIENCY THEOREMS TESTS + +# 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 + end + + @test Catalyst.satisfiesdeficiencyone(rn) == false +end + +# Fails because linkage deficiencies do not sum to total deficiency +let + rn = @reaction_network begin + (k1, k2), A <--> 2A + (k3, k4), A + B <--> C + (k5, k6), C <--> B + end + + @test Catalyst.satisfiesdeficiencyone(rn) == false +end + +# Fails because a linkage class has deficiency two +let + 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 + +### 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, 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 Catalyst.satisfiesdeficiencyzero(rn) == true + @test Catalyst.satisfiesdeficiencyzero(rn2) == false + @test Catalyst.satisfiesdeficiencyzero(rn3) == false + @test Catalyst.satisfiesdeficiencyzero(rn4) == false +end + +### 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)) + + 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 + (k1, k2), A <--> B + (k3, k4), B <--> C + (k5, k6), C <--> A + end + + 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. +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 + + rcs, D = reactioncomplexes(rn) + k = rand(rng, numparams(rn)) + p = parameters(rn) + rates = Dict(zip(parameters(rn), k)) + @test Catalyst.isdetailedbalanced(rn, rates) == false + + # 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 +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 + + 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 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 diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 21f7cb402e..7ffbc166ed 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -1,4 +1,5 @@ ### Prepares Tests ### +# Tests for network structural information: associated matrices, graphs, linkage classes, etc. # Fetch packages. using Catalyst, LinearAlgebra, Test, SparseArrays @@ -7,6 +8,8 @@ using Catalyst, LinearAlgebra, Test, SparseArrays using StableRNGs rng = StableRNG(514) +### Defining reaction networks. + ### Basic Tests ### # Tests basic `ReactionComplex` properties. @@ -59,21 +62,7 @@ let k = rand(rng, numparams(MAPK)) rates = Dict(zip(parameters(MAPK), k)) @test Catalyst.iscomplexbalanced(MAPK, rates) == false - # i=0; - # for lcs in linkageclasses(MAPK) - # i=i+1 - # println("Linkage no ",i) - # for comps in rcs[lcs] - # if comps.speciesids ≠ Int64[] - # println(sum(species(rn2)[comps.speciesids])) - # else - # println("0") - # end - # end - # println("-----------") - # end - - # Testing if cycles identifies reversible reactions as cycles (one forward, one reverse) + cyclemat = Catalyst.cycles(MAPK) S = netstoichmat(MAPK) for i in 1:size(S, 2)-1 @@ -110,19 +99,6 @@ let k = rand(rng, numparams(rn2)) rates = Dict(zip(parameters(rn2), k)) @test Catalyst.iscomplexbalanced(rn2, rates) == false - # i=0; - # for lcs in linkageclasses(rn2) - # i=i+1 - # println("Linkage no ",i) - # for comps in rcs[lcs] - # if comps.speciesids ≠ Int64[] - # println(sum(species(rn2)[comps.speciesids])) - # else - # println("0") - # end - # end - # println("-----------") - # end end # Tests network analysis functions on third network (by comparing to manually computed outputs). @@ -154,208 +130,6 @@ let k = rand(rng, numparams(rn3)) rates = Dict(zip(parameters(rn3), k)) @test Catalyst.iscomplexbalanced(rn3, rates) == false - # i=0; - # for lcs in linkageclasses(rn3) - # i=i+1 - # println("Linkage no ",i) - # for comps in rcs[lcs] - # if comps.speciesids ≠ Int64[] - # println(sum(species(rn3)[comps.speciesids])) - # else - # println("0") - # end - # end - # println("-----------") - # end -end - -let - rn4 = @reaction_network begin - (k1, k2), C1 <--> C2 - (k3, k4), C2 <--> C3 - (k5, k6), C3 <--> C1 - end - - k = rand(rng, numparams(rn4)) - rates = Dict(zip(parameters(rn4), k)) - @test Catalyst.iscomplexbalanced(rn4, rates) == true -end - -### Tests Reversibility ### - -# Test function. -function testreversibility(rn, B, rev, weak_rev) - @test isreversible(rn) == rev - subrn = subnetworks(rn) - @test isweaklyreversible(rn, subrn) == weak_rev -end - -# Tests reversibility for networks with known reversibility. -let - rn = @reaction_network begin - (k2, k1), A1 <--> A2 + A3 - k3, A2 + A3 --> A4 - k4, A4 --> A5 - (k6, k5), A5 <--> 2A6 - k7, 2A6 --> A4 - k8, A4 + A5 --> A7 - end - rev = false - weak_rev = false - testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) - - k = rand(rng, numparams(rn)) - rates = Dict(zip(parameters(rn), k)) - @test Catalyst.iscomplexbalanced(rn, rates) == false -end - -let - rn = @reaction_network begin - (k2, k1), A1 <--> A2 + A3 - k3, A2 + A3 --> A4 - k4, A4 --> A5 - (k6, k5), A5 <--> 2A6 - k7, A4 --> 2A6 - (k9, k8), A4 + A5 <--> A7 - end - rev = false - weak_rev = false - testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) - - k = rand(rng, numparams(rn)) - rates = Dict(zip(parameters(rn), k)) - @test Catalyst.iscomplexbalanced(rn, rates) == false -end - -let - rn = @reaction_network begin - k1, A --> B - k2, A --> C - end - rev = false - weak_rev = false - testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) - k = rand(rng, numparams(rn)) - rates = Dict(zip(parameters(rn), k)) - @test Catalyst.iscomplexbalanced(rn, rates) == false -end - -let - rn = @reaction_network begin - k1, A --> B - k2, A --> C - k3, B + C --> 2A - end - rev = false - weak_rev = false - testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) - - k = rand(rng, numparams(rn)) - rates = Dict(zip(parameters(rn), k)) - @test Catalyst.iscomplexbalanced(rn, rates) == false -end - -let - rn = @reaction_network begin - (k2, k1), A <--> 2B - (k4, k3), A + C <--> D - k5, D --> B + E - k6, B + E --> A + C - end - rev = false - weak_rev = true - testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) - - k = rand(rng, numparams(rn)) - rates = Dict(zip(parameters(rn), k)) - @test Catalyst.iscomplexbalanced(rn, rates) == true -end - -let - rn = @reaction_network begin - (k2, k1), A + E <--> AE - k3, AE --> B + E - end - rev = false - weak_rev = false - testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) - - k = rand(rng, numparams(rn)) - rates = Dict(zip(parameters(rn), k)) - @test Catalyst.iscomplexbalanced(rn, rates) == false -end - -let - rn = @reaction_network begin - (k2, k1), A + E <--> AE - (k4, k3), AE <--> B + E - end - rev = true - weak_rev = true - testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) - - k = rand(rng, numparams(rn)) - rates = Dict(zip(parameters(rn), k)) - @test Catalyst.iscomplexbalanced(rn, rates) == true -end - -let - rn = @reaction_network begin (k2, k1), A + B <--> 2A end - rev = true - weak_rev = true - testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) - - k = rand(rng, numparams(rn)) - rates = Dict(zip(parameters(rn), k)) - @test Catalyst.iscomplexbalanced(rn, rates) == true -end - -let - rn = @reaction_network begin - k1, A + B --> 3A - k2, 3A --> 2A + C - k3, 2A + C --> 2B - k4, 2B --> A + B - end - rev = false - weak_rev = true - testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) - - k = rand(rng, numparams(rn)) - rates = Dict(zip(parameters(rn), k)) - @test Catalyst.iscomplexbalanced(rn, rates) == true -end - -let - rn = @reaction_network begin - (k2, k1), A + E <--> AE - (k4, k3), AE <--> B + E - k5, B --> 0 - k6, 0 --> A - end - rev = false - weak_rev = false - testreversibility(rn, reactioncomplexes(rn)[2], rev, weak_rev) - - k = rand(rng, numparams(rn)) - rates = Dict(zip(parameters(rn), k)) - @test Catalyst.iscomplexbalanced(rn, rates) == false -end - -let - rn = @reaction_network begin - k1, 3A + 2B --> 3C - k2, B + 4D --> 2E - k3, 2E --> 3C - (k4, k5), B + 4D <--> 3A + 2B - k6, F --> B + 4D - k7, 3C --> F - end - - 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 @@ -499,106 +273,6 @@ let @test isempty(cyclemat) 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)) - - 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 - (k1, k2), A <--> B - (k3, k4), B <--> C - (k5, k6), C <--> A - end - - 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. -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 - - rcs, D = reactioncomplexes(rn) - k = rand(rng, numparams(rn)) - p = parameters(rn) - rates = Dict(zip(parameters(rn), k)) - @test Catalyst.isdetailedbalanced(rn, rates) == false - - # 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 -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 - - 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 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 ### Other Network Properties Tests ### @@ -624,6 +298,7 @@ let complexoutgoingmat(rs; sparse = true) == sparse(cmplx_out_mat) end + # Tests outgoing complexes matrices (2). # Checks using dense and sparse representation. let @@ -649,286 +324,103 @@ let complexoutgoingmat(rs; sparse = true) == sparse(cmplx_out_mat) end -# Tests that `iscomplexbalanced` works for different rate inputs. -# Tests that non-valid rate input yields and error -let - # Declares network. - rn = @reaction_network begin - k1, 3A + 2B --> 3C - k2, B + 4D --> 2E - k3, 2E --> 3C - (k4, k5), B + 4D <--> 3A + 2B - k6, F --> B + 4D - k7, 3C --> F - end - - # Declares rate alternatives. - k = rand(rng, numparams(rn)) - rates_vec = Pair.(parameters(rn), k) - rates_tup = Tuple(rates_vec) - rates_dict = Dict(rates_vec) - rates_invalid = k - - # Tests that inputs are handled correctly. - @test Catalyst.iscomplexbalanced(rn, rates_vec) == Catalyst.iscomplexbalanced(rn, rates_tup) - @test Catalyst.iscomplexbalanced(rn, rates_tup) == Catalyst.iscomplexbalanced(rn, rates_dict) - @test_throws Exception Catalyst.iscomplexbalanced(rn, k) -end - -# Tests rate matrix computation for various input types. -let - # Declares network and its known rate matrix. - rn = @reaction_network begin - (k2, k1), A1 <--> A2 + A3 - k3, A2 + A3 --> A4 - k4, A4 --> A5 - (k6, k5), A5 <--> 2A6 - k7, 2A6 --> A4 - k8, A4 + A5 --> A7 - end - rate_mat = [ - 0.0 1.0 0.0 0.0 0.0 0.0 0.0; - 2.0 0.0 3.0 0.0 0.0 0.0 0.0; - 0.0 0.0 0.0 4.0 0.0 0.0 0.0; - 0.0 0.0 0.0 0.0 5.0 0.0 0.0; - 0.0 0.0 7.0 6.0 0.0 0.0 0.0; - 0.0 0.0 0.0 0.0 0.0 0.0 8.0; - 0.0 0.0 0.0 0.0 0.0 0.0 0.0; - ] - - # Declares rate alternatives. - rate_vals = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0] - rates_vec = Pair.(parameters(rn), rate_vals) - rates_tup = Tuple(rates_vec) - rates_dict = Dict(rates_vec) - rates_invalid = reshape(rate_vals, 1, 8) - - # Tests that all input types generates the correct rate matrix. - Catalyst.ratematrix(rn, rate_vals) == rate_mat - Catalyst.ratematrix(rn, rates_vec) == rate_mat - Catalyst.ratematrix(rn, rates_tup) == rate_mat - Catalyst.ratematrix(rn, rates_dict) == rate_mat - @test_throws Exception Catalyst.iscomplexbalanced(rn, rates_invalid) -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 - 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 - - -### 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)) - - 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 - (k1, k2), A <--> B - (k3, k4), B <--> C - (k5, k6), C <--> A - end - - 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. -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 - - rcs, D = reactioncomplexes(rn) - k = rand(rng, numparams(rn)) - p = parameters(rn) - rates = Dict(zip(parameters(rn), k)) - @test Catalyst.isdetailedbalanced(rn, rates) == false - - # 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 +### Tests for the matrices and vectors that are in the species-formation rate function +# ẋ = S * K * Φ(t) +# ẋ = Y * A_K * Φ(t) -# 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 - - 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 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 - -# 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 + MAPK = @reaction_network MAPK begin + (k₁, k₂),KKK + E1 <--> KKKE1 + k₃, KKKE1 --> KKK_ + E1 + (k₄, k₅), KKK_ + E2 <--> KKKE2 + k₆, KKKE2 --> KKK + E2 + (k₇, k₈), KK + KKK_ <--> KK_KKK_ + k₉, KK_KKK_ --> KKP + KKK_ + (k₁₀, k₁₁), KKP + KKK_ <--> KKPKKK_ + k₁₂, KKPKKK_ --> KKPP + KKK_ + (k₁₃, k₁₄), KKP + KKPase <--> KKPKKPase + k₁₅, KKPPKKPase --> KKP + KKPase + k₁₆,KKPKKPase --> KK + KKPase + (k₁₇, k₁₈), KKPP + KKPase <--> KKPPKKPase + (k₁₉, k₂₀), KKPP + K <--> KKPPK + k₂₁, KKPPK --> KKPP + KP + (k₂₂, k₂₃), KKPP + KP <--> KPKKPP + k₂₄, KPKKPP --> KPP + KKPP + (k₂₅, k₂₆), KP + KPase <--> KPKPase + k₂₇, KKPPKPase --> KP + KPase + k₂₈, KPKPase --> K + KPase + (k₂₉, k₃₀), KPP + KPase <--> KKPPKPase end - @test Catalyst.satisfiesdeficiencyone(rn) == false -end + Φ = Catalyst.massactionvector(MAPK) + specs = species(MAPK) + @test isequal(Φ[1], MAPK.KKK * MAPK.E1) + @test isequal(Φ[6], MAPK.KKK * MAPK.E2) + @test isequal(Φ[8], MAPK.KK_KKK_) + + K = fluxmat(MAPK) + @test isequal(K[1, 1], MAPK.k₁) + @test all(==(0), K[1, 2:end]) + @test isequal(K[2, 2], MAPK.k₂) + @test all(==(0), vcat(K[2,1], K[2,3:end])) + @test isequal(K[3, 2], MAPK.k₃) + @test all(==(0), vcat(K[3,1], K[3,3:end])) + @test count(k -> !isequal(k, 0), K) == length(reactions(MAPK)) + K = fluxmat(MAPK; sparse = true) + @test Catalyst.issparse(K) + + A_k = Catalyst.laplacianmat(MAPK) + @test all(col -> sum(col) == 0, eachcol(A_k)) + A_k = Catalyst.laplacianmat(MAPK; sparse = true) + @test Catalyst.issparse(A_k) -# Fails because linkage deficiencies do not sum to total deficiency -let - rn = @reaction_network begin - (k1, k2), A <--> 2A - (k3, k4), A + B <--> C - (k5, k6), C <--> B - end + S = netstoichmat(MAPK); Y = complexstoichmat(MAPK) + @test isequal(S*K, Y*A_k) - @test Catalyst.satisfiesdeficiencyone(rn) == false -end + eqs = Catalyst.assemble_oderhs(MAPK, specs) + @test all(iszero, simplify(eqs - S*K*Φ)) + @test all(iszero, simplify(eqs - Y*A_k*Φ)) -# Fails because a linkage class has deficiency two -let - rn = @reaction_network begin - k1, 3A --> A + 2B - k2, A + 2B --> 3B - k3, 3B --> 2A + B - k4, 2A + B --> 3A - end + # Test using numbers + k = rand(rng, numparams(MAPK)) + ratevec = collect(zip(parameters(MAPK), k)) + ratemap = Dict(ratemap) + ratetup = Tuple(ratemap) - @test Catalyst.satisfiesdeficiencyone(rn) == false -end + @test Catalyst.fluxmat(MAPK, ratemap) == Catalyst.fluxmat(MAPK, ratevec) == Catalyst.fluxmat(MAPK, ratetup) + K = Catalyst.fluxmat(MAPK, ratemap; sparse = true) + @test Catalyst.issparse(K) + + A_k = Catalyst.laplacianmat(MAPK, ratemap) + @test all(col -> sum(col) == 0, eachcol(A_k)) + A_k = Catalyst.laplacianmat(MAPK; sparse = true) + @test Catalyst.issparse(A_k) -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 + numeqs = similar(eqs) + for i in 1:length(eqs) + numeqs[i] = substitute(eqs[i], ratemap) end - - @test Catalyst.satisfiesdeficiencyone(rn) == true + @test all(iszero, simplify(numeqs - S*K*Φ)) + @test all(iszero, simplify(numeqs - Y*A_k*Φ)) end -### Some tests for deficiency zero networks. - +# Test handling for weird complexes. let rn = @reaction_network begin - (k1, k2), A <--> 2B - (k3, k4), A + C <--> D - k5, D --> B + E - k6, B + E --> A + C + k1, 2X + Y + 3Z --> ∅ end + Φ = Catalyst.massactionvector(rn) + @test isequal(Φ[1], rn.X^2 * rn.Y * rn.Z^3) + @test isequal(Φ[2], 1) - # 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 + u0vec = [:X => 3., :Y => .5, :Z => 2.] + u0map = Dict(u0vec) + u0tup = Tuple(u0vec) - # No longer weakly reversible - rn3 = @reaction_network begin - k1, A --> 2B - (k3, k4), A + C <--> 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 Catalyst.satisfiesdeficiencyzero(rn) == true - @test Catalyst.satisfiesdeficiencyzero(rn2) == false - @test Catalyst.satisfiesdeficiencyzero(rn3) == false - @test Catalyst.satisfiesdeficiencyzero(rn4) == false + Φ = Catalyst.massactionvector(rn, u0vec) + @test isequal(Φ[1], 36.) + Φ = Catalyst.massactionvector(rn, u0tup) + @test isequal(Φ[1], 36.) + Φ = Catalyst.massactionvector(rn, u0map) + @test isequal(Φ[1], 36.) end - diff --git a/test/runtests.jl b/test/runtests.jl index cd552c3157..216f7a8112 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -42,6 +42,7 @@ end # Tests reaction network analysis features. @time @safetestset "Conservation Laws" begin include("network_analysis/conservation_laws.jl") end @time @safetestset "Network Properties" begin include("network_analysis/network_properties.jl") end + @time @safetestset "CRN Theory" begin include("network_analysis/crn_theory.jl") end # Tests ODE, SDE, jump simulations, nonlinear solving, and steady state simulations. @time @safetestset "ODE System Simulations" begin include("simulation_and_solving/simulate_ODEs.jl") end From f8df438ac95fe77fe62c27d686242eca6bdb36bf Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 22 Nov 2024 22:35:07 -0500 Subject: [PATCH 04/29] exports --- src/Catalyst.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 73fbc7586f..1d809cca60 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -128,8 +128,8 @@ export @reaction_network, @network_component, @reaction, @species # Network analysis functionality. include("network_analysis.jl") export reactioncomplexmap, reactioncomplexes, incidencemat -export complexstoichmat -export complexoutgoingmat, incidencematgraph, linkageclasses, stronglinkageclasses, +export complexstoichmat, laplacianmat, fluxmat, massactionvector, complexoutgoingmat +export incidencematgraph, linkageclasses, stronglinkageclasses, terminallinkageclasses, deficiency, subnetworks export linkagedeficiencies, isreversible, isweaklyreversible export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants From 389d09a1170726e0b42540311fad6c41aee0f6e9 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 22 Nov 2024 22:38:53 -0500 Subject: [PATCH 05/29] cleanup --- test/network_analysis/network_properties.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 7ffbc166ed..a123ccf546 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -8,8 +8,6 @@ using Catalyst, LinearAlgebra, Test, SparseArrays using StableRNGs rng = StableRNG(514) -### Defining reaction networks. - ### Basic Tests ### # Tests basic `ReactionComplex` properties. @@ -62,7 +60,6 @@ let k = rand(rng, numparams(MAPK)) rates = Dict(zip(parameters(MAPK), k)) @test Catalyst.iscomplexbalanced(MAPK, rates) == false - cyclemat = Catalyst.cycles(MAPK) S = netstoichmat(MAPK) for i in 1:size(S, 2)-1 @@ -298,7 +295,6 @@ let complexoutgoingmat(rs; sparse = true) == sparse(cmplx_out_mat) end - # Tests outgoing complexes matrices (2). # Checks using dense and sparse representation. let From 7d6cbf58a7afc2167970e2dda3f3833e076a53ba Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 22 Nov 2024 22:51:08 -0500 Subject: [PATCH 06/29] fix test --- test/network_analysis/network_properties.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index a123ccf546..fe1676e119 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -380,8 +380,8 @@ let # Test using numbers k = rand(rng, numparams(MAPK)) ratevec = collect(zip(parameters(MAPK), k)) - ratemap = Dict(ratemap) - ratetup = Tuple(ratemap) + ratemap = Dict(ratevec) + ratetup = Tuple(ratevec) @test Catalyst.fluxmat(MAPK, ratemap) == Catalyst.fluxmat(MAPK, ratevec) == Catalyst.fluxmat(MAPK, ratetup) K = Catalyst.fluxmat(MAPK, ratemap; sparse = true) From 49225618bca6b1a5301eaeb29bee8271556362b3 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 23 Nov 2024 09:31:51 -0500 Subject: [PATCH 07/29] fix tesst --- 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 fe1676e119..5fa8a77c09 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -389,7 +389,7 @@ let A_k = Catalyst.laplacianmat(MAPK, ratemap) @test all(col -> sum(col) == 0, eachcol(A_k)) - A_k = Catalyst.laplacianmat(MAPK; sparse = true) + A_k = Catalyst.laplacianmat(MAPK, ratemap; sparse = true) @test Catalyst.issparse(A_k) numeqs = similar(eqs) From 536a50d5e33a87e1f06f1d3c24d1bf579da02dcf Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 23 Nov 2024 09:49:22 -0500 Subject: [PATCH 08/29] fix test --- test/network_analysis/crn_theory.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/network_analysis/crn_theory.jl b/test/network_analysis/crn_theory.jl index f2f40040c9..deb53b21b5 100644 --- a/test/network_analysis/crn_theory.jl +++ b/test/network_analysis/crn_theory.jl @@ -1,4 +1,6 @@ # Tests for properties from chemical reaction network theory: deficiency theorems, complex/detailed balance, etc. +using Catalyst, StableRNGs, LinearAlgebra, Test +rng = StableRNG(514) # Tests that `iscomplexbalanced` works for different rate inputs. # Tests that non-valid rate input yields and error From 835b4fff43dcff19872f760bb80535861baba361 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 23 Nov 2024 09:52:38 -0500 Subject: [PATCH 09/29] rm from NetworkProperties --- src/reactionsystem.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 6d24d28315..83ed069ea0 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -92,7 +92,6 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re incidencemat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0) complexstoichmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0) complexoutgoingmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0) - networklaplacian::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0) incidencegraph::Graphs.SimpleDiGraph{Int} = Graphs.DiGraph() linkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0) stronglinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0) From 4a691b4af5be34be54767f796e288164dd362f62 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 23 Nov 2024 10:38:55 -0500 Subject: [PATCH 10/29] comment out stability tests --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 216f7a8112..45357ff72e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -93,7 +93,7 @@ end # @time @safetestset "Structural Identifiability Extension" begin include("extensions/structural_identifiability.jl") end # Tests stability computation (but requires the HomotopyContinuation extension). - @time @safetestset "Steady State Stability Computations" begin include("extensions/stability_computation.jl") end + # @time @safetestset "Steady State Stability Computations" begin include("extensions/stability_computation.jl") end # Test spatial plotting, using CarioMakie and GraphMakie @time @safetestset "Lattice Simulation Plotting" begin include("extensions/lattice_simulation_plotting.jl") end From 211ac81db9d370d295327df05d8c47fa8b8ce89d Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 23 Nov 2024 11:22:53 -0500 Subject: [PATCH 11/29] update docstrings --- src/network_analysis.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 873922be66..bc3bfb014f 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -193,7 +193,7 @@ function complexstoichmat(::Type{Matrix{Int}}, rn::ReactionSystem, rcs) end """ - laplacianmat(rn::ReactionSystem, parammap=nothing; sparse=false) + laplacianmat(rn::ReactionSystem, pmap=Dict(); sparse=false) Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the graph Laplacian, and Φ is the [`massactionvector`](@ref). A_k is an n-by-n matrix, where n is the number of complexes, where A_{ij} = k_{ij} if a reaction exists between the two complexes and 0 otherwise. Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via parammap. @@ -204,11 +204,12 @@ function laplacianmat(rn::ReactionSystem, pmap = Dict(); sparse = false) end """ - fluxmat(rn::ReactionSystem, parammap = nothing; sparse=true) + fluxmat(rn::ReactionSystem, pmap = Dict(); sparse=false) Return an r×c matrix K such that, if complex j is the substrate complex of reaction i, then K_{ij} = k, the rate constant for this reaction. Mostly a helper function for the network Laplacian, [`networklaplacianmat`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix. + Returns a symbolic matrix by default, but will return a numerical matrix if rate constants are specified as a `Tuple`, `Vector`, or `Dict` of symbol-value pairs via `pmap`. """ -function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=true) +function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) rates = reactionrates(rn) if !isempty(pmap) @@ -265,7 +266,7 @@ function fluxmat(rn::ReactionSystem, pmap::Tuple; sparse = false) end """ - massactionvector(rn::ReactionSystem, scmap = nothing) + massactionvector(rn::ReactionSystem, scmap = Dict()) Return the vector whose entries correspond to the "mass action products" of each complex. For example, given the complex A + B, the corresponding entry of the vector would be A*B, and for the complex 2X + Y, the corresponding entry would be X^2*Y. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the [`networklaplacian`](@ref). This utility returns Φ(x). Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap. @@ -275,7 +276,7 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict()) sm = speciesmap(rn); specs = species(rn) if !all(r -> ismassaction(r, rn), rxs) - error("The supplied ReactionSystem has reactions that are not ismassaction. The monomial vector is only defined for pure mass action networks.") + error("The supplied ReactionSystem has reactions that are not ismassaction. The mass action vector is only defined for pure mass action networks.") end if !isempty(scmap) From 3bdf298f2ac9723d6b2142c077d7aa454064d57b Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 23 Nov 2024 11:24:39 -0500 Subject: [PATCH 12/29] up --- 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 bc3bfb014f..98211a358a 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -206,7 +206,7 @@ end """ fluxmat(rn::ReactionSystem, pmap = Dict(); sparse=false) - Return an r×c matrix K such that, if complex j is the substrate complex of reaction i, then K_{ij} = k, the rate constant for this reaction. Mostly a helper function for the network Laplacian, [`networklaplacianmat`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix. + Return an r×c matrix K such that, if complex j is the substrate complex of reaction i, then K_{ij} = k, the rate constant for this reaction. Mostly a helper function for the network Laplacian, [`networklaplacianmat`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix and Φ(x) is the [`massactionvector`](@ref). Returns a symbolic matrix by default, but will return a numerical matrix if rate constants are specified as a `Tuple`, `Vector`, or `Dict` of symbol-value pairs via `pmap`. """ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) From 312862885ebfe6079fd69bc53dcb647506b0c569 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 25 Nov 2024 17:57:00 -0500 Subject: [PATCH 13/29] up --- src/network_analysis.jl | 50 +++++----- test/network_analysis/network_properties.jl | 105 ++++++++++++++++++-- 2 files changed, 120 insertions(+), 35 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 98211a358a..f73838e047 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -195,8 +195,8 @@ end """ laplacianmat(rn::ReactionSystem, pmap=Dict(); sparse=false) - Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the graph Laplacian, and Φ is the [`massactionvector`](@ref). A_k is an n-by-n matrix, where n is the number of complexes, where A_{ij} = k_{ij} if a reaction exists between the two complexes and 0 otherwise. - Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via parammap. + Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the graph Laplacian, and Φ is the [`massactionvector`](@ref). A_k is an n-by-n matrix, where n is the number of complexes, where ``A_{ij} = k_{ij}`` if a reaction exists between the two complexes and 0 otherwise. + Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. """ function laplacianmat(rn::ReactionSystem, pmap = Dict(); sparse = false) D = incidencemat(rn; sparse); K = fluxmat(rn, pmap; sparse) @@ -212,19 +212,13 @@ end function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) rates = reactionrates(rn) - if !isempty(pmap) - length(pmap) != length(parameters(rn)) && error("Incorrect number of parameters were specified.") - pmap = symmap_to_varmap(rn, pmap) - pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap) - rates = [substitute(rate, pmap) for rate in rates] - end - + !isempty(pmap) && (rates = substitutevals(rn, pmap, parameters(rn), rates)) rcmap = reactioncomplexmap(rn); nc = length(rcmap); nr = length(rates) mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates) - K = if sparse - fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates) + if sparse + return fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates) else - fluxmat(Matrix{mtype}, rcmap, rates) + return fluxmat(Matrix{mtype}, rcmap, rates) end end @@ -265,13 +259,22 @@ function fluxmat(rn::ReactionSystem, pmap::Tuple; sparse = false) fluxmat(rn, pdict; sparse) end +# Helper to substitute values into a (vector of) symbolic expressions. The syms are the symbols to substitute and the symexprs are the expressions to substitute into. +function substitutevals(rn::ReactionSystem, map::Dict, syms, symexprs) + length(map) != length(syms) && error("Incorrect number of parameter-value pairs were specified.") + map = symmap_to_varmap(rn, map) + map = Dict(ModelingToolkit.value(k) => v for (k, v) in map) + vals = [substitute(expr, map) for expr in symexprs] +end + """ - massactionvector(rn::ReactionSystem, scmap = Dict()) + massactionvector(rn::ReactionSystem, scmap = Dict(); combinatoric_ratelaws = true) Return the vector whose entries correspond to the "mass action products" of each complex. For example, given the complex A + B, the corresponding entry of the vector would be A*B, and for the complex 2X + Y, the corresponding entry would be X^2*Y. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the [`networklaplacian`](@ref). This utility returns Φ(x). Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap. + If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws). """ -function massactionvector(rn::ReactionSystem, scmap::Dict = Dict()) +function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = true) r = numreactions(rn); rxs = reactions(rn) sm = speciesmap(rn); specs = species(rn) @@ -279,12 +282,7 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict()) error("The supplied ReactionSystem has reactions that are not ismassaction. The mass action vector is only defined for pure mass action networks.") end - if !isempty(scmap) - length(scmap) != length(specs) && error("Incorrect number of species concentrations were specified.") - scmap = symmap_to_varmap(rn, scmap) - scmap = Dict(ModelingToolkit.value(k) => v for (k, v) in scmap) - specs = [substitute(s, scmap) for s in specs] - end + !isempty(scmap) && (specs = substitutevals(rn, scmap, specs, specs)) vtype = eltype(specs) <: Symbolics.BasicSymbolic ? Num : eltype(specs) Φ = Vector{vtype}() @@ -292,20 +290,22 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict()) for comp in keys(reactioncomplexmap(rn)) subs = map(ce -> getfield(ce, :speciesid), comp) stoich = map(ce -> getfield(ce, :speciesstoich), comp) - push!(Φ, prod(vtype[specs[s]^α for (s, α) in zip(subs, stoich)])) + maprod = prod(vtype[specs[s]^α for (s, α) in zip(subs, stoich)]) + combinatoric_ratelaws && (maprod /= prod(map(factorial, stoich))) + push!(Φ, maprod) end Φ end -function massactionvector(rn::ReactionSystem, scmap::Tuple) +function massactionvector(rn::ReactionSystem, scmap::Tuple; combinatoric_ratelaws = true) sdict = Dict(scmap) - massactionvector(rn, sdict) + massactionvector(rn, sdict; combinatoric_ratelaws) end -function massactionvector(rn::ReactionSystem, scmap::Vector) +function massactionvector(rn::ReactionSystem, scmap::Vector; combinatoric_ratelaws = true) sdict = Dict(scmap) - massactionvector(rn, sdict) + massactionvector(rn, sdict; combinatoric_ratelaws) end @doc raw""" diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 5fa8a77c09..4a70b6e970 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -350,11 +350,44 @@ let Φ = Catalyst.massactionvector(MAPK) specs = species(MAPK) - @test isequal(Φ[1], MAPK.KKK * MAPK.E1) - @test isequal(Φ[6], MAPK.KKK * MAPK.E2) - @test isequal(Φ[8], MAPK.KK_KKK_) + truevec = [MAPK.KKK * MAPK.E1, + MAPK.KKKE1, + MAPK.KKK_ * MAPK.E1, + MAPK.KKK_ * MAPK.E2, + MAPK.KKKE2, + MAPK.KKK * MAPK.E2, + MAPK.KK * MAPK.KKK_, + MAPK.KK_KKK_, + MAPK.KKP * MAPK.KKK_, + MAPK.KKPKKK_, + MAPK.KKPP * MAPK.KKK_, + MAPK.KKP * MAPK.KKPase, + MAPK.KKPKKPase, + MAPK.KKPPKKPase, + MAPK.KK * MAPK.KKPase, + MAPK.KKPP * MAPK.KKPase, + MAPK.KKPP * MAPK.K, + MAPK.KKPPK, + MAPK.KKPP * MAPK.KP, + MAPK.KPKKPP, + MAPK.KPP * MAPK.KKPP, + MAPK.KP * MAPK.KPase, + MAPK.KPKPase, + MAPK.KKPPKPase, + MAPK.K * MAPK.KPase, + MAPK.KPP * MAPK.KPase, + ] + @test isequal(Φ, truevec) K = fluxmat(MAPK) + # Construct matrix from incidence matrix + mat = zeros(Num, 30, 26); D = incidencemat(MAPK) + rates = reactionrates(MAPK) + for (i, col) in enumerate(eachcol(D)) + sub = findfirst(==(-1), col) + mat[i, sub] = rates[i] + end + @test isequal(K, mat) @test isequal(K[1, 1], MAPK.k₁) @test all(==(0), K[1, 2:end]) @test isequal(K[2, 2], MAPK.k₂) @@ -400,23 +433,75 @@ let @test all(iszero, simplify(numeqs - Y*A_k*Φ)) end -# Test handling for weird complexes. +# Test handling for weird complexes and combinatoric rate laws. let rn = @reaction_network begin k1, 2X + Y + 3Z --> ∅ + (k2, k3), 2Y + 2Z <--> 3X + end + + Φ = Catalyst.massactionvector(rn); specs = species(rn) + crvec = [rn.X^2/2 * rn.Y * rn.Z^3/6, + 1., + rn.Y^2/2 * rn.Z^2/2, + rn.X^3/6] + @test isequal(Φ, crvec) + ncrvec = [rn.X^2 * rn.Y * rn.Z^3, + 1., + rn.Y^2 * rn.Z^2, + rn.X^3] + Φ_2 = Catalyst.massactionvector(rn; combinatoric_ratelaws = false) + @test isequal(Φ_2, ncrvec) + + # Test that the ODEs generated are the same. + eqs = Catalyst.assemble_oderhs(rn, specs) + S = netstoichmat(rn); Y = complexstoichmat(rn); K = fluxmat(rn); A_k = laplacianmat(rn) + @test all(iszero, simplify(eqs - S*K*Φ)) + @test all(iszero, simplify(eqs - Y*A_k*Φ)) + + eq_ncr = Catalyst.assemble_oderhs(rn, specs; combinatoric_ratelaws = false) + @test all(iszero, simplify(eq_ncr - S*K*Φ_2)) + @test all(iszero, simplify(eq_ncr - Y*A_k*Φ_2)) + + # Test that the ODEs with rate constants are the same. + k = rand(rng, numparams(rn)) + ratevec = collect(zip(parameters(rn), k)) + ratemap = Dict(ratevec) + K = fluxmat(rn, ratemap); A_k = laplacianmat(rn, ratemap) + + numeqs = similar(eqs) + for i in 1:length(eqs) + numeqs[i] = substitute(eqs[i], ratemap) + end + @test_broken all(iszero, simplify(numeqs - S*K*Φ)) + @test_broken all(iszero, simplify(numeqs - Y*A_k*Φ)) + + numeqs_ncr = similar(eq_ncr) + for i in 1:length(eq_ncr) + numeqs_ncr[i] = substitute(eq_ncr[i], ratemap) end - Φ = Catalyst.massactionvector(rn) - @test isequal(Φ[1], rn.X^2 * rn.Y * rn.Z^3) - @test isequal(Φ[2], 1) + @test all(iszero, simplify(numeqs_ncr - S*K*Φ_2)) + @test all(iszero, simplify(numeqs_ncr - Y*A_k*Φ_2)) + # Test that handling of species concentrations is correct. u0vec = [:X => 3., :Y => .5, :Z => 2.] u0map = Dict(u0vec) u0tup = Tuple(u0vec) Φ = Catalyst.massactionvector(rn, u0vec) - @test isequal(Φ[1], 36.) - Φ = Catalyst.massactionvector(rn, u0tup) + @test isequal(Φ[1], 3.) + Φ = Catalyst.massactionvector(rn, u0tup; combinatoric_ratelaws = false) @test isequal(Φ[1], 36.) Φ = Catalyst.massactionvector(rn, u0map) - @test isequal(Φ[1], 36.) + @test isequal(Φ[1], 3.) + + # Test full simplification. + u0map = symmap_to_varmap(rn, u0map) + numeqs = [substitute(eq, u0map) for eq in numeqs] + @test isapprox(numeqs, S*K*Φ) + @test isapprox(numeqs, Y*A_k*Φ) + + numeqs_ncr = [substitute(eq, u0map) for eq in numeqs_ncr] + @test isapprox(numeqs_ncr - S*K*Φ_2) + @test isapprox(numeqs_ncr - Y*A_k*Φ_2) end From b921fed28df3526a6013d3bb14057e93bad643e9 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 25 Nov 2024 17:59:49 -0500 Subject: [PATCH 14/29] Fix docstring --- 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 f73838e047..5b8a92db60 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -192,10 +192,10 @@ function complexstoichmat(::Type{Matrix{Int}}, rn::ReactionSystem, rcs) Z end -""" +@doc raw""" laplacianmat(rn::ReactionSystem, pmap=Dict(); sparse=false) - Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the graph Laplacian, and Φ is the [`massactionvector`](@ref). A_k is an n-by-n matrix, where n is the number of complexes, where ``A_{ij} = k_{ij}`` if a reaction exists between the two complexes and 0 otherwise. + Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the graph Laplacian, and ``Φ`` is the [`massactionvector`](@ref). ``A_k`` is an n-by-n matrix, where n is the number of complexes, where ``A_{ij} = k_{ij}`` if a reaction exists between the two complexes and 0 otherwise. Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. """ function laplacianmat(rn::ReactionSystem, pmap = Dict(); sparse = false) @@ -203,10 +203,10 @@ function laplacianmat(rn::ReactionSystem, pmap = Dict(); sparse = false) D*K end -""" +@doc raw""" fluxmat(rn::ReactionSystem, pmap = Dict(); sparse=false) - Return an r×c matrix K such that, if complex j is the substrate complex of reaction i, then K_{ij} = k, the rate constant for this reaction. Mostly a helper function for the network Laplacian, [`networklaplacianmat`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix and Φ(x) is the [`massactionvector`](@ref). + Return an r×c matrix ``K`` such that, if complex ``j`` is the substrate complex of reaction ``i``, then ``K_{ij} = k``, the rate constant for this reaction. Mostly a helper function for the network Laplacian, [`laplacianmat`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix and ``Φ(x)`` is the [`massactionvector`](@ref). Returns a symbolic matrix by default, but will return a numerical matrix if rate constants are specified as a `Tuple`, `Vector`, or `Dict` of symbol-value pairs via `pmap`. """ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) @@ -270,7 +270,7 @@ end """ massactionvector(rn::ReactionSystem, scmap = Dict(); combinatoric_ratelaws = true) - Return the vector whose entries correspond to the "mass action products" of each complex. For example, given the complex A + B, the corresponding entry of the vector would be A*B, and for the complex 2X + Y, the corresponding entry would be X^2*Y. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where Y is the [`complexstoichmat`](@ref) and A_k is the negative of the [`networklaplacian`](@ref). This utility returns Φ(x). + Return the vector whose entries correspond to the "mass action products" of each complex. For example, given the complex A + B, the corresponding entry of the vector would be ``A*B``, and for the complex 2X + Y, the corresponding entry would be ``X^2*Y``. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the [`laplacianmat`](@ref). This utility returns ``Φ(x)``. Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap. If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws). """ From c82c80e26d273c076d78246be6db17e46e516778 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 25 Nov 2024 18:00:46 -0500 Subject: [PATCH 15/29] add comment --- test/network_analysis/network_properties.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 4a70b6e970..e963630f9f 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -473,6 +473,7 @@ let for i in 1:length(eqs) numeqs[i] = substitute(eqs[i], ratemap) end + # Broken but the difference is just numerical, something on the order of 1e-17 times a term @test_broken all(iszero, simplify(numeqs - S*K*Φ)) @test_broken all(iszero, simplify(numeqs - Y*A_k*Φ)) From cf89a29ee670b8119667fba2d9372890c28ae822 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 25 Nov 2024 18:22:29 -0500 Subject: [PATCH 16/29] test fix --- 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 e963630f9f..ca91e175b5 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -474,8 +474,8 @@ let numeqs[i] = substitute(eqs[i], ratemap) end # Broken but the difference is just numerical, something on the order of 1e-17 times a term - @test_broken all(iszero, simplify(numeqs - S*K*Φ)) - @test_broken all(iszero, simplify(numeqs - Y*A_k*Φ)) + @test all(iszero, simplify(numeqs - S*K*Φ)) + @test all(iszero, simplify(numeqs - Y*A_k*Φ)) numeqs_ncr = similar(eq_ncr) for i in 1:length(eq_ncr) @@ -503,6 +503,6 @@ let @test isapprox(numeqs, Y*A_k*Φ) numeqs_ncr = [substitute(eq, u0map) for eq in numeqs_ncr] - @test isapprox(numeqs_ncr - S*K*Φ_2) - @test isapprox(numeqs_ncr - Y*A_k*Φ_2) + @test isapprox(numeqs_ncr, S*K*Φ_2) + @test isapprox(numeqs_ncr, Y*A_k*Φ_2) end From 6e6843658aa6bab772106d8c5fafd8af2af6ff8e Mon Sep 17 00:00:00 2001 From: vyudu Date: Tue, 26 Nov 2024 00:04:11 -0500 Subject: [PATCH 17/29] fix test --- test/network_analysis/network_properties.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index ca91e175b5..988c5e2365 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -491,8 +491,8 @@ let Φ = Catalyst.massactionvector(rn, u0vec) @test isequal(Φ[1], 3.) - Φ = Catalyst.massactionvector(rn, u0tup; combinatoric_ratelaws = false) - @test isequal(Φ[1], 36.) + Φ_2 = Catalyst.massactionvector(rn, u0tup; combinatoric_ratelaws = false) + @test isequal(Φ_2[1], 36.) Φ = Catalyst.massactionvector(rn, u0map) @test isequal(Φ[1], 3.) From b3dce887962ffc36da6b5d3b4c8b2006ee5f760f Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 4 Dec 2024 14:17:50 -0500 Subject: [PATCH 18/29] format --- src/network_analysis.jl | 26 +++++++++++++-------- test/network_analysis/network_properties.jl | 17 ++++++++++---- 2 files changed, 28 insertions(+), 15 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 5b8a92db60..7b3f817416 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -199,7 +199,8 @@ end Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. """ function laplacianmat(rn::ReactionSystem, pmap = Dict(); sparse = false) - D = incidencemat(rn; sparse); K = fluxmat(rn, pmap; sparse) + D = incidencemat(rn; sparse) + K = fluxmat(rn, pmap; sparse) D*K end @@ -213,7 +214,9 @@ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) rates = reactionrates(rn) !isempty(pmap) && (rates = substitutevals(rn, pmap, parameters(rn), rates)) - rcmap = reactioncomplexmap(rn); nc = length(rcmap); nr = length(rates) + rcmap = reactioncomplexmap(rn) + nc = length(rcmap) + nr = length(rates) mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates) if sparse return fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates) @@ -239,7 +242,8 @@ function fluxmat(::Type{SparseMatrixCSC{T, Int}}, rcmap, rates) where T end function fluxmat(::Type{Matrix{T}}, rcmap, rates) where T - nr = length(rates); nc = length(rcmap) + nr = length(rates) + nc = length(rcmap) K = zeros(T, nr, nc) for (i, (complex, rxs)) in enumerate(rcmap) for (rx, dir) in rxs @@ -260,7 +264,7 @@ function fluxmat(rn::ReactionSystem, pmap::Tuple; sparse = false) end # Helper to substitute values into a (vector of) symbolic expressions. The syms are the symbols to substitute and the symexprs are the expressions to substitute into. -function substitutevals(rn::ReactionSystem, map::Dict, syms, symexprs) +function substitutevals(rn::ReactionSystem, map::Dict, syms, symexprs) length(map) != length(syms) && error("Incorrect number of parameter-value pairs were specified.") map = symmap_to_varmap(rn, map) map = Dict(ModelingToolkit.value(k) => v for (k, v) in map) @@ -272,11 +276,13 @@ end Return the vector whose entries correspond to the "mass action products" of each complex. For example, given the complex A + B, the corresponding entry of the vector would be ``A*B``, and for the complex 2X + Y, the corresponding entry would be ``X^2*Y``. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the [`laplacianmat`](@ref). This utility returns ``Φ(x)``. Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap. - If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws). + If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws). Will default to the default for the system. """ -function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = true) - r = numreactions(rn); rxs = reactions(rn) - sm = speciesmap(rn); specs = species(rn) +function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = rn.combinatoric_ratelaws) + r = numreactions(rn) + rxs = reactions(rn) + sm = speciesmap(rn) + specs = species(rn) if !all(r -> ismassaction(r, rn), rxs) error("The supplied ReactionSystem has reactions that are not ismassaction. The mass action vector is only defined for pure mass action networks.") @@ -298,12 +304,12 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric Φ end -function massactionvector(rn::ReactionSystem, scmap::Tuple; combinatoric_ratelaws = true) +function massactionvector(rn::ReactionSystem, scmap::Tuple; combinatoric_ratelaws = rn.combinatoric_ratelaws) sdict = Dict(scmap) massactionvector(rn, sdict; combinatoric_ratelaws) end -function massactionvector(rn::ReactionSystem, scmap::Vector; combinatoric_ratelaws = true) +function massactionvector(rn::ReactionSystem, scmap::Vector; combinatoric_ratelaws = rn.combinatoric_ratelaws) sdict = Dict(scmap) massactionvector(rn, sdict; combinatoric_ratelaws) end diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 988c5e2365..2fd7ff0f35 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -381,7 +381,8 @@ let K = fluxmat(MAPK) # Construct matrix from incidence matrix - mat = zeros(Num, 30, 26); D = incidencemat(MAPK) + mat = zeros(Num, 30, 26) + D = incidencemat(MAPK) rates = reactionrates(MAPK) for (i, col) in enumerate(eachcol(D)) sub = findfirst(==(-1), col) @@ -403,7 +404,8 @@ let A_k = Catalyst.laplacianmat(MAPK; sparse = true) @test Catalyst.issparse(A_k) - S = netstoichmat(MAPK); Y = complexstoichmat(MAPK) + S = netstoichmat(MAPK) + Y = complexstoichmat(MAPK) @test isequal(S*K, Y*A_k) eqs = Catalyst.assemble_oderhs(MAPK, specs) @@ -440,7 +442,8 @@ let (k2, k3), 2Y + 2Z <--> 3X end - Φ = Catalyst.massactionvector(rn); specs = species(rn) + Φ = Catalyst.massactionvector(rn) + specs = species(rn) crvec = [rn.X^2/2 * rn.Y * rn.Z^3/6, 1., rn.Y^2/2 * rn.Z^2/2, @@ -455,7 +458,10 @@ let # Test that the ODEs generated are the same. eqs = Catalyst.assemble_oderhs(rn, specs) - S = netstoichmat(rn); Y = complexstoichmat(rn); K = fluxmat(rn); A_k = laplacianmat(rn) + S = netstoichmat(rn) + Y = complexstoichmat(rn) + K = fluxmat(rn) + A_k = laplacianmat(rn) @test all(iszero, simplify(eqs - S*K*Φ)) @test all(iszero, simplify(eqs - Y*A_k*Φ)) @@ -467,7 +473,8 @@ let k = rand(rng, numparams(rn)) ratevec = collect(zip(parameters(rn), k)) ratemap = Dict(ratevec) - K = fluxmat(rn, ratemap); A_k = laplacianmat(rn, ratemap) + K = fluxmat(rn, ratemap) + A_k = laplacianmat(rn, ratemap) numeqs = similar(eqs) for i in 1:length(eqs) From b777518c0301bd7e46e3045932b1614e12c27bd8 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 4 Dec 2024 14:52:52 -0500 Subject: [PATCH 19/29] up --- src/network_analysis.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 7b3f817416..dd2d762182 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -211,9 +211,12 @@ end Returns a symbolic matrix by default, but will return a numerical matrix if rate constants are specified as a `Tuple`, `Vector`, or `Dict` of symbol-value pairs via `pmap`. """ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) - rates = reactionrates(rn) + rates = if isempty(pmap) + reactionrates(rn) + else + substitutevals(rn, pmap, parameters(rn), reactionrates(rn)) + end - !isempty(pmap) && (rates = substitutevals(rn, pmap, parameters(rn), rates)) rcmap = reactioncomplexmap(rn) nc = length(rcmap) nr = length(rates) @@ -282,14 +285,16 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric r = numreactions(rn) rxs = reactions(rn) sm = speciesmap(rn) - specs = species(rn) + specs = if isempty(scmap) + species(rn) + else + substitutevals(rn, scmap, species(rn), species(rn)) + end if !all(r -> ismassaction(r, rn), rxs) error("The supplied ReactionSystem has reactions that are not ismassaction. The mass action vector is only defined for pure mass action networks.") end - !isempty(scmap) && (specs = substitutevals(rn, scmap, specs, specs)) - vtype = eltype(specs) <: Symbolics.BasicSymbolic ? Num : eltype(specs) Φ = Vector{vtype}() rcmap = reactioncomplexmap(rn) From a782851d618f753178a62c3485fa8394ac4a6ef0 Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 4 Dec 2024 15:37:08 -0500 Subject: [PATCH 20/29] fix --- src/network_analysis.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index dd2d762182..36d4d2ea95 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -281,10 +281,11 @@ end Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap. If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws). Will default to the default for the system. """ -function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = rn.combinatoric_ratelaws) +function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) r = numreactions(rn) rxs = reactions(rn) sm = speciesmap(rn) + specs = if isempty(scmap) species(rn) else @@ -309,12 +310,12 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric Φ end -function massactionvector(rn::ReactionSystem, scmap::Tuple; combinatoric_ratelaws = rn.combinatoric_ratelaws) +function massactionvector(rn::ReactionSystem, scmap::Tuple; combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) sdict = Dict(scmap) massactionvector(rn, sdict; combinatoric_ratelaws) end -function massactionvector(rn::ReactionSystem, scmap::Vector; combinatoric_ratelaws = rn.combinatoric_ratelaws) +function massactionvector(rn::ReactionSystem, scmap::Vector; combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) sdict = Dict(scmap) massactionvector(rn, sdict; combinatoric_ratelaws) end From 32deeea148dcf55e7354e592b598f1035faebff7 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 13 Dec 2024 09:16:55 +0900 Subject: [PATCH 21/29] 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 36d4d2ea95..be1c6ef41b 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -220,7 +220,7 @@ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) rcmap = reactioncomplexmap(rn) nc = length(rcmap) nr = length(rates) - mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates) + mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Any : eltype(rates) if sparse return fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates) else @@ -296,7 +296,7 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric error("The supplied ReactionSystem has reactions that are not ismassaction. The mass action vector is only defined for pure mass action networks.") end - vtype = eltype(specs) <: Symbolics.BasicSymbolic ? Num : eltype(specs) + vtype = eltype(specs) <: Symbolics.BasicSymbolic ? Any : eltype(specs) Φ = Vector{vtype}() rcmap = reactioncomplexmap(rn) for comp in keys(reactioncomplexmap(rn)) From 3f1d1ca8556faf0d825fb55fba13a6f356f4fe2b Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 13 Dec 2024 09:51:36 +0900 Subject: [PATCH 22/29] up --- src/network_analysis.jl | 13 +++++++------ test/runtests.jl | 2 +- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index be1c6ef41b..6874969ac1 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -220,12 +220,13 @@ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) rcmap = reactioncomplexmap(rn) nc = length(rcmap) nr = length(rates) - mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Any : eltype(rates) - if sparse - return fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates) + mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates) + fluxmat = if sparse + fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates) else - return fluxmat(Matrix{mtype}, rcmap, rates) + fluxmat(Matrix{mtype}, rcmap, rates) end + mtype == Num ? Matrix{Any}(fluxmat) : fluxmat end function fluxmat(::Type{SparseMatrixCSC{T, Int}}, rcmap, rates) where T @@ -296,7 +297,7 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric error("The supplied ReactionSystem has reactions that are not ismassaction. The mass action vector is only defined for pure mass action networks.") end - vtype = eltype(specs) <: Symbolics.BasicSymbolic ? Any : eltype(specs) + vtype = eltype(specs) <: Symbolics.BasicSymbolic ? Num : eltype(specs) Φ = Vector{vtype}() rcmap = reactioncomplexmap(rn) for comp in keys(reactioncomplexmap(rn)) @@ -307,7 +308,7 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric push!(Φ, maprod) end - Φ + vtype == Num ? Vector{Any}(Φ) : Φ end function massactionvector(rn::ReactionSystem, scmap::Tuple; combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) diff --git a/test/runtests.jl b/test/runtests.jl index f251cd36d3..55e1d5d7c8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -49,7 +49,7 @@ end @time @safetestset "Automatic Jacobian Construction" begin include("simulation_and_solving/jacobian_construction.jl") end @time @safetestset "SDE System Simulations" begin include("simulation_and_solving/simulate_SDEs.jl") end @time @safetestset "Jump System Simulations" begin include("simulation_and_solving/simulate_jumps.jl") end - @time @safetestset "Nonlinear and SteadyState System Solving" begin include("simulation_and_solving/solve_nonlinear.jl") end + # @time @safetestset "Nonlinear and SteadyState System Solving" begin include("simulation_and_solving/solve_nonlinear.jl") end # Tests upstream SciML and DiffEq stuff. @time @safetestset "MTK Structure Indexing" begin include("upstream/mtk_structure_indexing.jl") end From ce8e1109ddbc08d7a79c65d5600ca6638f72033a Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 13 Dec 2024 10:15:09 +0900 Subject: [PATCH 23/29] fixes --- src/network_analysis.jl | 7 +++---- test/network_analysis/network_properties.jl | 4 ++-- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 6874969ac1..5e3373eb8e 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -221,12 +221,12 @@ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) nc = length(rcmap) nr = length(rates) mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates) - fluxmat = if sparse + mat = if sparse fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates) else fluxmat(Matrix{mtype}, rcmap, rates) end - mtype == Num ? Matrix{Any}(fluxmat) : fluxmat + mtype == Num ? Matrix{Any}(mat) : mat end function fluxmat(::Type{SparseMatrixCSC{T, Int}}, rcmap, rates) where T @@ -250,8 +250,7 @@ function fluxmat(::Type{Matrix{T}}, rcmap, rates) where T nc = length(rcmap) K = zeros(T, nr, nc) for (i, (complex, rxs)) in enumerate(rcmap) - for (rx, dir) in rxs - dir == -1 && (K[rx, i] = rates[rx]) + for (rx, dir) in rxs dir == -1 && (K[rx, i] = rates[rx]) end end K diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 2fd7ff0f35..4844007938 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -379,7 +379,7 @@ let ] @test isequal(Φ, truevec) - K = fluxmat(MAPK) + K = Catalyst.fluxmat(MAPK) # Construct matrix from incidence matrix mat = zeros(Num, 30, 26) D = incidencemat(MAPK) @@ -396,7 +396,7 @@ let @test isequal(K[3, 2], MAPK.k₃) @test all(==(0), vcat(K[3,1], K[3,3:end])) @test count(k -> !isequal(k, 0), K) == length(reactions(MAPK)) - K = fluxmat(MAPK; sparse = true) + K = Catalyst.fluxmat(MAPK; sparse = true) @test Catalyst.issparse(K) A_k = Catalyst.laplacianmat(MAPK) From 855145415dabdcad6c1116a56ac12dcc0b54e344 Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 14 Dec 2024 16:18:38 +0800 Subject: [PATCH 24/29] up --- src/network_analysis.jl | 67 ++++++++++----------- test/network_analysis/network_properties.jl | 13 +--- test/runtests.jl | 2 +- 3 files changed, 37 insertions(+), 45 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 5e3373eb8e..841551125a 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -198,19 +198,22 @@ end Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the graph Laplacian, and ``Φ`` is the [`massactionvector`](@ref). ``A_k`` is an n-by-n matrix, where n is the number of complexes, where ``A_{ij} = k_{ij}`` if a reaction exists between the two complexes and 0 otherwise. Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. """ -function laplacianmat(rn::ReactionSystem, pmap = Dict(); sparse = false) - D = incidencemat(rn; sparse) - K = fluxmat(rn, pmap; sparse) - D*K +function laplacianmat(rn::ReactionSystem, pmap::Dict = Dict()) + D = incidencemat(rn) + K = fluxmat(rn, pmap) + D * K end +Base.zero(::Type{Union{R, Symbolics.BasicSymbolic{Real}}}) where R <: Real = zero(R) +Base.one(::Type{Union{R, Symbolics.BasicSymbolic{Real}}}) where R <: Real = one(R) + @doc raw""" fluxmat(rn::ReactionSystem, pmap = Dict(); sparse=false) Return an r×c matrix ``K`` such that, if complex ``j`` is the substrate complex of reaction ``i``, then ``K_{ij} = k``, the rate constant for this reaction. Mostly a helper function for the network Laplacian, [`laplacianmat`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix and ``Φ(x)`` is the [`massactionvector`](@ref). Returns a symbolic matrix by default, but will return a numerical matrix if rate constants are specified as a `Tuple`, `Vector`, or `Dict` of symbol-value pairs via `pmap`. """ -function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) +function fluxmat(rn::ReactionSystem, pmap::Dict = Dict()) rates = if isempty(pmap) reactionrates(rn) else @@ -221,49 +224,45 @@ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) nc = length(rcmap) nr = length(rates) mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates) - mat = if sparse - fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates) - else - fluxmat(Matrix{mtype}, rcmap, rates) - end - mtype == Num ? Matrix{Any}(mat) : mat -end - -function fluxmat(::Type{SparseMatrixCSC{T, Int}}, rcmap, rates) where T - Is = Int[] - Js = Int[] - Vs = T[] - for (i, (complex, rxs)) in enumerate(rcmap) - for (rx, dir) in rxs - dir == -1 && begin - push!(Is, rx) - push!(Js, i) - push!(Vs, rates[rx]) - end - end - end - Z = sparse(Is, Js, Vs, length(rates), length(rcmap)) -end + Matrix{Any}(fluxmat(Matrix{mtype}, rcmap, rates)) +end + +# function fluxmat(::Type{SparseMatrixCSC{T, Int}}, rcmap, rates) where T +# Is = Int[] +# Js = Int[] +# Vs = T[] +# for (i, (complex, rxs)) in enumerate(rcmap) +# for (rx, dir) in rxs +# dir == -1 && begin +# push!(Is, rx) +# push!(Js, i) +# push!(Vs, rates[rx]) +# end +# end +# end +# Z = sparse(Is, Js, Vs, length(rates), length(rcmap)) +# end function fluxmat(::Type{Matrix{T}}, rcmap, rates) where T nr = length(rates) nc = length(rcmap) K = zeros(T, nr, nc) for (i, (complex, rxs)) in enumerate(rcmap) - for (rx, dir) in rxs dir == -1 && (K[rx, i] = rates[rx]) + for (rx, dir) in rxs + dir == -1 && (K[rx, i] = rates[rx]) end end K end -function fluxmat(rn::ReactionSystem, pmap::Vector; sparse = false) +function fluxmat(rn::ReactionSystem, pmap::Vector) pdict = Dict(pmap) - fluxmat(rn, pdict; sparse) + fluxmat(rn, pdict) end -function fluxmat(rn::ReactionSystem, pmap::Tuple; sparse = false) +function fluxmat(rn::ReactionSystem, pmap::Tuple) pdict = Dict(pmap) - fluxmat(rn, pdict; sparse) + fluxmat(rn, pdict) end # Helper to substitute values into a (vector of) symbolic expressions. The syms are the symbols to substitute and the symexprs are the expressions to substitute into. @@ -307,7 +306,7 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric push!(Φ, maprod) end - vtype == Num ? Vector{Any}(Φ) : Φ + Vector{Any}(Φ) end function massactionvector(rn::ReactionSystem, scmap::Tuple; combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl index 4844007938..0fa8808f7a 100644 --- a/test/network_analysis/network_properties.jl +++ b/test/network_analysis/network_properties.jl @@ -380,8 +380,8 @@ let @test isequal(Φ, truevec) K = Catalyst.fluxmat(MAPK) - # Construct matrix from incidence matrix - mat = zeros(Num, 30, 26) + # Construct flux matrix from incidence matrix + mat = Matrix{Any}(zeros(30, 26)) D = incidencemat(MAPK) rates = reactionrates(MAPK) for (i, col) in enumerate(eachcol(D)) @@ -396,13 +396,9 @@ let @test isequal(K[3, 2], MAPK.k₃) @test all(==(0), vcat(K[3,1], K[3,3:end])) @test count(k -> !isequal(k, 0), K) == length(reactions(MAPK)) - K = Catalyst.fluxmat(MAPK; sparse = true) - @test Catalyst.issparse(K) A_k = Catalyst.laplacianmat(MAPK) @test all(col -> sum(col) == 0, eachcol(A_k)) - A_k = Catalyst.laplacianmat(MAPK; sparse = true) - @test Catalyst.issparse(A_k) S = netstoichmat(MAPK) Y = complexstoichmat(MAPK) @@ -419,13 +415,10 @@ let ratetup = Tuple(ratevec) @test Catalyst.fluxmat(MAPK, ratemap) == Catalyst.fluxmat(MAPK, ratevec) == Catalyst.fluxmat(MAPK, ratetup) - K = Catalyst.fluxmat(MAPK, ratemap; sparse = true) - @test Catalyst.issparse(K) + K = Catalyst.fluxmat(MAPK, ratemap) A_k = Catalyst.laplacianmat(MAPK, ratemap) @test all(col -> sum(col) == 0, eachcol(A_k)) - A_k = Catalyst.laplacianmat(MAPK, ratemap; sparse = true) - @test Catalyst.issparse(A_k) numeqs = similar(eqs) for i in 1:length(eqs) diff --git a/test/runtests.jl b/test/runtests.jl index 55e1d5d7c8..2e2763e994 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,7 +15,7 @@ end ### Run Tests ### @time begin if GROUP == "All" || GROUP == "Core" - # Tests the `ReactionSystem` structure and its properties. + Tests the `ReactionSystem` structure and its properties. @time @safetestset "Reaction Structure" begin include("reactionsystem_core/reaction.jl") end @time @safetestset "ReactionSystem Structure" begin include("reactionsystem_core/reactionsystem.jl") end @time @safetestset "Higher Order Reactions" begin include("reactionsystem_core/higher_order_reactions.jl") end From a6d924a949d97ac3b756261f2b27cec05b2309ab Mon Sep 17 00:00:00 2001 From: vyudu Date: Sat, 14 Dec 2024 16:19:07 +0800 Subject: [PATCH 25/29] up --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 2e2763e994..55e1d5d7c8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,7 +15,7 @@ end ### Run Tests ### @time begin if GROUP == "All" || GROUP == "Core" - Tests the `ReactionSystem` structure and its properties. + # Tests the `ReactionSystem` structure and its properties. @time @safetestset "Reaction Structure" begin include("reactionsystem_core/reaction.jl") end @time @safetestset "ReactionSystem Structure" begin include("reactionsystem_core/reactionsystem.jl") end @time @safetestset "Higher Order Reactions" begin include("reactionsystem_core/higher_order_reactions.jl") end From b8e3255ebd2817ff97219a6f0f855c8aea5aef15 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Mon, 6 Jan 2025 11:07:05 -0500 Subject: [PATCH 26/29] Update test/runtests.jl --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 55e1d5d7c8..f251cd36d3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -49,7 +49,7 @@ end @time @safetestset "Automatic Jacobian Construction" begin include("simulation_and_solving/jacobian_construction.jl") end @time @safetestset "SDE System Simulations" begin include("simulation_and_solving/simulate_SDEs.jl") end @time @safetestset "Jump System Simulations" begin include("simulation_and_solving/simulate_jumps.jl") end - # @time @safetestset "Nonlinear and SteadyState System Solving" begin include("simulation_and_solving/solve_nonlinear.jl") end + @time @safetestset "Nonlinear and SteadyState System Solving" begin include("simulation_and_solving/solve_nonlinear.jl") end # Tests upstream SciML and DiffEq stuff. @time @safetestset "MTK Structure Indexing" begin include("upstream/mtk_structure_indexing.jl") end From aa12d240e52160b8ddce54a16d718c61f6977964 Mon Sep 17 00:00:00 2001 From: vyudu Date: Mon, 6 Jan 2025 12:22:21 -0500 Subject: [PATCH 27/29] revert to return , add docstrings --- src/network_analysis.jl | 56 ++++++++++++++++++++++++----------------- 1 file changed, 33 insertions(+), 23 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 841551125a..5d78c2b3ce 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -197,10 +197,12 @@ end Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the graph Laplacian, and ``Φ`` is the [`massactionvector`](@ref). ``A_k`` is an n-by-n matrix, where n is the number of complexes, where ``A_{ij} = k_{ij}`` if a reaction exists between the two complexes and 0 otherwise. Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. + + **Warning**: Unlike other Catalyst functions, the `laplacianmat` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs. """ -function laplacianmat(rn::ReactionSystem, pmap::Dict = Dict()) - D = incidencemat(rn) - K = fluxmat(rn, pmap) +function laplacianmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false) + D = incidencemat(rn; sparse) + K = fluxmat(rn, pmap; sparse) D * K end @@ -212,8 +214,10 @@ Base.one(::Type{Union{R, Symbolics.BasicSymbolic{Real}}}) where R <: Real = one( Return an r×c matrix ``K`` such that, if complex ``j`` is the substrate complex of reaction ``i``, then ``K_{ij} = k``, the rate constant for this reaction. Mostly a helper function for the network Laplacian, [`laplacianmat`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix and ``Φ(x)`` is the [`massactionvector`](@ref). Returns a symbolic matrix by default, but will return a numerical matrix if rate constants are specified as a `Tuple`, `Vector`, or `Dict` of symbol-value pairs via `pmap`. + + **Warning**: Unlike other Catalyst functions, the `fluxmat` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs. """ -function fluxmat(rn::ReactionSystem, pmap::Dict = Dict()) +function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) rates = if isempty(pmap) reactionrates(rn) else @@ -224,24 +228,28 @@ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict()) nc = length(rcmap) nr = length(rates) mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates) - Matrix{Any}(fluxmat(Matrix{mtype}, rcmap, rates)) -end - -# function fluxmat(::Type{SparseMatrixCSC{T, Int}}, rcmap, rates) where T -# Is = Int[] -# Js = Int[] -# Vs = T[] -# for (i, (complex, rxs)) in enumerate(rcmap) -# for (rx, dir) in rxs -# dir == -1 && begin -# push!(Is, rx) -# push!(Js, i) -# push!(Vs, rates[rx]) -# end -# end -# end -# Z = sparse(Is, Js, Vs, length(rates), length(rcmap)) -# end + if sparse + return fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates) + else + return fluxmat(Matrix{mtype}, rcmap, rates) + end +end + +function fluxmat(::Type{SparseMatrixCSC{T, Int}}, rcmap, rates) where T + Is = Int[] + Js = Int[] + Vs = T[] + for (i, (complex, rxs)) in enumerate(rcmap) + for (rx, dir) in rxs + dir == -1 && begin + push!(Is, rx) + push!(Js, i) + push!(Vs, rates[rx]) + end + end + end + Z = sparse(Is, Js, Vs, length(rates), length(rcmap)) +end function fluxmat(::Type{Matrix{T}}, rcmap, rates) where T nr = length(rates) @@ -279,6 +287,8 @@ end Return the vector whose entries correspond to the "mass action products" of each complex. For example, given the complex A + B, the corresponding entry of the vector would be ``A*B``, and for the complex 2X + Y, the corresponding entry would be ``X^2*Y``. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the [`laplacianmat`](@ref). This utility returns ``Φ(x)``. Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap. If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws). Will default to the default for the system. + + **Warning**: Unlike other Catalyst functions, the `massactionvector` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs. """ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) r = numreactions(rn) @@ -306,7 +316,7 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric push!(Φ, maprod) end - Vector{Any}(Φ) + Φ end function massactionvector(rn::ReactionSystem, scmap::Tuple; combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) From bfe849f09cc2b5ef68e19cc6565f5708f671702f Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 8 Jan 2025 13:09:18 -0500 Subject: [PATCH 28/29] modify docstring --- Project.toml | 4 ++-- src/network_analysis.jl | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 36eaf5d9e2..5e954e9a5a 100644 --- a/Project.toml +++ b/Project.toml @@ -58,13 +58,13 @@ JumpProcesses = "9.13.2" LaTeXStrings = "1.3.0" Latexify = "0.16.5" MacroTools = "0.5.5" -ModelingToolkit = "9.32" +ModelingToolkit = "< 9.60" NetworkLayout = "0.4.7" Parameters = "0.12" Reexport = "0.2, 1.0" Requires = "1.0" RuntimeGeneratedFunctions = "0.5.12" -SciMLBase = "2.57.2" +SciMLBase = "< 2.66.1" Setfield = "1" # StructuralIdentifiability = "0.5.8" SymbolicUtils = "2.1.2, 3.3.0" diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 5d78c2b3ce..18bbcdb6f9 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -198,7 +198,7 @@ end Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the graph Laplacian, and ``Φ`` is the [`massactionvector`](@ref). ``A_k`` is an n-by-n matrix, where n is the number of complexes, where ``A_{ij} = k_{ij}`` if a reaction exists between the two complexes and 0 otherwise. Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. - **Warning**: Unlike other Catalyst functions, the `laplacianmat` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs. + **Warning**: Unlike other Catalyst functions, the `laplacianmat` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs, and to ensure that multiplying the sparse form of the matrix will work. """ function laplacianmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false) D = incidencemat(rn; sparse) @@ -215,7 +215,7 @@ Base.one(::Type{Union{R, Symbolics.BasicSymbolic{Real}}}) where R <: Real = one( Return an r×c matrix ``K`` such that, if complex ``j`` is the substrate complex of reaction ``i``, then ``K_{ij} = k``, the rate constant for this reaction. Mostly a helper function for the network Laplacian, [`laplacianmat`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix and ``Φ(x)`` is the [`massactionvector`](@ref). Returns a symbolic matrix by default, but will return a numerical matrix if rate constants are specified as a `Tuple`, `Vector`, or `Dict` of symbol-value pairs via `pmap`. - **Warning**: Unlike other Catalyst functions, the `fluxmat` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs. + **Warning**: Unlike other Catalyst functions, the `fluxmat` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs, and to ensure that multiplying the sparse form of the matrix will work. """ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) rates = if isempty(pmap) @@ -288,7 +288,7 @@ end Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap. If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws). Will default to the default for the system. - **Warning**: Unlike other Catalyst functions, the `massactionvector` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs. + **Warning**: Unlike other Catalyst functions, the `massactionvector` function will return a `Vector{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs. """ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn)) r = numreactions(rn) From 6c259316c72796175acf1ec6ada6e81698e32f0e Mon Sep 17 00:00:00 2001 From: vyudu Date: Wed, 8 Jan 2025 13:26:20 -0500 Subject: [PATCH 29/29] update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 5e954e9a5a..6850c2f2a9 100644 --- a/Project.toml +++ b/Project.toml @@ -64,7 +64,7 @@ Parameters = "0.12" Reexport = "0.2, 1.0" Requires = "1.0" RuntimeGeneratedFunctions = "0.5.12" -SciMLBase = "< 2.66.1" +SciMLBase = "2.57.2" Setfield = "1" # StructuralIdentifiability = "0.5.8" SymbolicUtils = "2.1.2, 3.3.0"