change mmdf structure a little bit

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
@@ -108,8 +106,8 @@ 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;
     concentration_ratios = Dict(
         "atp" => ("atp_c", "adp_c", 10.0),
         "nadh" => ("nadh_c", "nad_c", 0.13),
@@ -124,4 +122,4 @@ 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.79911, atol = TEST_TOLERANCE) #src

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

src/frontend/mmdf.jl

@@ -18,23 +18,24 @@
 $(TYPEDSIGNATURES)
 
 Perform a max-min driving force analysis using `optimizer` on the `model` with
-supplied reaction standard Gibbs energies in `reaction_standard_gibbs_free_energies`.
+supplied reaction standard Gibbs energies in
+`reaction_standard_gibbs_free_energies`.
 
-The method was described by by: Noor, et al., "Pathway thermodynamics highlights
+The method is described by by: Noor, et al., "Pathway thermodynamics highlights
 kinetic obstacles in central metabolism.", PLoS computational biology, 2014.
 
 `reference_flux` sets the directions of each reaction in `model`. The scale of
 the values is not important, only the direction is examined (w.r.t.
 `reference_flux_atol` tolerance). Ideally, the supplied `reference_flux` should
 be completely free of internal cycles, which enables the thermodynamic
 consistency. To get the cycle-free flux, you can use
-[`loopless_flux_balance_analysis`](@ref) (computationally complex but gives
-very good fluxes), [`parsimonious_flux_balance_analysis`](@ref) or
+[`loopless_flux_balance_analysis`](@ref) (computationally expensive, but
+consistent), [`parsimonious_flux_balance_analysis`](@ref) or
 [`linear_parsimonious_flux_balance_analysis`](@ref) (computationally simplest
 but the consistency is not guaranteed).
 
-Internally, [`log_concentration_constraints`](@ref) is used to lay out the
-base structure of the problem.
+Internally, [`log_concentration_constraints`](@ref) is used to lay out the base
+structure of the problem.
 
 Following arguments are set optionally:
 - `water_metabolites`, `proton_metabolites` and `ignored_metabolites` allow to
@@ -46,8 +47,8 @@ Following arguments are set optionally:
 - `concentration_lower_bound` and `concentration_upper_bound` set the default
   concentration bounds for all other metabolites
 - `concentration ratios` is a dictionary that assigns a tuple of
