`Linkage classes`, `deficiency`, Part 1

[yewalenikhil65]

These are a nice set of methods that could be useful in evaluating

• linkage classes in the network
• reversibility(weakly or strongly reversible)
• deficiency of the reaction network(once we know the number of linkage classes)

But would mean dependency on LightGraphs.jl (we already have such dependency in ModelingToolkit i believe). Is it advisable to pursue this with a view of PR for Catalyst? Or is there existing functionality related to this in Graphviz ?

[yewalenikhil65]

Rough idea for Linkage class

using LightGraphs
using Catalyst

function create_graph(incidence::Matrix{Int64})
           n = size(incidence)[1]
           graph = LightGraphs.DiGraph(n)    # find nodes where edge goes out and edge goes

           v_out = findall(x -> x == -1, incidence) #  edge goes out of vertex
           v_in = findall(x -> x == 1, incidence)  # edge goes in of vertex
           # match the Cartesian index between the two lists, add edge
           for i in v_in
               for o in v_out
                   if o[2] == i[2]
                       add_edge!(graph, o[1], i[1])
                   end
               end
           end
           return graph
       end


rn = @reaction_network begin
    k₁, 2A --> B
    k₂, A --> C
    k₃, C --> D
    k₄, B + D --> E
end k₁ k₂ k₃ k₄
B = reactioncomplexes(rn)[2]
g = create_graph(B)
julia> connected_components(g)
3-element Vector{Vector{Int64}}:
 [1, 2]
 [3, 4, 5]
 [6, 7]

where connected_components return the vectors of indices of complexes [2A, B, A,C,D,B+D,E]
This is apparent from complexgraph(rn) too

Meaning, we have number of linkage classes = length(connected_components(g)) = 3.
Thus , as per definition,deficiency of the network = number of complexes - number of linkage classes - rank of stoichiometry matrix

[yewalenikhil65]

To avoid dependency on LightGraphs.jl , we could referr to some code from https://github.com/JuliaGraphs/LightGraphs.jl/blob/master/src/SimpleGraphs/simpledigraph.jl

[isaacsas]

Adding a LightGraphs dependency is fine! (ModelingToolkit already has one, so it is going to be loaded there anyways...) I agree that getting these calculations in would be nice; I'd like to build to having some basic summarize type function that can report lots of info about the system just from the network properties. Knowing the deficiency is one piece in understanding certain steady-state properties.

Why do you need all the findall stuff in create_graph? Can't you just loop over the columns of the incidence matrix, and search through the rows saving the one -1 and the one 1 entry indices, and then adding an edge between them? That would be more memory efficient and likely faster. We'd now also need a version that works with a sparse incidence matrix.

Maybe a good name would be incidencegraph or incidencematgraph?

[yewalenikhil65]

Why do you need all the findall stuff in create_graph? Can't you just loop over the columns of the incidence matrix, and search through the rows saving the one -1 and the one 1 entry indices, and then adding an edge between them? That would be more memory efficient and likely faster. We'd now also need a version that works with a sparse incidence matrix.

findall was just my lazy way copying from some tutorial of graphs I was looking into. I will create a more efificent form in PR. both for dense and sparse versions

Maybe a good name would be incidencegraph or incidencematgraph?

incidencematgraph sounds cool. But I want to highlight that this graph and complexgraph Graphviz based function that we added are not different in essence.

[isaacsas]

Yeah, that is a potential problem with various names. I guess we can clean that up later; there aren't so many graph options at the moment that it should be a problem...

[yewalenikhil65]

Thus , as per definition,deficiency of the network = number of complexes - number of linkage classes - rank of stoichiometry matrix

rank of stoichiometry matrix in case of sparse version case, demands the matrix type to be SparseMatrixCSC{Float,Int} , as i think rank of sparse matrix in determined using qr algorithm. requiring Float entries

using plain convert(SparseMatrixCSC{Float,Int},A) is fine ?

[isaacsas]

Ultimately we probably need a sparse integer matrix method to directly calculate this. Maybe throw a warning and then convert to a dense matrix in the function? We can improve it later on.

I'd suggest calculating the rank using AbstractAlgebra.jl since it has methods that should work over the integers directly. i.e. something like

using AbstractAlgebra
S = your stoichiometry matrix
SI = matrix(zz, S)
r = rank(SI)

This might not be the most efficient, but it should be an integer-based method.

[yewalenikhil65]

I'd suggest calculating the rank using AbstractAlgebra.jl since it has methods that should work over the integers directly. i.e. something like

using AbstractAlgebra
S = your stoichiometry matrix
SI = matrix(zz, S)
r = rank(SI)

This might not be the most efficient, but it should be an integer-based method.

Thanks for AbstractAlgebra hint.. I am almost done with this issue n will raise PR tonight Just thinking about if I could add something on reversibility(weak or strong) too

[yewalenikhil65]

Closed with #394