Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

network_analysis cleanup: Updating ratematrix to be more similar to new matrix API functions #1172

Merged
merged 6 commits into from
Jan 17, 2025
Merged

network_analysis cleanup: Updating ratematrix to be more similar to new matrix API functions #1172

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
30cf870
simplify iscomplexbalanced, rename adjacencymat, adjust implementatio…
vyudu Jan 10, 2025
288041e
Merge branch 'master' into network-analysis-touchup
vyudu Jan 10, 2025
21f7d6e
Merge branch 'master' of https://github.com/SciML/Catalyst.jl into ne…
vyudu Jan 10, 2025
94a27c9
fix
vyudu Jan 10, 2025
12e7b2b
add validation
vyudu Jan 17, 2025
49a56fa
Update network_analysis.jl
vyudu Jan 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
125 changes: 59 additions & 66 deletions src/network_analysis.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -206,9 +206,6 @@ function laplacianmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false)
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)
vyudu marked this conversation as resolved.
Show resolved Hide resolved

@doc raw"""
fluxmat(rn::ReactionSystem, pmap = Dict(); sparse=false)

Expand All @@ -225,8 +222,6 @@ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false)
end

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)
Expand Down Expand Up @@ -936,7 +931,7 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict; abstol=0, re
complexes, D = reactioncomplexes(rs)
img = incidencematgraph(rs)
undir_img = SimpleGraph(incidencematgraph(rs))
K = ratematrix(rs, pmap)
K = adjacencymat(rs, pmap)

spanning_forest = Graphs.kruskal_mst(undir_img)
outofforest_edges = setdiff(collect(edges(undir_img)), spanning_forest)
Expand Down Expand Up @@ -993,52 +988,26 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap)
end

"""
iscomplexbalanced(rs::ReactionSystem, parametermap)
iscomplexbalanced(rs::ReactionSystem, parametermap; sparse = false)

Constructively compute whether a network will have complex-balanced equilibrium
solutions, following the method in van der Schaft et al., [2015](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3).
Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...].
"""
function iscomplexbalanced(rs::ReactionSystem, parametermap::Dict)
if length(parametermap) != numparams(rs)
error("Incorrect number of parameters specified.")
end

pmap = symmap_to_varmap(rs, parametermap)
pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap)

sm = speciesmap(rs)
cm = reactioncomplexmap(rs)
complexes, D = reactioncomplexes(rs)
Requires the specification of values for the parameters/rate constants. Accepts a dictionary, vector, or tuple of parameter-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...].
"""
function iscomplexbalanced(rs::ReactionSystem, parametermap::Dict; sparse = false)
rxns = reactions(rs)
nc = length(complexes)
nr = numreactions(rs)
nm = numspecies(rs)

if !all(r -> ismassaction(r, rs), rxns)
error("The supplied ReactionSystem has reactions that are not ismassaction. Testing for being complex balanced is currently only supported for pure mass action networks.")
end

rates = [substitute(rate, pmap) for rate in reactionrates(rs)]

# Construct kinetic matrix, K
K = zeros(nr, nc)
for c in 1:nc
complex = complexes[c]
for (r, dir) in cm[complex]
rxn = rxns[r]
if dir == -1
K[r, c] = rates[r]
end
end
end

L = -D * K
S = netstoichmat(rs)
L = laplacianmat(rs, parametermap; sparse)
vyudu marked this conversation as resolved.
Show resolved Hide resolved
D = incidencemat(rs; sparse)
S = netstoichmat(rs; sparse)

# Compute ρ using the matrix-tree theorem
g = incidencematgraph(rs)
R = ratematrix(rs, rates)
R = adjacencymat(rs, parametermap; sparse)
ρ = matrixtree(g, R)

# Determine if 1) ρ is positive and 2) D^T Ln ρ lies in the image of S^T
Expand Down Expand Up @@ -1069,50 +1038,74 @@ function iscomplexbalanced(rs::ReactionSystem, parametermap)
end

