up

Project.toml

@@ -4,6 +4,7 @@ authors = ["Vincent Du"]
 version = "0.1.0"
 
 [deps]
+CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60"
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
@@ -24,9 +25,11 @@ ReactionNetworkImporters = "b4db0fb7-de2a-5028-82bf-5021f5cfa881"
 SBMLToolkit = "86080e66-c8ac-44c2-a1a0-9adaadfe4a4e"
 Satisfiability = "160ab843-0bc6-4ba4-9585-b7478b70f443"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 
 [compat]
+CDDLib = "0.9.4"
 Catalyst = "14.4.1"
 DynamicPolynomials = "0.6.0"
 Graphs = "1.12.0"
@@ -40,20 +43,22 @@ ModelingToolkit = "9.46.1"
 MultivariatePolynomials = "0.5.7"
 Nemo = "0.47"
 Oscar = "1.1.1"
+Polyhedra = "0.7.8"
 PolynomialRoots = "1.0.0"
 ReactionNetworkImporters = "0.15.1"
 SBMLToolkit = "0.1.29"
 SafeTestsets = "0.1.0"
 Satisfiability = "0.2"
 SparseArrays = "1"
+StaticArrays = "1.9.8"
 SymbolicUtils = "3.7.2"
 julia = "1.10"
 
 [extras]
 SBMLImporter = "210efffb-c3c8-456d-a807-6f55560b12fe"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [targets]

src/CatalystNetworkAnalysis.jl

@@ -11,6 +11,9 @@ using SparseArrays
 using IterTools
 using StaticArrays 
 
+using Polyhedra
+import CDDLib
+
 import ModelingToolkit as MT
 
 const M::Float64 = 1E4

src/lp_utils.jl

@@ -113,3 +113,21 @@ function is_extreme_idxset(S::M, idxs::Vector{Int}) where {M <: AbstractMatrix}
     cone_mat_eq = @view cone_mat[rowidxs, :]
     return rank(cone_mat_eq) == n - 1
 end
+
+"""
+    elementaryfluxmodes(rn::ReactionSystem)
+
+    Given a reaction network, return the set of elementary flux modes of the reaction network. 
+"""
+function elementary_flux_modes(rn::ReactionSystem)
+    S = netstoichmat(rn)
+    m, n = size(S)
+    hyperplanes = [Polyhedra.HyperPlane(S[i, :], 0) for i in 1:m]
+    halfspaces = [Polyhedra.HalfSpace(-I(n)[i, :], 0) for i in 1:n]
+    polycone = Polyhedra.polyhedron(hrep(hyperplanes, halfspaces))
+
+    Polyhedra.vrep(polycone)
+    Polyhedra.removevredundancy!(polycone)
+
+    EFMs = reduce(hcat, map(x->x.a, polycone.vrep.rays.rays))
+end

src/persistence.jl

@@ -4,7 +4,6 @@
 #   permanence in networks
 #   improve efficiency of siphon detection
 
-
 """
     ispersistent(rs::ReactionSystem)
 
@@ -94,9 +93,9 @@ function minimalsiphons_smt(rs::ReactionSystem)
     return removesupersets(siphons)
 end
 
-function removesupersets(indexsets::AbstractArray{<:Array})
+function removesupersets(indexsets)
     indexsets = sort(indexsets, by=x->length(x))
-    minimalsets = Array[]
+    minimalsets = Array{Int64}[]
 
     for s in indexsets
         if !any(ms->issubset(ms, s), minimalsets)
@@ -148,5 +147,5 @@ function iscritical(siphon::Vector, S::Matrix)
 
     # If there is a non-negative vector in the nullspace of S_red', then there is a positive conservation law with a support that is the subset of the siphon, and the siphon is not critical 
     conslaws_r = conservationlaws(S_r)
-    !haspositivesolution(copy(conslaws_r'), nonneg = true)
+    !has_positive_solution(copy(conslaws_r'), nonneg = true)
 end

src/translated.jl

@@ -7,8 +7,8 @@
 struct Translation{T<:Int} 
     rn::ReactionSystem
     Y_T::Matrix{T} # The new complex-composition matrix.
-    δ_K::T
-    δ_E::T
+    δ_K::T # Kinetic deficiency 
+    δ_E::T # Effective deficiency
 end
 
 """
@@ -35,10 +35,12 @@ function WRDZ_translation(rn::ReactionSystem)
 
         # Update reactant complexes
         @label step_3
-        rxs = findall(==(1), rr_adj[i, :])
+        rxs = findall(==(1), @view rr_adj[i, :])
         for j in rxs
             p, s = (dst(e[i]), src(e[j]))
-            (@view Y_T[:, j]) .= Y_K[:, p] - Y_K[:, s] + Y_T[:, i]
+            for k in 1:size(Y_T, 1)
+                Y_T[k, j] = Y_K[k, p] - Y_K[k, s] + Y_T[k, i]
+            end
             P3 = setdiff(push!(P3, j), P2)
         end
 
@@ -113,7 +115,7 @@ function rrgraph(rn::ReactionSystem)
     end
 
     optimize!(model)
-    is_solved_and_feasible(model) ? return edge : return
+    is_solved_and_feasible(model) && return edge 
 end
 
 # Generate partitions of the reaction-to-reaction graph

src/utils.jl

@@ -18,7 +18,7 @@ end
 """
     removespec(rn, spec)
 
-    Return the net stoichiometric matrix and incidence matrix obtained when removing a species at index spec from the reaction network. 
+    Return the net stoichiometric matrix and incidence matrix obtained when removing a species at index spec from the reaction network.
 """
 function removespec(rn::ReactionSystem, spec::Int64) 
     s, r = size(netstoichmat(rn))
@@ -39,9 +39,8 @@ end
 """
     transitiveclosure(arr, relation)
 
-    Given an iterable array and a (symmetric, reflexive) relation (a function that takes two objects, and returns true if they belong to the relation), return the partitions of the set under the transitive closure of the given relation. 
+    Given an iterable array and a (symmetric, reflexive) relation (a function that takes two objects, and returns true if they belong to the relation), return the partitions of the set under the transitive closure of the given relation.
 """
-
 function transitiveclosure(arr, relation) 
     adjmat = SparseArrays.spzeros(Bool, length(arr), length(arr))
     for idx in CartesianIndices(adjmat)

test/lp_utils.jl

@@ -0,0 +1,29 @@
+# Testing of the linear programming utilities in CatalystNetworkAnalysis
+
+import CatalystNetworkAnalysis as C
+# Test that matrices with positive element in their image space are identified.
+let
+    a = rand(1:10, 10, 10)
+    @test has_positive_solution(a)
+
+    b = -rand(1:10, 10, 10)
+    @test has_positive_solution(b)
+
+    a = [1 -1; -1 1]
+    @test !has_positive_solution(b)
+
+    b = [-1 1 1 1;
+         1 -1 1 1;
+         1 1 -1 1
+         1 1 1 -1]
+    @test !has_positive_solution(b)
+    @test has_positive_solution(b, nonneg=true)
+end
+
+# Testing that extreme rays are properly identified.
+let
+end
+
+# Testing elementary flux modes are found.
+let
+end