Merge pull request #438 from LCSB-BioCore/mo-more-modifications

Update modifications

docs/src/notebooks/2_finding_balance.jl

@@ -70,7 +70,7 @@ dict_soln = flux_balance_analysis_dict(
         change_objective("R_BIOMASS_Ecoli_core_w_GAM"),
 
         ## this fixes a specific rate of the glucose exchange
-        change_constraint("R_EX_glc__D_e", -12, -12),
+        change_constraint("R_EX_glc__D_e"; lb = -12, ub = -12),
 
         ## this knocks out two genes, i.e. constrains their associated reactions to zero.
         knockout(["b0978", "b0734"]), ## the gene IDs are cytochrome oxidase (CYTBD)
@@ -103,8 +103,8 @@ fva_mins, fva_maxs = flux_variability_analysis_dict(
     bounds = objective_bounds(0.99), # the objective function is allowed to vary by ~1% from the FBA optimum
     modifications = [
         change_optimizer_attribute("IPM_IterationsLimit", 500),
-        change_constraint("R_EX_glc__D_e", -10, -10),
-        change_constraint("R_EX_o2_e", 0.0, 0.0),
+        change_constraint("R_EX_glc__D_e"; lb = -10, ub = -10),
+        change_constraint("R_EX_o2_e"; lb = 0.0, ub = 0.0),
     ],
 )
 
@@ -128,8 +128,8 @@ vs = flux_variability_analysis(
     bounds = objective_bounds(0.50), # biomass can vary up to 50% less than optimum
     modifications = [
         change_optimizer_attribute("IPM_IterationsLimit", 500),
-        change_constraint("R_EX_glc__D_e", -10, -10),
-        change_constraint("R_EX_o2_e", 0.0, 0.0),
+        change_constraint("R_EX_glc__D_e"; lb = -10, ub = -10),
+        change_constraint("R_EX_o2_e"; lb = 0.0, ub = 0.0),
     ],
     ret = m ->
         (COBREXA.JuMP.objective_value(m), COBREXA.JuMP.value(m[:x][biomass_idx])), # m is the model and m[:x] extracts the fluxes from the model
@@ -153,7 +153,7 @@ dict_soln = parsimonious_flux_balance_analysis_dict(
     OSQP.Optimizer;
     modifications = [
         silence, # silence the optimizer (OSQP is very verbose by default)
-        change_constraint("R_EX_glc__D_e", -12, -12),
+        change_constraint("R_EX_glc__D_e"; lb = -12, ub = -12),
     ],
 )
 
@@ -168,7 +168,7 @@ vec_soln = parsimonious_flux_balance_analysis_vec(
     model,
     Tulip.Optimizer; # start with Tulip
     modifications = [
-        change_constraint("R_EX_glc__D_e", -12, -12),
+        change_constraint("R_EX_glc__D_e"; lb = -12, ub = -12),
         change_optimizer_attribute("IPM_IterationsLimit", 500), # we may change Tulip-specific attributes here
     ],
     qp_modifications = [
@@ -176,3 +176,97 @@ vec_soln = parsimonious_flux_balance_analysis_vec(
         silence, # and make it quiet.
     ],
 )
+
+# ## Flux balance analysis with molecular crowding (FBAwMC)
+
+# Flux balance with molecular crowding incorporates enzyme capacity constraints into the
+# classic flux balance analysis algorithm. Essentially, an extra constraint is added to
+# the optimization problem: `∑ wᵢ × vᵢ ≤ 1` where `wᵢ` weights each internal flux `vᵢ`. See
+# `Beg, Qasim K., et al. "Intracellular crowding defines the mode and sequence of
+# substrate uptake by Escherichia coli and constrains its metabolic activity." Proceedings of
+# the National Academy of Sciences 104.31 (2007): 12663-12668.` for more details.
+
+# First load the model
+model = load_model("e_coli_core.json")
+
+# Next, simulate the model over a range of substrate uptake rates.
+without_crowding = Dict{Float64,Vector{Float64}}()
+with_crowding = Dict{Float64,Vector{Float64}}()
+glucose_uptakes = collect(-(1.0:0.5:20))
+
+for glc in glucose_uptakes
+    no_crowding = flux_balance_analysis_dict( # classic FBA
+        model,
+        Tulip.Optimizer;
+        modifications = [
+            change_optimizer_attribute("IPM_IterationsLimit", 1000),
+            change_constraint("EX_glc__D_e"; lb = glc),
+        ],
+    )
+
+    without_crowding[glc] =
+        [no_crowding["BIOMASS_Ecoli_core_w_GAM"], no_crowding["EX_ac_e"]]
+
+    crowding = flux_balance_analysis_dict( # FBAwMC
+        model,
+        Tulip.Optimizer;
+        modifications = [
+            change_optimizer_attribute("IPM_IterationsLimit", 1000),
+            add_crowding_constraint(0.004), # crowding constraint gets added here
+            change_constraint("EX_glc__D_e"; lb = glc),
+        ],
+    )
+
+    with_crowding[glc] = [crowding["BIOMASS_Ecoli_core_w_GAM"], crowding["EX_ac_e"]]
+end
+
+# Finally, plot the results to compare classic FBA with FBAwMC.
+using CairoMakie
+fig = Figure();
+ax1 = Axis(fig[1, 1]);
+lines!(
+    ax1,
+    -glucose_uptakes,
+    [without_crowding[glc][1] for glc in glucose_uptakes],
+    label = "no crowding",
+    linewidth = 5,
+    linestyle = :dash,
+)
+lines!(
+    ax1,
+    -glucose_uptakes,
+    [with_crowding[glc][1] for glc in glucose_uptakes],
+    label = "with crowding",
+    linewidth = 5,
+    linestyle = :dot,
+)
+ax1.ylabel = "Growth rate [1/h]"
+ax2 = Axis(fig[2, 1])
+lines!(
+    ax2,
+    -glucose_uptakes,
+    [without_crowding[glc][2] for glc in glucose_uptakes],
+    label = "no crowding",
+    linewidth = 5,
+    linestyle = :dash,
+)
+lines!(
+    ax2,
+    -glucose_uptakes,
+    [with_crowding[glc][2] for glc in glucose_uptakes],
+    label = "with crowding",
+    linewidth = 5,
+    linestyle = :dot,
+)
+fig[1:2, 2] = Legend(fig, ax1, "Constraint")
+ax2.xlabel = "Glucose uptake rate [mmol/gDW/h]"
+ax2.ylabel = "Acetate production rate\n[mmol/gDW/h]"
+fig
+
+# Notice that overflow metabolism is modeled by incorporating the crowding constraint.
+
+#md # !!! tip "Tip: code your own modification like [`add_crowding`](@ref)"
+#md #       Making custom problem modification functions is really simple due to the
+#md #       tight intergration between COBREXA and JuMP. Look at the source code for
+#md #       the implemented modifications in `src\analysis\modifications` to get a flavour
+#md #       for it.

docs/src/notebooks/3_basic_stdmodel_usage.jl

@@ -46,8 +46,8 @@ fluxes = flux_balance_analysis_dict(
     Tulip.Optimizer;
     modifications = [
         change_objective("BIOMASS_Ecoli_core_w_GAM"),
-        change_constraint("EX_glc__D_e", -12, -12),
-        change_constraint("EX_o2_e", 0, 0),
+        change_constraint("EX_glc__D_e"; lb = -12, ub = -12),
+        change_constraint("EX_o2_e"; lb = 0, ub = 0),
     ],
 )
 

docs/src/notebooks/4_basic_core_coupled_usage.jl

@@ -32,8 +32,8 @@ dict_sol = flux_balance_analysis_dict(
     Tulip.Optimizer;
     modifications = [
         change_objective("R_BIOMASS_Ecoli_core_w_GAM"),
-        change_constraint("R_EX_glc__D_e", -12, -12),
-        change_constraint("R_EX_o2_e", 0, 0),
+        change_constraint("R_EX_glc__D_e"; lb = -12, ub = -12),
+        change_constraint("R_EX_o2_e"; lb = 0, ub = 0),
     ],
 )
 

src/analysis/flux_variability_analysis.jl

@@ -125,8 +125,8 @@ mins, maxs = flux_variability_analysis_dict(
     bounds = objective_bounds(0.99),
     modifications = [
         change_optimizer_attribute("IPM_IterationsLimit", 500),
-        change_constraint("EX_glc__D_e", -10, -10),
-        change_constraint("EX_o2_e", 0.0, 0.0),
+        change_constraint("EX_glc__D_e"; lb = -10, ub = -10),
+        change_constraint("EX_o2_e"; lb = 0, ub = 0),
     ],
 )
 ```

src/analysis/modifications/crowding.jl

+
+"""
+    add_crowding_constraint(weights::Dict{Int64, Float64})
+
+Adds a molecular crowding constraint to the optimization problem: `∑ wᵢ × vᵢ ≤ 1` where `wᵢ`
+is a weight and `vᵢ` is a flux index in the model's reactions specified in `weights` as `vᵢ
+=> wᵢ` pairs.
+
+See Beg, Qasim K., et al. "Intracellular crowding defines the mode and sequence of
+substrate uptake by Escherichia coli and constrains its metabolic activity." Proceedings of
+the National Academy of Sciences 104.31 (2007) for more details.
+"""
+add_crowding_constraint(weights::Dict{Int64,Float64}) =
+    (model, opt_model) -> begin
+        idxs = collect(keys(weights)) # order of keys and values is the same
+        ws = values(weights)
+        # since fluxes can be positive or negative, need absolute value: ∑ wᵢ × |vᵢ| ≤ 1
+        # introduce slack variables to handle this
+        @variable(opt_model, crowding_slack[1:length(weights)])
+        @constraint(opt_model, crowding_slack .>= opt_model[:x][idxs])
+        @constraint(opt_model, crowding_slack .>= -opt_model[:x][idxs])
+        @constraint(opt_model, sum(w * crowding_slack[i] for (i, w) in enumerate(ws)) <= 1)
+    end
+
+"""
+    add_crowding_constraint(weight::Float64; kwargs)
+
+Variant of [`add_crowding_constraint`](@ref) that takes a single weight and assigns it to
+each internal reaction flux, where internal reactions are identified with
+[`find_internal_reactions`](@ref) and `kwargs` are passed to this function.
+"""
+add_crowding_constraint(weight::Float64; kwargs...) =
+    (model, opt_model) -> begin
+        idxs = find_internal_reactions(model; kwargs...)
+        add_crowding_constraint(Dict(zip(idxs, fill(weight, length(idxs)))))(
+            model,
+            opt_model,
+        )
+    end
+
+"""
+    add_crowding_constraint(weights::Dict{String, Float64})
+
+Variant of [`add_crowding_constraint`](@ref) that takes a dictinary of reactions `ids`
+instead of reaction indices mapped to weights.
+"""
+add_crowding_constraint(weights::Dict{String,Float64}) =
+    (model, opt_model) -> begin
+        idxs = indexin(keys(weights), reactions(model))
+        nothing in idxs && throw(ArgumentError("Reaction id not found in model."))
+        add_crowding_constraint(Dict(zip(Int.(idxs), values(weights))))(model, opt_model)
+    end

src/analysis/modifications/generic.jl

@@ -13,11 +13,11 @@ constrain_objective_value(tolerance) =
     end
 
 """
-    change_constraint(id::String, lb, ub)
+    change_constraint(id::String; lb=nothing, ub=nothing)
 
-Change the lower and upper bounds (`lb` and `ub` respectively) of reaction `id`.
+Change the lower and upper bounds (`lb` and `ub` respectively) of reaction `id` if supplied.
 """
-change_constraint(id::String, lb, ub) =
+change_constraint(id::String; lb = nothing, ub = nothing) =
     (model, opt_model) -> begin
         ind = first(indexin([id], reactions(model)))
         isnothing(ind) && throw(DomainError(id, "No matching reaction was found."))

src/base/utils/looks_like.jl

@@ -79,6 +79,7 @@ function looks_like_biomass_reaction(
     any(occursin(x, rxn_id) for x in biomass_strings) &&
         !(exclude_exchanges && any(startswith(rxn_id, x) for x in exchange_prefixes))
 end
+
 """
     find_biomass_reactions(m::MetabolicModel; kwargs...)
 
@@ -97,6 +98,66 @@ forwarded to [`looks_like_biomass_reaction`](@ref).
 find_biomass_reaction_ids(m::MetabolicModel; kwargs...) =
     filter(id -> looks_like_biomass_reaction(id, kwargs...), reactions(m))
 
+"""
+    looks_like_internal_reaction(
+        rxn_id::String;
+        exchange_prefixes = _constants.exchange_prefixes,
+        biomass_strings = _constants.biomass_strings,
+    )::Bool
+
+A predicate that matches reaction identifiers that look like
+internal reactions, i.e. reactions that are neither exchange nor biomass reactions.
+Exchange reactions are identified based on matching prefixes in
+the set `exchange_prefixes` and biomass reactions are identified by looking for
+occurences of `biomass_strings` in the reaction id.
+
+Also see [`find_internal_reactions`](@ref).
+
+# Example
+```
+findall(looks_like_internal_reaction, reactions(model)) # returns indices
+filter(looks_like_internal_reaction, reactions(model)) # returns Strings
+
+# to use the optional arguments you need to expand the function's arguments
+# using an anonymous function
+findall(x -> looks_like_internal_reaction(x; exchange_prefixes=["EX_", "R_EX_"]), reactions(model)) # returns indices
+filter(x -> looks_like_internal_reaction(x; exchange_prefixes=["EX_", "R_EX_"]), reactions(model)) # returns Strings
+```
+"""
+function looks_like_internal_reaction(
+    rxn_id::String;
+    exchange_prefixes = _constants.exchange_prefixes,
+    biomass_strings = _constants.biomass_strings,
+)::Bool
+    !looks_like_biomass_reaction(
+        rxn_id;
+        exchange_prefixes = exchange_prefixes,
+        biomass_strings = biomass_strings,
+    ) &&
+        !looks_like_exchange_reaction(
+            rxn_id;
+            exchange_prefixes = exchange_prefixes,
+            biomass_strings = biomass_strings,
+        )
+end
+
+"""
+    find_internal_reactions(m::MetabolicModel; kwargs...)
+
+Shortcut for finding internal reaction indices in a model; arguments are
+forwarded to [`looks_like_internal_reaction`](@ref).
+"""
+find_internal_reactions(m::MetabolicModel; kwargs...) =
+    findall(id -> looks_like_internal_reaction(id; kwargs...), reactions(m))
+
+"""
+    find_internal_reaction_ids(m::MetabolicModel; kwargs...)
+
+Shortcut for finding internal reaction identifiers in a model; arguments are
+forwarded to [`looks_like_internal_reaction`](@ref).
+"""
+find_internal_reaction_ids(m::MetabolicModel; kwargs...) =
+    filter(id -> looks_like_internal_reaction(id, kwargs...), reactions(m))
 
 """
     looks_like_exchange_metabolite(rxn_id::String;

test/analysis/fba_variants.jl

+@testset "FBAwMC" begin
+    model = load_model(model_paths["e_coli_core.json"])
+    idxs = find_internal_reactions(model)
+
+    sol = flux_balance_analysis_dict(
+        model,
+        Tulip.Optimizer;
+        modifications = [
+            change_optimizer_attribute("IPM_IterationsLimit", 1000),
+            add_crowding_constraint(0.004),
+            change_constraint("EX_glc__D_e"; lb = -6),
+        ],
+    )
+
+    @test isapprox(
+        sol["BIOMASS_Ecoli_core_w_GAM"],
+        0.491026987015203,
+        atol = TEST_TOLERANCE,
+    )
+
+    @test isapprox(sol["EX_ac_e"], 0.7084745257320869, atol = TEST_TOLERANCE)
+end

test/analysis/flux_balance_analysis.jl

@@ -39,7 +39,7 @@ end
         Tulip.Optimizer;
         modifications = [
             change_objective("BIOMASS_Ecoli_core_w_GAM"),
-            change_constraint("EX_glc__D_e", -12, -12),
+            change_constraint("EX_glc__D_e"; lb = -12, ub = -12),
             change_sense(MAX_SENSE),
             change_optimizer_attribute("IPM_IterationsLimit", 110),
         ],
@@ -71,7 +71,7 @@ end
     @test_throws DomainError flux_balance_analysis_dict(
         model,
         Tulip.Optimizer;
-        modifications = [change_constraint("gbbrsh", -12, -12)],
+        modifications = [change_constraint("gbbrsh"; lb = -12, ub = -12)],
     )
     @test_throws DomainError flux_balance_analysis_dict(
         model,

test/analysis/flux_variability_analysis.jl

@@ -82,8 +82,8 @@ end
         bounds = objective_bounds(0.99),
         modifications = [
             change_optimizer_attribute("IPM_IterationsLimit", 500),
-            change_constraint("EX_glc__D_e", -10, -10),
-            change_constraint("EX_o2_e", 0.0, 0.0),
+            change_constraint("EX_glc__D_e"; lb = -10, ub = -10),
+            change_constraint("EX_o2_e"; lb = 0.0, ub = 0.0),
         ],
     )
 

test/analysis/knockouts.jl

@@ -104,7 +104,7 @@ end
             Tulip.Optimizer;
             modifications = [
                 change_objective("BIOMASS_Ecoli_core_w_GAM"),
-                change_constraint("EX_glc__D_e", -12, -12),
+                change_constraint("EX_glc__D_e"; lb = -12, ub = -12),
                 change_sense(MAX_SENSE),
                 change_optimizer_attribute("IPM_IterationsLimit", 110),
                 knockout(["b0978", "b0734"]), # knockouts out cytbd
@@ -121,7 +121,7 @@ end
             Tulip.Optimizer;
             modifications = [
                 change_objective("BIOMASS_Ecoli_core_w_GAM"),
-                change_constraint("EX_glc__D_e", -12, -12),
+                change_constraint("EX_glc__D_e"; lb = -12, ub = -12),
                 change_sense(MAX_SENSE),
                 change_optimizer_attribute("IPM_IterationsLimit", 110),
                 knockout("b2779"), # knockouts out enolase

test/analysis/parsimonious_flux_balance_analysis.jl

@@ -4,7 +4,7 @@
         model,
         Tulip.Optimizer;
         modifications = [
-            change_constraint("EX_glc__D_e", -12, -12),
+            change_constraint("EX_glc__D_e"; lb = -12, ub = -12),
             change_optimizer_attribute("IPM_IterationsLimit", 500),
         ],
         qp_modifications = [change_optimizer(OSQP.Optimizer), silence],
@@ -13,7 +13,7 @@
         model,
         Tulip.Optimizer;
         modifications = [
-            change_constraint("EX_glc__D_e", -12, -12),
+            change_constraint("EX_glc__D_e"; lb = -12, ub = -12),
             change_optimizer_attribute("IPM_IterationsLimit", 500),
         ],
         qp_modifications = [change_optimizer(OSQP.Optimizer), silence],

test/analysis/sampling/warmup_variability.jl

@@ -5,7 +5,7 @@
         model,
         Tulip.Optimizer,
         100;
-        modifications = [change_constraint("EX_glc__D_e", -2, 2)],
+        modifications = [change_constraint("EX_glc__D_e"; lb = -2, ub = 2)],
         workers = W,
     )
 

test/base/utils/looks_like.jl

@@ -15,6 +15,8 @@
         ["m1"; "m2"; "m3"],
     )
     @test filter(looks_like_exchange_reaction, reactions(cp)) == ["EX_m1"]
+    @test find_internal_reactions(cp) == [2, 3]
+    @test find_internal_reaction_ids(cp) == ["r2", "r3"]
 
     cp = CoreModel(
         [-1.0 0 0; 0 0 -1; 0 -1 0],
@@ -55,25 +57,30 @@ end
     @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20
     @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20
     @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1
+    @test length(filter(looks_like_internal_reaction, reactions(model))) == 74
 
     model = load_model(model_paths["e_coli_core.xml"])
     @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20
     @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20
     @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1
+    @test length(filter(looks_like_internal_reaction, reactions(model))) == 74
 
     model = load_model(model_paths["e_coli_core.mat"])
     @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20
     @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20
     @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1
+    @test length(filter(looks_like_internal_reaction, reactions(model))) == 74
 
     model = convert(StandardModel, model)
     @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20
     @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20
     @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1
+    @test length(filter(looks_like_internal_reaction, reactions(model))) == 74
 
     model = convert(CoreModelCoupled, model)
     @test length(filter(looks_like_exchange_reaction, reactions(model))) == 20
     @test length(filter(looks_like_exchange_metabolite, metabolites(model))) == 20
     @test length(filter(looks_like_biomass_reaction, reactions(model))) == 1
+    @test length(filter(looks_like_internal_reaction, reactions(model))) == 74
 
 end