implement extended SimplifiedEnzymeConstrainedModel

src/analysis.jl

@@ -24,7 +24,7 @@ using ..Internal: constants
 using ..Log.Internal: @models_log
 using ..Solver
 using ..Types
-using ..Types: _EnzymeConstrainedReactionColumn, _SimplifiedEnzymeConstrainedColumn
+using ..Types: _EnzymeConstrainedReactionColumn, SimplifiedEnzymeConstrainedColumn
 
 using Distributed
 using DistributedData

src/reconstruction/simplified_enzyme_constrained.jl

@@ -1,25 +1,59 @@
 """
 $(TYPEDSIGNATURES)
 
-Wrap a `model` with a structure given by sMOMENT algorithm; returns a
+Wrap a `model` with a structure given by sMOMENT algorithm. Returns a
 [`SimplifiedEnzymeConstrainedModel`](@ref) (see the documentation for details).
 The sMOMENT algorithm only uses one [`Isozyme`](@ref) per reaction. If multiple
 isozymes are present in `model`, the "fastest" isozyme will be used. This is
 determined based on maximum kcat (forward or backward) divided by mass of the
-isozyme. The `total_gene_product_mass_bound` is the maximum "enzyme capacity" in
-the model.
+isozyme. Multiple enzyme capacity constraint can be placed on the model using
+the keyword arguments.
+
+Parameters `reaction_mass_groups` and `reaction_mass_group_bounds` specify
+groups of reactions (by their IDs), and their respective total mass limit.
+Reactions that are not listed in any reaction mass group are ignored (likewise
+if they don't have isozymes).
+
+For simplicity, specifying the `total_reaction_mass_bound` argument overrides
+the above arguments by internally specifying a single group on all reactions
+that acts like a maximum "enzyme capacity" in the model.
 
 # Example
 ```
 ecmodel = make_simplified_enzyme_constrained_model(
     model;
-    total_gene_product_mass_bound = 0.5
+    reaction_mass_groups = Dict(
+        "membrane" => ["r1", "r2"],
+        "total" => ["r1", "r2", "r3"],
+    ),
+    reaction_mass_group_bounds = Dict(
+        "membrane" => 0.2,
+        "total" => 0.5,
+    ),
+)
+
+ecmodel2 = make_simplified_enzyme_constrained_model(
+    model;
+    total_reaction_mass_bound = 0.5
 )
+```
 """
 function make_simplified_enzyme_constrained_model(
     model::AbstractMetabolicModel;
-    total_gene_product_mass_bound::Float64 = 0.5,
+    total_reaction_mass_bound::Maybe{Float64} = nothing,
+    reaction_mass_groups::Maybe{Dict{String,Vector{String}}} = nothing,
+    reaction_mass_group_bounds::Maybe{Dict{String,Float64}} = nothing,
 )
