Extend core functionality

Project.toml

@@ -9,16 +9,20 @@ ConstraintTrees = "5515826b-29c3-47a5-8849-8513ac836620"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 DistributedData = "f6a0035f-c5ac-4ad0-b410-ad102ced35df"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SBML = "e5567a89-2604-4b09-9718-f5f78e97c3bb"
+SBMLFBCModels = "3e8f9d1a-ffc1-486d-82d6-6c7276635980"
+SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 
 [compat]
 Clarabel = "0.3"
-ConstraintTrees = "0.2"
 DistributedData = "0.2"
 DocStringExtensions = "0.8, 0.9"
 JuMP = "1"

src/COBREXA.jl

@@ -4,12 +4,17 @@ $(README)
 module COBREXA
 
 using DocStringExtensions
-
+import AbstractFBCModels as A
 import ConstraintTrees as C
+import JuMP as J
+import SparseArrays: sparse
 
-include("types.jl")
+include("types/maybe.jl")
+include("types/isozyme.jl")
 include("solver.jl")
-include("builders/core.jl") #TODO more stuff
+include("builders/core.jl")
+include("builders/enzyme.jl")
+
 include("utils/downloads.jl")
 
 end # module COBREXA

src/builders/core.jl

@@ -1,29 +1,25 @@
 
-import AbstractFBCModels as F
-import SparseArrays: sparse
-
 """
 $(TYPEDSIGNATURES)
 
 A constraint tree that models the content of the given instance of
 `AbstractFBCModel`.
 """
-metabolic_model(model::F.AbstractFBCModel) =
-    let
-        rxns =
-            Symbol.(F.reactions(model)), mets =
-                Symbol.(F.metabolites(model)), lbs, ubs =
-                    F.bounds(model), stoi =
-                        F.stoichiometry(model), bal =
-                            F.balance(model), obj = F.objective(model)
+function metabolic_model(model::A.AbstractFBCModel)
+    rxns = Symbol.(A.reactions(model))
+    mets = Symbol.(A.metabolites(model))
+    lbs, ubs = A.bounds(model)
+    stoi = A.stoichiometry(model)
+    bal = A.balance(model)
+    obj = A.objective(model)
 
-        :fluxes^C.variables(keys = rxns, bounds = zip(lbs, ubs)) *
-        :balance^C.ConstraintTree(
-            m => Constraint(value = Value(sparse(row)), bound = b) for
-            (m, row, b) in zip(mets, eachrow(stoi), bals)
-        ) *
-        :objective^C.Constraint(C.Value(sparse(obj)))
-    end
+    :fluxes^C.variables(keys = rxns, bounds = zip(lbs, ubs)) *
+    :flux_stoichiometry^C.ConstraintTree(
+        met => C.Constraint(value = C.LinearValue(sparse(row)), bound = b) for
+        (met, row, b) in zip(mets, eachrow(stoi), bal)
+    ) *
+    :objective^C.Constraint(C.LinearValue(sparse(obj)))
+end
 
 """
 $(TYPEDSIGNATURES)
@@ -33,6 +29,18 @@ Shortcut for allocation non-negative ("unsigned") variables. The argument
 """
 unsigned_variables(; keys) = C.variables(; keys, bounds = Ref((0.0, Inf)))
 
