Improve MMDF

docs/Project.toml

@@ -9,6 +9,7 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 JSONFBCModels = "475c1105-d6ed-49c1-9b32-c11adca6d3e8"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
+UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [compat]
 Documenter = "1"

docs/src/examples/06-mmdf.jl

@@ -37,7 +37,6 @@ import Downloads: download
 # Additionally to COBREXA, and the model format package, we will need a solver
 # -- let's use GLPK here:
 
-using COBREXA
 import JSONFBCModels
 import GLPK
 
@@ -64,7 +63,6 @@ reaction_standard_gibbs_free_energies = Dict{String,Float64}(
     "TPI" => 5.621932460512994,
 )
 
-
 # (The units of the energies are kJ/mol.)
 
 # ## Running basic max min driving force analysis
@@ -85,16 +83,16 @@ reaction_standard_gibbs_free_energies = Dict{String,Float64}(
 # other reactions.
 
 reference_flux = Dict(
-    "ENO" => 1.0,
+    "ENO" => 2.0,
     "FBA" => 1.0,
-    "GAPD" => 1.0,
+    "GAPD" => 2.0,
     "GLCpts" => 1.0,
-    "LDH_D" => -1.0,
+    "LDH_D" => -2.0,
     "PFK" => 1.0,
     "PGI" => 1.0,
-    "PGK" => -1.0,
-    "PGM" => -1.0,
-    "PYK" => 1.0,
+    "PGK" => -2.0,
+    "PGM" => -2.0,
+    "PYK" => 2.0,
     "TPI" => 1.0,
 )
 
@@ -108,11 +106,12 @@ reference_flux = Dict(
 
 mmdf_solution = max_min_driving_force_analysis(
     model,
-    reaction_standard_gibbs_free_energies;
-    reference_flux,
+    reaction_standard_gibbs_free_energies,
+    reference_flux;
+    constant_concentrations = Dict("lac__D_c" => 1e-1, "nad_c" => 2.6e-3),
     concentration_ratios = Dict(
         "atp" => ("atp_c", "adp_c", 10.0),
-        "nadh" => ("nadh_c", "nad_c", 0.13),
+        "nadh" => ("nadh_c", "nad_c", 0.1),
     ),
     proton_metabolites = ["h_c", "h_e"],
     water_metabolites = ["h2o_c", "h2o_e"],
@@ -124,4 +123,32 @@ mmdf_solution = max_min_driving_force_analysis(
 )
 
 # TODO verify correctness
-@test isapprox(mmdf_solution.min_driving_force, 2.79911, atol = TEST_TOLERANCE) #src
+@test isapprox(mmdf_solution.min_driving_force, -2.4739129, atol = TEST_TOLERANCE) #src
+
+# ## Plot the results
+# We can see that the ΔG bottleneck is 2.5 kJ/mol, and that there is a
+# precipitous increase in driving force near the end of glycolysis. The overall
+# ΔG for the optimized pathway, under the restrictions in the model, is -158
+# kJ/mol, which compares favourably with the estimated ΔG under standard
+# biological conditions: -133 kJ/mol.
+
+using CairoMakie
+
+glycolysis_reaction_order =
+    ["GLCpts", "PGI", "PFK", "FBA", "TPI", "GAPD", "PGK", "PGM", "ENO", "PYK", "LDH_D"]
+
+glycolysis_thermo = cumsum(
+    reference_flux[rid] * mmdf_solution.reaction_gibbs_free_energies[Symbol(rid)] for
+    rid in glycolysis_reaction_order
+)
+
+lines(
+    1:length(glycolysis_reaction_order),
+    glycolysis_thermo,
+    axis = (
+        xlabel = "Reactions",
+        ylabel = "Cumulative ΔG [kJ/mol]",
+        xticks = (1:length(glycolysis_reaction_order), glycolysis_reaction_order),
+    ),
+)
+

src/COBREXA.jl

@@ -70,6 +70,7 @@ include("builders/loopless.jl")
 include("builders/objectives.jl")
 include("builders/scale.jl")
 include("builders/unsigned.jl")
+include("builders/mmdf.jl")
 
 # simplified front-ends for the above
 include("frontend/balance.jl")

src/builders/fbc.jl

@@ -109,46 +109,3 @@ function flux_balance_constraints(
 end
 
 export flux_balance_constraints
-
-"""
-$(TYPEDSIGNATURES)
-
-Build log-concentration-stoichiometry constraints for the `model`, as used e.g.
-by [`max_min_driving_force_analysis`](@ref).
-
-The output constraint tree contains a log-concentration variable for each
-metabolite in subtree `log_concentrations`. Individual reactions' total
-reactant log concentrations (i.e., all log concentrations of actual reactants
-minus all log concentrations of products) have their own variables in
-`reactant_log_concentrations`.
-
-Function `concentration_bound` may return a bound for the log-concentration of
-a given metabolite (compatible with `ConstraintTrees.Bound`), or `nothing`.
-"""
-function log_concentration_constraints(
-    model::A.AbstractFBCModel;
-    concentration_bound = _ -> nothing,
-)
-    rxns = Symbol.(A.reactions(model))
-    mets = Symbol.(A.metabolites(model))
-    stoi = A.stoichiometry(model)
-
-    constraints =
-        :log_concentrations^C.variables(keys = mets, bounds = concentration_bound.(mets))
-
-    cs = C.ConstraintTree()
-
-    for (midx, ridx, coeff) in zip(SparseArrays.findnz(stoi)...)
-        rid = rxns[ridx]
-        value = constraints.log_concentrations[mets[midx]].value * coeff
-        if haskey(cs, rid)
-            cs[rid].value += value
-        else
-            cs[rid] = C.Constraint(; value)
-        end
-    end
-
-    return constraints * :reactant_log_concentrations^cs
-end
-
-export log_concentration_constraints

src/builders/mmdf.jl

@@ -0,0 +1,70 @@
+
+# Copyright (c) 2021-2024, University of Luxembourg
+# Copyright (c) 2021-2024, Heinrich-Heine University Duesseldorf
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""
+$(TYPEDSIGNATURES)
+
+Allocate log-concentration-stoichiometry constraints for the `reaction_subset`
+and `metabolite_subset` in `model`, as used by e.g.
+[`max_min_driving_force_analysis`](@ref).
+
+The output constraint tree contains a log-concentration variable for each
+metabolite in subtree `log_concentrations`. The reaction quotient in log
+variables for each reaction is stored in `log_concentration_stoichiometry`.
+
+Function `concentration_bound` may return a bound for the log-concentration of a
+given metabolite (compatible with `ConstraintTrees.Bound`), or `nothing`.
+"""
+function log_concentration_constraints(
+    model::A.AbstractFBCModel;
+    reaction_subset = A.reactions(model),
+    metabolite_subset = A.metabolites(model),
+    concentration_bound = _ -> nothing,
+)
+    stoi = A.stoichiometry(model)
+
+    constraints =
+        :log_concentrations^C.variables(
+            keys = Symbol.(metabolite_subset),
+            bounds = concentration_bound.(metabolite_subset),
+        )
+
+    # map old idxs of metabolites to new order of smaller system
+    midx_lookup = Dict(
+        indexin(metabolite_subset, A.metabolites(model)) .=>
+            1:length(metabolite_subset),
+    )
+
+    cs = C.ConstraintTree()
+    for (ridx_new, ridx_old) in enumerate(indexin(reaction_subset, A.reactions(model)))
+        midxs_old, stoich_coeffs = SparseArrays.findnz(stoi[:, ridx_old])
+        rid = Symbol(reaction_subset[ridx_new])
+        for (midx_old, stoich_coeff) in zip(midxs_old, stoich_coeffs)
+            midx_new = midx_lookup[midx_old]
+            mid = Symbol(metabolite_subset[midx_new])
+            value = constraints.log_concentrations[mid].value * stoich_coeff
+            if haskey(cs, rid)
+                cs[rid].value += value
+            else
+                cs[rid] = C.Constraint(; value)
+            end
+        end
+    end
+
+    return constraints * :log_concentration_stoichiometries^cs
+end
+
+export log_concentration_constraints