-  metabolite-metabolite-concentration ratio constraint names; e.g. ATP/ADP
-  ratio can be fixed to five-times-more-ATP by setting `concentration_ratios =
+  metabolite-metabolite-concentration ratio constraint names; e.g. ATP/ADP ratio
+  can be fixed to five-times-more-ATP by setting `concentration_ratios =
   Dict("adenosin_ratio" => ("atp", "adp", 5.0))`
 - `T` and `R` default to the usual required thermodynamic constraints in the
   usual units (K and kJ/K/mol, respectively). These multiply the
@@ -56,93 +57,72 @@ Following arguments are set optionally:
 """
 function max_min_driving_force_analysis(
     model::A.AbstractFBCModel,
-    reaction_standard_gibbs_free_energies::Dict{String,Float64};
-    reference_flux = Dict{String,Float64}(),
+    reaction_standard_gibbs_free_energies::Dict{String,Float64},
+    reference_flux = Dict{String,Float64}();
     concentration_ratios = Dict{String,Tuple{String,String,Float64}}(),
     constant_concentrations = Dict{String,Float64}(),
     ignored_metabolites = [],
     proton_metabolites = [],
     water_metabolites = [],
-    concentration_lower_bound = 1e-9, # M
-    concentration_upper_bound = 1e-1, # M
+    concentration_lower_bound = 1e-9, # Molar
+    concentration_upper_bound = 1e-1, # Molar
     T = 298.15, # Kelvin
     R = 8.31446261815324e-3, # kJ/K/mol
     reference_flux_atol = 1e-6,
-    check_ignored_reactions = missing,
     settings = [],
     optimizer,
 )
 
+    #=
+    For MMDF, one is almost always interested in running the analysis on a
+    subset of the full model, because thermodynamic data is usually not
+    available for all the reactions/metabolites. Missing data may cause subtle
+    problems, and is usually best to restrict your model to reactions where
+    thermodynamic data exists.  
+    =#
+
+    reaction_subset =
+        isempty(reference_flux) ? A.reactions(model) :
+        collect(k for (k, v) in reference_flux if abs(v) > reference_flux_atol)
+
+    metabolite_subset =
+        isempty(reference_flux) ? A.metabolites(model) :
+        unique(
+            mid for rid in reaction_subset for
+            mid in keys(A.reaction_stoichiometry(model, string(rid)))
+        )
+
+
     # First let's just check if all the identifiers are okay because we use
     # quite a lot of these; the rest of the function may be a bit cleaner with
     # this checked properly.
 
-    model_reactions = Set(A.reactions(model))
-    model_metabolites = Set(A.metabolites(model))
-
-    all(in(model_reactions), keys(reaction_standard_gibbs_free_energies)) || throw(
+    all(in(Set(keys(reaction_standard_gibbs_free_energies))), reaction_subset) || throw(
         DomainError(
             reaction_standard_gibbs_free_energies,
-            "unknown reactions referenced by reaction_standard_gibbs_free_energies",
+            "some reactions are not accounted for in reaction_standard_gibbs_free_energies.",
         ),
     )
-    all(x -> haskey(reaction_standard_gibbs_free_energies, x), keys(reference_flux)) ||
-        throw(DomainError(reference_flux, "some reactions have no reference flux"))
-    all(in(model_reactions), keys(reference_flux)) || throw(
-        DomainError(
-            reaction_standard_gibbs_free_energies,
-            "unknown reactions referenced by reference_flux",
-        ),
+    all(in(Set(A.reactions(model))), reaction_subset) || throw(
+        DomainError(reference_flux, "unknown reactions referenced by reference_flux."),
     )
-    all(in(model_metabolites), keys(constant_concentrations)) || throw(
+    all(in(metabolite_subset), keys(constant_concentrations)) || throw(
         DomainError(
             constant_concentrations,
-            "unknown metabolites referenced by constant_concentrations",
+            "unknown metabolites referenced by constant_concentrations.",
         ),
     )
     all(
-        in(model_metabolites),
+        in(metabolite_subset),
         (m for (_, (x, y, _)) in concentration_ratios for m in (x, y)),
     ) || throw(
         DomainError(
             concentration_ratios,
-            "unknown metabolites referenced by concentration_ratios",
-        ),
-    )
-    all(in(model_metabolites), proton_metabolites) || throw(
-        DomainError(
-            concentration_ratios,
-            "unknown metabolites referenced by proton_metabolites",
-        ),
-    )
-    all(in(model_metabolites), water_metabolites) || throw(
-        DomainError(
-            concentration_ratios,
-            "unknown metabolites referenced by water_metabolites",
+            "unknown metabolites referenced by concentration_ratios.",
         ),
     )
-    all(in(model_metabolites), ignored_metabolites) || throw(
-        DomainError(
-            concentration_ratios,
-            "unknown metabolites referenced by ignored_metabolites",
-        ),
-    )
-    if !ismissing(check_ignored_reactions) && (
-        all(
-            x -> !haskey(reaction_standard_gibbs_free_energies, x),
-            check_ignored_reactions,
-        ) || (
-            union(
-                Set(check_ignored_reactions),
-                Set(keys(reaction_standard_gibbs_free_energies)),
-            ) != model_reactions
-        )
-    )
-        throw(AssertionError("check_ignored_reactions validation failed"))
-    end
-
-    # that was a lot of checking.
 
+    # setup allocations
     default_concentration_bound =
         C.Between(log(concentration_lower_bound), log(concentration_upper_bound))
 
@@ -152,13 +132,16 @@ function max_min_driving_force_analysis(
         Set(Symbol.(ignored_metabolites)),
     )
 
+    # allocate variables
     constraints =
         log_concentration_constraints(
-            model,
-            concentration_bound = met -> if met in no_concentration_metabolites
+            model;
+            reaction_subset,
+            metabolite_subset,
+            concentration_bound = met -> if Symbol(met) in no_concentration_metabolites
                 C.EqualTo(0)
             else
-                mid = String(met)
+                mid = met
                 if haskey(constant_concentrations, mid)
                     C.EqualTo(log(constant_concentrations[mid]))
                 else
@@ -167,25 +150,30 @@ function max_min_driving_force_analysis(
             end,
         ) + :min_driving_force^C.variable()
 
+    # constrain variables
     driving_forces = C.ConstraintTree(
+        let r = Symbol(rid),
+            dGr0 = reaction_standard_gibbs_free_energies[rid],
+            dGr = dGr0 + R * T * constraints.log_concentration_stoichiometries[r].value
+
+            r => C.Constraint(dGr, C.Between(-Inf, Inf))
+
+        end for rid in reaction_subset
+    )
+
+    driving_force_thresholds = C.ConstraintTree(
         let r = Symbol(rid),
             rf = reference_flux[rid],
-            df = dGr0 + R * T * constraints.reactant_log_concentrations[r].value
+            dGr = (rf > 0 ? 1 : -1) * driving_forces[r].value
 
-            r => if isapprox(rf, 0.0, atol = reference_flux_atol)
-                C.Constraint(df, C.EqualTo(0))
-            else
-                C.Constraint(rf > 0 ? -df : df, C.Between(0, Inf))
-            end
-        end for (rid, dGr0) in reaction_standard_gibbs_free_energies
+            r => less_or_equal_constraint(dGr, constraints.min_driving_force)
+        end for rid in reaction_subset
     )
 
     constraints =
         constraints *
-        :driving_forces^driving_forces *
-        :min_driving_force_thresholds^C.map(driving_forces) do c
-            less_or_equal_constraint(constraints.min_driving_force, c)
-        end *
+        :reaction_gibbs_free_energies^driving_forces *
+        :min_driving_force_thresholds^driving_force_thresholds *
         :concentration_ratio_constraints^C.ConstraintTree(
             Symbol(cid) => difference_constraint(
                 constraints.log_concentrations[Symbol(m1)],
@@ -197,6 +185,7 @@ function max_min_driving_force_analysis(
     optimized_values(
         constraints;
         objective = constraints.min_driving_force.value,
+        sense = Minimal,
         optimizer,
         settings,
     )