Merge pull request #893 from vyudu/cycles

Adding cycles/flux mode calculation

src/network_analysis.jl

@@ -639,22 +639,9 @@ end
 Given the net stoichiometry matrix of a reaction system, computes a matrix of
 conservation laws, each represented as a row in the output.
 """
-function conservationlaws(nsm::T; col_order = nothing) where {T <: AbstractMatrix}
-
-    # compute the left nullspace over the integers
-    N = MT.nullspace(nsm'; col_order)
-
-    # if all coefficients for a conservation law are negative, make positive
-    for Nrow in eachcol(N)
-        all(r -> r <= 0, Nrow) && (Nrow .*= -1)
-    end
-
-    # check we haven't overflowed
-    iszero(N' * nsm) || error("Calculation of the conservation law matrix was inaccurate, "
-          * "likely due to numerical overflow. Please use a larger integer "
-          * "type like Int128 or BigInt for the net stoichiometry matrix.")
-
-    T(N')
+function conservationlaws(nsm::Matrix; col_order = nothing)
+    conslaws = positive_nullspace(nsm'; col_order = col_order)
+    Matrix(conslaws)
 end
 
 # Used in the subsequent function.
@@ -905,3 +892,49 @@ function treeweight(t::SimpleDiGraph, g::SimpleDiGraph, distmx::Matrix)
     end
     prod
 end
+
+"""
+    cycles(rs::ReactionSystem)
+
+    Returns the matrix of a basis of cycles (or flux vectors), or a basis for reaction fluxes for which the system is at steady state. 
+    These correspond to right eigenvectors of the stoichiometric matrix. Equivalent to [`fluxmodebasis`](@ref). 
+"""
+
+function cycles(rs::ReactionSystem)
+    nps = get_networkproperties(rs)
+    nsm = netstoichmat(rs)
+    !isempty(nps.cyclemat) && return nps.cyclemat
+    nps.cyclemat = cycles(nsm; col_order = nps.col_order)
+    nps.cyclemat
+end
+
+function cycles(nsm::Matrix; col_order = nothing)
+    positive_nullspace(nsm; col_order)
+end
+
+function positive_nullspace(M::T; col_order = nothing) where {T <: AbstractMatrix}
+    # compute the left nullspace over the integers
+    N = MT.nullspace(M; col_order)
+
+    # if all coefficients for a cycle are negative, make positive
+    for Ncol in eachcol(N)
+        all(r -> r <= 0, Ncol) && (Ncol .*= -1)
+    end
+
+    # check we haven't overflowed
+    iszero(M * N) || error("Calculation of the cycle matrix was inaccurate, "
+          * "likely due to numerical overflow. Please use a larger integer "
+          * "type like Int128 or BigInt for the net stoichiometry matrix.")
+
+    T(N)
+end
+
+"""
+    fluxvectors(rs::ReactionSystem)
+
+    See documentation for [`cycles`](@ref). 
+"""
+
+function fluxvectors(rs::ReactionSystem)
+    cycles(rs)
+end

test/network_analysis/network_properties.jl

@@ -72,6 +72,16 @@ let
     #     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
+        if S[:,i] == -S[:,i+1]
+           cycle = [(j == i) || (j == i+1) ? 1 : 0 for j in 1:size(S,2)]
+           @test rank(cyclemat) == rank(hcat(cyclemat, cycle))
+        end
+    end
 end
 
 # Tests network analysis functions on a second network (by comparing to manually computed outputs).
@@ -349,7 +359,6 @@ end
 
 ### STRONG LINKAGE CLASS TESTS
 
-
 # a) Checks that strong/terminal linkage classes are correctly found. Should identify the (A, B+C) linkage class as non-terminal, since B + C produces D
 let
     rn = @reaction_network begin
@@ -431,6 +440,64 @@ let
     @test issubset([[3,4], [5,6,7]], tslcs) 
 end
 
+# Cycle Test: Open Reaction Network
+let
+    rn = @reaction_network begin
+        k1, 0 --> X1
+        k2, X1 --> 0
+        k3, X1 --> X2
+        (k4, k5), X2 <--> X3
+        (k6, k7), X3 <--> 0
+    end
+
+    # 0 --> X1 --> X2 --> X3 --> 0
+    cycle = [1, 0, 1, 1, 0, 1, 0]
+    cyclemat = Catalyst.cycles(rn)
+    @test rank(cyclemat) == rank(hcat(cyclemat, cycle))
+end
+
+# From stoichiometric matrix. Reference: Trinh, 2008, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2909134/
+let
+   S = [1 -1 0 0 -1 0 0 0 0;
+        0 0 0 0 1 -1 -1 -1 0;
+        0 1 -1 0 0 1 0 0 0;
+        0 0 1 0 0 0 0 0 -1;
+        0 0 1 -1 0 0 2 0 0]
+
+   EFMs = [1 0 1 1 0 1 1 1;
+           1 0 0 1 0 0 1 0;
+           0 1 0 1 0 0 0 1;
+           0 1 0 1 2 2 2 1;
+           0 0 1 0 0 1 0 1;
+           -1 1 0 0 0 0 -1 1;
+           0 0 0 0 1 1 1 0;
+           1 -1 1 0 -1 0 0 0;
+           0 1 0 1 0 0 0 1]
+
+   cyclemat = Catalyst.cycles(S)
+   for i in 1:size(EFMs, 2)
+       EFM = EFMs[:, i]
+       @test rank(cyclemat) == rank(hcat(cyclemat, EFM))
+   end
+end
+
+# No cycles should exist in the following network (the graph is treelike and irreversible)
+
+let
+    rn = @reaction_network begin
+        k1, A + B --> C + D
+        k2, C + D --> E + F
+        k3, C + D --> 2G + H
+        k4, 2G + H --> 3I
+        k5, E + F --> J 
+        k6, 3I --> K
+    end
+
+    S = netstoichmat(rn)
+    cyclemat = Catalyst.cycles(S)
+    @test isempty(cyclemat)
+end
+
 ### Other Network Properties Tests ###
 
 # Tests outgoing complexes matrices (1).