+    # fix kwarg inputs
+    if !isnothing(total_reaction_mass_bound)
+        reaction_mass_groups = Dict("uncategorized" => variables(model)) # TODO should be reactions
+        reaction_mass_group_bounds = Dict("uncategorized" => total_reaction_mass_bound)
+    end
+    isnothing(reaction_mass_groups) &&
+        throw(ArgumentError("missing reaction mass group specification"))
+    isnothing(reaction_mass_group_bounds) &&
+        throw(ArgumentError("missing reaction mass group bounds"))
+
     # helper function to rank the isozymes by relative speed
     speed_enzyme(model, isozyme) =
         max(isozyme.kcat_forward, isozyme.kcat_backward) / sum(
@@ -34,20 +68,23 @@ function make_simplified_enzyme_constrained_model(
         argmax(isozyme -> speed_enzyme(model, isozyme), isozymes)
     end
 
-    columns = Vector{Types._SimplifiedEnzymeConstrainedColumn}()
+    columns = Vector{Types.SimplifiedEnzymeConstrainedColumn}()
 
-    (lbs, ubs) = bounds(model)
-    rids = variables(model)
+    bound_ids = keys(reaction_mass_group_bounds)
+    total_reaction_mass_bounds = collect(values(reaction_mass_group_bounds))
 
-    for i = 1:n_variables(model)
+    (lbs, ubs) = bounds(model) # TODO need a reaction_bounds accessor for full generality
+    rids = variables(model) # TODO needs to be reactions
+
+    for i = 1:n_variables(model) # TODO this should be reactions 
 
         isozyme = ris_(model, rids[i])
 
         if isnothing(isozyme)
             # non-enzymatic reaction (or a totally ignored one)
             push!(
                 columns,
-                Types._SimplifiedEnzymeConstrainedColumn(i, 0, lbs[i], ubs[i], 0),
+                Types.SimplifiedEnzymeConstrainedColumn(i, 0, lbs[i], ubs[i], 0, Int64[]),
             )
             continue
         end
@@ -57,16 +94,22 @@ function make_simplified_enzyme_constrained_model(
             (gid, ps) in isozyme.gene_product_stoichiometry
         )
 
+        bidxs = [
+            idx for
+            (idx, bid) in enumerate(bound_ids) if rids[i] in reaction_mass_groups[bid]
+        ]
+
         if min(lbs[i], ubs[i]) < 0 && isozyme.kcat_backward > constants.tolerance
             # reaction can run in reverse
             push!(
                 columns,
-                Types._SimplifiedEnzymeConstrainedColumn(
+                Types.SimplifiedEnzymeConstrainedColumn(
                     i,
                     -1,
                     max(-ubs[i], 0),
                     -lbs[i],
                     mw / isozyme.kcat_backward,
+                    bidxs,
                 ),
             )
         end
@@ -75,16 +118,17 @@ function make_simplified_enzyme_constrained_model(
             # reaction can run forward
             push!(
                 columns,
-                Types._SimplifiedEnzymeConstrainedColumn(
+                Types.SimplifiedEnzymeConstrainedColumn(
                     i,
                     1,
                     max(lbs[i], 0),
                     ubs[i],
                     mw / isozyme.kcat_forward,
+                    bidxs,
                 ),
             )
         end
     end
 
-    return SimplifiedEnzymeConstrainedModel(columns, total_gene_product_mass_bound, model)
+    return SimplifiedEnzymeConstrainedModel(columns, total_reaction_mass_bounds, model)
 end

src/types/wrappers/SimplifiedEnzymeConstrainedModel.jl

@@ -7,12 +7,13 @@ A helper type that describes the contents of [`SimplifiedEnzymeConstrainedModel`
 # Fields
 $(TYPEDFIELDS)
 """
-struct _SimplifiedEnzymeConstrainedColumn
+struct SimplifiedEnzymeConstrainedColumn
     reaction_idx::Int # number of the corresponding reaction in the inner model
     direction::Int # 0 if "as is" and unique, -1 if reverse-only part, 1 if forward-only part
     lb::Float64 # must be 0 if the reaction is unidirectional (if direction!=0)
     ub::Float64
-    capacity_required::Float64 # must be 0 for bidirectional reactions (if direction==0)
+    capacity_contribution::Float64 # must be 0 for bidirectional reactions (if direction==0)
+    capacity_bound_idxs::Vector{Int64} # index of associated bound(s)
 end
 
 """
@@ -29,24 +30,24 @@ The model is constructed as follows:
 - stoichiometry of the original model is retained as much as possible, but
   enzymatic reations are split into forward and reverse parts (marked by a
   suffix like `...#forward` and `...#reverse`),
-- coupling is added to simulate a virtual metabolite "enzyme capacity", which is
-  consumed by all enzymatic reactions at a rate given by enzyme mass divided by
-  the corresponding kcat,
-- the total consumption of the enzyme capacity is constrained to a fixed
-  maximum.
+- coupling is added to simulate lumped virtual metabolites that act like enzyme
+  capacities. These are consumed by enzymatic reactions at a rate given by
+  enzyme mass divided by the corresponding kcat,
+- the total consumption of the enzyme capacity bounds is constrained to be less
+  than some fixed values.
 
 The `SimplifiedEnzymeConstrainedModel` structure contains a worked-out
 representation of the optimization problem atop a wrapped
-[`AbstractMetabolicModel`](@ref), in particular the separation of certain
-reactions into unidirectional forward and reverse parts (which changes the
-stoichiometric matrix), an "enzyme capacity" required for each reaction, and the
-value of the maximum capacity constraint. Original coupling in the inner model
-is retained.
+[`AbstractMetabolicModel`](@ref). The internal representation of the model
+splits reactions into unidirectional forward and reverse parts (which changes
+the stoichiometric matrix). 
 
 In the structure, the field `columns` describes the correspondence of
-stoichiometry columns to the stoichiometry and data of the internal wrapped
-model, and `total_gene_product_mass_bound` is the total bound on the enzyme
-capacity consumption as specified in sMOMENT algorithm.
+stoichiometry columns to the stoichiometry, and data of the internal wrapped
+model. Multiple capacity bounds may be added through
+`total_reaction_mass_bounds`. These bounds are connected to the model through
+`columns`. Since this algorithm is reaction centered, no enzymes directly appear
+in the formulation.
 
 This implementation allows easy access to fluxes from the split reactions
 (available in `variables(model)`), while the original "simple" reactions from
@@ -64,8 +65,8 @@ enzymes, or those that should be ignored need to return `nothing` when
 $(TYPEDFIELDS)
 """
 struct SimplifiedEnzymeConstrainedModel <: AbstractModelWrapper
-    columns::Vector{_SimplifiedEnzymeConstrainedColumn}
-    total_gene_product_mass_bound::Float64
+    columns::Vector{SimplifiedEnzymeConstrainedColumn}
+    total_reaction_mass_bounds::Vector{Float64}
 
     inner::AbstractMetabolicModel
 end
@@ -118,15 +119,37 @@ Accessors.reaction_variables(model::SimplifiedEnzymeConstrainedModel) =
         reaction_variables_matrix(model),
     ) # TODO currently inefficient
 
-Accessors.coupling(model::SimplifiedEnzymeConstrainedModel) = vcat(
-    coupling(model.inner) * simplified_enzyme_constrained_column_reactions(model),
-    [col.capacity_required for col in model.columns]',
-)
+function Accessors.coupling(model::SimplifiedEnzymeConstrainedModel)
+    inner_coupling =
+        coupling(model.inner) * simplified_enzyme_constrained_column_reactions(model)
+
+    I = Int64[]
+    J = Int64[]
+    V = Float64[]
+    for (col_idx, col) in enumerate(model.columns)
+        for row_idx in col.capacity_bound_idxs
+            push!(J, col_idx)
+            push!(I, row_idx)
+            push!(V, col.capacity_contribution)
+        end
+    end
+
+    capacity_coupling =
+        sparse(I, J, V, length(model.total_reaction_mass_bounds), length(model.columns))
+
+    return [
+        inner_coupling
+        capacity_coupling
+    ]
+end
 
 Accessors.n_coupling_constraints(model::SimplifiedEnzymeConstrainedModel) =
-    n_coupling_constraints(model.inner) + 1
+    n_coupling_constraints(model.inner) + length(model.total_reaction_mass_bounds)
 
 Accessors.coupling_bounds(model::SimplifiedEnzymeConstrainedModel) =
     let (iclb, icub) = coupling_bounds(model.inner)
-        (vcat(iclb, [0.0]), vcat(icub, [model.total_gene_product_mass_bound]))
+        (
+            vcat(iclb, zeros(length(model.total_reaction_mass_bounds))),
+            vcat(icub, model.total_reaction_mass_bounds),
+        )
     end

test/reconstruction/simplified_enzyme_constrained.jl

@@ -34,7 +34,7 @@
             lower_bounds = [-1000.0, -1.0],
             upper_bounds = [nothing, 12.0],
         ) |>
-        with_simplified_enzyme_constraints(total_gene_product_mass_bound = 100.0)
+        with_simplified_enzyme_constraints(total_reaction_mass_bound = 100.0)
 
     rxn_fluxes =
         flux_balance_analysis(
@@ -51,7 +51,7 @@
 end
 
 @testset "Small SMOMENT" begin
-    
+
     m = ObjectModel()
     m1 = Metabolite("m1")
     m2 = Metabolite("m2")
@@ -77,13 +77,15 @@ end
         Gene(id = "g4", product_upper_bound = 10.0, product_molar_mass = 4.0)
     ]
 
-    m.reactions["r3"].gene_associations = [Isozyme(["g1"]; kcat_forward = 1.0, kcat_backward = 1.0)]
-    m.reactions["r4"].gene_associations = [Isozyme(["g1"]; kcat_forward = 2.0, kcat_backward = 2.0)]
+    m.reactions["r3"].gene_associations =
+        [Isozyme(["g1"]; kcat_forward = 1.0, kcat_backward = 1.0)]
+    m.reactions["r4"].gene_associations =
+        [Isozyme(["g1"]; kcat_forward = 2.0, kcat_backward = 2.0)]
     m.reactions["r5"].gene_associations = [
         Isozyme(;
             gene_product_stoichiometry = Dict("g3" => 1, "g4" => 2),
-            kcat_forward = 3.0,
-            kcat_backward = 3.0,
+            kcat_forward = 2.0,
+            kcat_backward = 2.0,
         ),
     ]
     m.objective = Dict("r6" => 1.0)
@@ -93,18 +95,24 @@ end
 
     sm = make_simplified_enzyme_constrained_model(
         m;
-        total_gene_product_mass_bound = 0.5,
+        reaction_mass_groups = Dict("b1" => ["r3", "r4"], "b2" => ["r5", "r4"]),
+        reaction_mass_group_bounds = Dict("b2" => 0.5, "b1" => 0.2),
     )
 
-    # res = flux_balance_analysis(
-    #     gm,
-    #     Tulip.Optimizer;
-    #     modifications = [modify_optimizer_attribute("IPM_IterationsLimit", 1000)],
-    # )
+    stoichiometry(sm)
+
+    cpling = sum(coupling(sm), dims = 2)
+    @test 1.5 in cpling && 11.5 in cpling
 
-    # rxn_fluxes = values_dict(:reaction, res)
-    # gene_products = values_dict(:enzyme, res)
-    # mass_groups = values_dict(:enzyme_group, res)
+    lbs, ubs = coupling_bounds(sm)
+    @test all(lbs .== 0.0)
+    @test 0.5 in ubs && 0.2 in ubs
 
+    res = flux_balance_analysis(
+        sm,
+        Tulip.Optimizer;
+        modifications = [modify_optimizer_attribute("IPM_IterationsLimit", 1000)],
+    )
 
+    @test isapprox(solved_objective_value(res), 0.21212120975836252; atol = TEST_TOLERANCE)
 end