diff --git a/docs/Project.toml b/docs/Project.toml index 2661856c..ff981410 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/src/examples/06-mmdf.jl b/docs/src/examples/06-mmdf.jl index df3e672c..f81b37ae 100644 --- a/docs/src/examples/06-mmdf.jl +++ b/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), + ), +) + diff --git a/src/COBREXA.jl b/src/COBREXA.jl index 9719c4b5..a9754c58 100644 --- a/src/COBREXA.jl +++ b/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") diff --git a/src/builders/fbc.jl b/src/builders/fbc.jl index 910895a8..50750ca7 100644 --- a/src/builders/fbc.jl +++ b/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 diff --git a/src/builders/mmdf.jl b/src/builders/mmdf.jl new file mode 100644 index 00000000..26e7b531 --- /dev/null +++ b/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 diff --git a/src/frontend/mmdf.jl b/src/frontend/mmdf.jl index 27cf4639..84730046 100644 --- a/src/frontend/mmdf.jl +++ b/src/frontend/mmdf.jl @@ -18,9 +18,10 @@ $(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 @@ -28,13 +29,13 @@ 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, )