+function fluxes_in_direction(fluxes::C.ConstraintTree, direction = :forward)
+    keys = Symbol[]
+    for (id, flux) in fluxes
+        if direction == :forward
+            last(flux.bound) > 0 && push!(keys, id)
+        else
+            first(flux.bound) < 0 && push!(keys, id)
+        end
+    end
+    C.variables(; keys, bounds = Ref((0.0, Inf)))
+end
+
 """
 $(TYPEDSIGNATURES)
 
@@ -44,7 +52,7 @@ Keys in the result are the same as the keys of `signed` constraints.
 Typically, this can be used to create "unidirectional" fluxes
 together with [`unsigned_variables`](@ref):
 ```
-uvars = unsigned_variables(keys(myModel.fluxes))
+uvars = unsigned_variables(keys = collect(keys(myModel.fluxes)))
 
 myModel = myModel +
     :fluxes_forward^uvars +
@@ -64,9 +72,10 @@ sign_split_constraints(;
     negative::C.ConstraintTree,
     signed::C.ConstraintTree,
 ) = C.ConstraintTree(
-    C.Constraint(
-        value = s + (haskey(negative, k) ? negative[k].value : zero(C.Value)) -
-                (haskey(positive, k) ? positive[k].value : zero(C.Value)),
+    k => C.Constraint(
+        value = s.value +
+                (haskey(negative, k) ? negative[k].value : zero(typeof(s.value))) -
+                (haskey(positive, k) ? positive[k].value : zero(typeof(s.value))),
         bound = 0.0,
-    ) for (k, s) in C.elems(signed)
+    ) for (k, s) in signed
 )

src/builders/enzyme.jl

@@ -0,0 +1,132 @@
+"""
+$(TYPEDSIGNATURES)
+
+Allocate enzyme variables (gene products in the model) to a constraint tree.
+"""
+enzyme_variables(model::A.AbstractFBCModel) =
+    C.variables(; keys = Symbol.(A.genes(model)), bounds = Ref((0.0, Inf)))
+
+"""
+$(TYPEDSIGNATURES)
+
+Create isozyme variables for reactions. A single reaction may be catalyzed by
+multiple enzymes (isozymes), and the total flux through a reaction is the sum
+through of all these isozymes. These variables are linked to fluxes through
+[`link_isozymes`](@ref).
+"""
+function isozyme_variables(
+    reaction_id::Symbol,
+    reaction_isozymes::Dict{Symbol,Dict{Symbol, Isozyme}},
+)
+    C.variables(;
+        keys = collect(keys(reaction_isozymes[reaction_id])),
+        bounds = Ref((0.0, Inf)),
+    )
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Link isozymes to fluxes. All forward (backward) fluxes are the sum of all the
+isozymes catalysing these fluxes.
+"""
+function link_isozymes(
+    fluxes_directional::C.ConstraintTree,
+    fluxes_isozymes::C.ConstraintTree,
+)
+    C.ConstraintTree(
+        k => C.Constraint(
+            value = s.value - sum(x.value for (_, x) in fluxes_isozymes[k]),
+            bound = 0.0,
+        ) for (k, s) in fluxes_directional if haskey(fluxes_isozymes, k)
+    )
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Create the enzyme "mass balance" matrix. In essence, the stoichiometric
+coefficient is subunit_stoichiometry / kcat for each directional, isozyme flux,
+and it must be balanced by the enzyme variable supply.
+"""
+function enzyme_stoichiometry(
+    enzymes::C.ConstraintTree,
+    fluxes_isozymes_forward::C.ConstraintTree,
+    fluxes_isozymes_backward::C.ConstraintTree,
+    reaction_isozymes::Dict{Symbol,Dict{Symbol, Isozyme}},
+)   
+    # map enzyme ids to reactions that use them (not all isozymes have to though)
+    enzyme_rid_lookup = Dict{Symbol,Vector{Symbol}}()
+    for (rid, isozymes) in reaction_isozymes
+        for isozyme in values(isozymes)
+            for gid in keys(isozyme.gene_product_stoichiometry)
+                rids = get!(enzyme_rid_lookup, Symbol(gid), Symbol[])
+                rid in rids || push!(rids, rid)
+            end
+        end
+    end
+
+    C.ConstraintTree(
+        gid => C.Constraint(
+            value = enz.value + # supply
+                    sum(
+                        enzyme_balance(
+                            gid,
+                            rid,
+                            fluxes_isozymes_forward,
+                            reaction_isozymes,
+                            :kcat_forward,
+                        ) for rid in enzyme_rid_lookup[gid] if
+                        haskey(fluxes_isozymes_forward, rid);
+                        init = zero(typeof(enz.value)),
+                    ) + # flux through positive isozymes
+                    sum(
+                        enzyme_balance(
+                            gid,
+                            rid,
+                            fluxes_isozymes_backward,
+                            reaction_isozymes,
+                            :kcat_backward,
+                        ) for rid in enzyme_rid_lookup[gid] if
+                        haskey(fluxes_isozymes_backward, rid);
+                        init = zero(typeof(enz.value)),
+                    ), # flux through negative isozymes
+            bound = 0.0,
+        ) for (gid, enz) in enzymes if gid in keys(enzyme_rid_lookup)
+    )
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Helper function to balance the forward or backward isozyme fluxes for a specific
+gene product.
+"""
+function enzyme_balance(
+    gid::Symbol,
+    rid::Symbol,
+    fluxes_isozymes::C.ConstraintTree, # direction
+    reaction_isozymes::Dict{Symbol,Dict{Symbol, Isozyme}},
+    direction = :kcat_forward,
+)   
+    isozyme_dict = reaction_isozymes[rid]
+
+    sum( # this is where the stoichiometry comes in
+        -isozyme_value.value * isozyme_dict[isozyme_id].gene_product_stoichiometry[string(gid)] /
+        getfield(isozyme_dict[isozyme_id], direction) for (isozyme_id, isozyme_value) in
+        fluxes_isozymes[rid] if gid in Symbol.(keys(isozyme_dict[isozyme_id].gene_product_stoichiometry));
+        init = zero(C.LinearValue),
+    )
+end
+
+function enzyme_capacity(
+    enzymes::C.ConstraintTree, 
+    enzyme_molar_mass::Dict{Symbol, Float64}, 
+    enzyme_ids::Vector{Symbol}, 
+    capacity::Float64,
+)   
+    C.Constraint(
+        value = sum(enzymes[gid].value * enzyme_molar_mass[gid] for gid in enzyme_ids),
+        bound = (0.0, capacity),
+    )
+end