"""
ratematrix(rs::ReactionSystem, parametermap)
adjacencymat(rs::ReactionSystem, pmap = Dict(); sparse = false)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this require ismassaction, and if so is that being checked anywhere?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

More generally, it would be good if any of these methods that only make sense with fixed rate constants (in the Catalyst sense) return an informative error if they get an invalid system.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Have added them. I think fluxmat and adjacencymat should work with reactions with fixed rate constants that aren't necessarily mass action (which is the case for generalized mass action systems e.g., the CatalystNetworkAnalysis stuff does some stuff with that) and the massactionvector should work with just mass action reactions

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, if they should work more generally then the checks aren't needed. I just want to make sure we don't have people using them for systems they aren't intended / expected to work with.

LMK how you want to proceed (I'm fine keeping the checks or dropping them if you feel they aren't actually needed here).

Copy link
Collaborator Author

@vyudu vyudu Jan 17, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the check currently is good, we do still wanna check that stuff like k*A, A --> B isn't allowed in the adjacencymat/fluxmat since you could no longer use it in the ODE system decomposition if it is time-dependent. But in theory something like k, A --> B that doesn't have mass action kinetics (e.g. the rate law is k*A^2 or something) is fine and constructing the adjacencymat/fluxmat for such a reaction is useful in some cases.


Given a reaction system with n complexes, outputs an n-by-n matrix where R_{ij} is the rate
constant of the reaction between complex i and complex j. Accepts a dictionary, vector, or tuple
of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...].

Equivalent to the adjacency matrix of the weighted graph whose weights are the
reaction rates.
Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap.

The difference between this matrix and [`fluxmat`](@ref) is quite subtle. The non-zero entries to both matrices are the rate constants. However:
- In `fluxmat`, the rows represent complexes and columns represent reactions, and an entry (c, r) is non-zero if reaction r's source complex is c.
- In `adjacencymat`, the rows and columns both represent complexes, and an entry (c1, c2) is non-zero if there is a reaction c1 --> c2.
"""
function ratematrix(rs::ReactionSystem, rates::Vector{Float64})
complexes, D = reactioncomplexes(rs)
n = length(complexes)
rxns = reactions(rs)
ratematrix = zeros(n, n)
function adjacencymat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false)
rates = if isempty(pmap)
reactionrates(rn)
else
substitutevals(rn, pmap, parameters(rn), reactionrates(rn))
end
mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates)

for r in 1:length(rxns)
rxn = rxns[r]
s = findfirst(==(-1), @view D[:, r])
p = findfirst(==(1), @view D[:, r])
ratematrix[s, p] = rates[r]
if sparse
return adjacencymat(SparseMatrixCSC{mtype, Int}, incidencemat(rn), rates)
else
return adjacencymat(Matrix{mtype}, incidencemat(rn), rates)
end
ratematrix
end

function ratematrix(rs::ReactionSystem, parametermap::Dict)
if length(parametermap) != numparams(rs)
error("Incorrect number of parameters specified.")
function adjacencymat(::Type{SparseMatrixCSC{T, Int}}, D, rates) where T
Is = Int[]
Js = Int[]
Vs = T[]
nc = size(D, 1)

for r in 1:length(rates)
s = findfirst(==(-1), @view D[:, r])
p = findfirst(==(1), @view D[:, r])
push!(Is, s)
push!(Js, p)
push!(Vs, rates[r])
end
A = sparse(Is, Js, Vs, nc, nc)
end

pmap = symmap_to_varmap(rs, parametermap)
pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap)
function adjacencymat(::Type{Matrix{T}}, D, rates) where T
nc = size(D, 1)
A = zeros(T, nc, nc)

rates = [substitute(rate, pmap) for rate in reactionrates(rs)]
ratematrix(rs, rates)
for r in 1:length(rates)
s = findfirst(==(-1), @view D[:, r])
p = findfirst(==(1), @view D[:, r])
A[s, p] = rates[r]
end
A
end

function ratematrix(rs::ReactionSystem, parametermap::Vector{<:Pair})
pdict = Dict(parametermap)
ratematrix(rs, pdict)
function adjacencymat(rn::ReactionSystem, pmap::Vector{<:Pair}; sparse=false)
pdict = Dict(pmap)
adjacencymat(rn, pdict; sparse)
end

function ratematrix(rs::ReactionSystem, parametermap::Tuple)
pdict = Dict(parametermap)
ratematrix(rs, pdict)
function adjacencymat(rn::ReactionSystem, pmap::Tuple; sparse=false)
pdict = Dict(pmap)
adjacencymat(rn, pdict; sparse)
end

function ratematrix(rs::ReactionSystem, parametermap)
function adjacencymat(rn::ReactionSystem, pmap)
error("Parameter map must be a dictionary, tuple, or vector of symbol/value pairs.")
end

Expand Down
16 changes: 9 additions & 7 deletions test/network_analysis/crn_theory.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +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
let
Expand Down Expand Up @@ -69,18 +68,21 @@ let
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 Catalyst.adjacencymat(rn, rates_vec) == rate_mat
@test Catalyst.adjacencymat(rn, rates_tup) == rate_mat
@test Catalyst.adjacencymat(rn, rates_dict) == rate_mat
@test_throws Exception Catalyst.adjacencymat(rn, rate_vals)

# Tests that throws error in rate matrix.
incorrect_param_dict = Dict(:k1 => 1.0)

@test_throws ErrorException Catalyst.ratematrix(rn, 123)
@test_throws ErrorException Catalyst.ratematrix(rn, incorrect_param_dict)
@test_throws ErrorException Catalyst.adjacencymat(rn, 123)
@test_throws ErrorException Catalyst.adjacencymat(rn, incorrect_param_dict)

@test_throws Exception Catalyst.iscomplexbalanced(rn, rates_invalid)

# Test sparse matrix
@test Catalyst.adjacencymat(rn, rates_vec; sparse = true) == rate_mat
end

### CONCENTRATION ROBUSTNESS TESTS
Expand Down
Loading