Thermodynamic models
using COBREXAHere we will solve the max min driving force analysis problem using the glycolysis pathway of E. coli. In essence, the method attempts to find metabolite concentrations (NB: not fluxes) that maximize the smallest thermodynamic driving force through each reaction. See Noor, et al., "Pathway thermodynamics highlights kinetic obstacles in central metabolism.", PLoS computational biology, 2014, for more details.
To do this, we will first need a model that includes glycolysis, which we can download if it is not already present.
using COBREXA
download_model(
"http://bigg.ucsd.edu/static/models/e_coli_core.json",
"e_coli_core.json",
"7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
)"e_coli_core.json"Additionally to COBREXA, and the model format package, we will need a solver – let's use HiGHS here:
import JSONFBCModels
import HiGHS
model = load_model("e_coli_core.json")JSONFBCModels.JSONFBCModel(#= 95 reactions, 72 metabolites =#)Thermodynamic data
We will need ΔᵣG⁰ data for each reaction we want to include in the thermodynamic model. To generate this data manually, use eQuilibrator. To generate automatically, it is possible to use the eQuilibrator.jl package.
reaction_standard_gibbs_free_energies = Dict{String,Float64}( # units of the energies are kJ/mol
"ENO" => -3.8108376097261782,
"FBA" => 23.376920310319235,
"GAPD" => 0.5307809794271634,
"GLCpts" => -45.42430981510088,
"LDH_D" => 20.04059765689044,
"PFK" => -18.546314942995934,
"PGI" => 2.6307087407442395,
"PGK" => 19.57192102020454,
"PGM" => -4.470553692565886,
"PYK" => -24.48733600711958,
"TPI" => 5.621932460512994,
)Dict{String, Float64} with 11 entries:
"PFK" => -18.5463
"PGI" => 2.63071
"GLCpts" => -45.4243
"PGK" => 19.5719
"TPI" => 5.62193
"LDH_D" => 20.0406
"ENO" => -3.81084
"PYK" => -24.4873
"FBA" => 23.3769
"GAPD" => 0.530781
"PGM" => -4.47055Running basic max min driving force analysis
If a reference flux is not specified, it is assumed that every reaction in the model should be included in the thermodynamic model, and that each reaction proceeds in the forward direction. This is usually not intended, and can be prevented by inputting a reference flux dictionary as shown below. This dictionary can be a flux solution. The sign of each flux is used to determine if the reaction runs forward or backward.
Using a reference solution
Frequently it is useful to check the max-min driving force of a specific FBA solution. In this case, one is usually only interested in a subset of all the reactions in a model. These reactions can be specified as a the reference_flux, to only compute the MMDF of these reactions, and ignore all other reactions.
reference_flux = Dict(
"ENO" => 1.0,
"FBA" => 1.0,
"GAPD" => 1.0,
"GLCpts" => 1.0,
"LDH_D" => -1.0,
"PFK" => 1.0,
"PGI" => 1.0,
"PGK" => -1.0,
"PGM" => 0.0,
"PYK" => 1.0,
"TPI" => 1.0,
)Dict{String, Float64} with 11 entries:
"PFK" => 1.0
"PGI" => 1.0
"GLCpts" => 1.0
"PGK" => -1.0
"TPI" => 1.0
"LDH_D" => -1.0
"ENO" => 1.0
"PYK" => 1.0
"FBA" => 1.0
"GAPD" => 1.0
"PGM" => 0.0It is most convenient to pass a flux solution into reference_flux, but take care about the fluxes with value near 0: Their desired sign may be a subject to floating-point robustness error. By default, max_min_driving_force_analysis considers everything that is approximately zero (via isapprox) to have zero flux, with the appropriate implications to concentration balance.
Solving the MMDF problem
mmdf_solution = max_min_driving_force_analysis(
model;
reaction_standard_gibbs_free_energies,
reference_flux,
constant_concentrations = Dict("g3p_c" => exp(-8.5)),
concentration_ratios = Dict(
"atp" => ("atp_c", "adp_c", 10.0),
"nadh" => ("nadh_c", "nad_c", 0.13),
),
proton_metabolites = ["h_c"],
water_metabolites = ["h2o_c"],
concentration_lower_bound = 1e-6, # mol/L
concentration_upper_bound = 1e-1, # mol/L
T = 298.15, # Kelvin
R = 8.31446261815324e-3, # kJ/K/mol
optimizer = HiGHS.Optimizer,
)ConstraintTrees.Tree{Float64} with 6 elements:
:concentration_ratio_constraints => ConstraintTrees.Tree{Float64}(#= 2 elemen…
:driving_forces => ConstraintTrees.Tree{Float64}(#= 11 eleme…
:log_concentration_stoichiometry => ConstraintTrees.Tree{Float64}(#= 11 eleme…
:log_concentrations => ConstraintTrees.Tree{Float64}(#= 17 eleme…
:min_driving_force => 1.88051
:min_driving_force_thresholds => ConstraintTrees.Tree{Float64}(#= 11 eleme…One may be also interested in seeing the FVA-like feasible concentration ranges in such model. The most straightforward way to find these is to use the associated constraint-system-building function max_min_driving_force_constraints together with constraints_variability as follows:
mmdf_system = max_min_driving_force_constraints(
model;
reaction_standard_gibbs_free_energies,
reference_flux,
constant_concentrations = Dict("g3p_c" => exp(-8.5)),
concentration_ratios = Dict(
"atp" => ("atp_c", "adp_c", 10.0),
"nadh" => ("nadh_c", "nad_c", 0.13),
),
proton_metabolites = ["h_c"],
water_metabolites = ["h2o_c"],
concentration_lower_bound = 1e-6, # mol/L
concentration_upper_bound = 1e-1, # mol/L
T = 298.15, # Kelvin
R = 8.31446261815324e-3, # kJ/K/mol
)
cva_solution = constraints_variability(
mmdf_system,
mmdf_system.log_concentrations,
objective = mmdf_system.min_driving_force.value,
optimizer = HiGHS.Optimizer,
)ConstraintTrees.Tree{Tuple{Union{Nothing, Float64}, Union{Nothing, Float64}}} with 17 elements:
Symbol("13dpg_c") => (-13.8155, -13.0569)
Symbol("2pg_c") => (-13.8155, -4.66251)
Symbol("3pg_c") => (-12.0121, -2.85911)
:adp_c => (-11.5129, -2.30259)
:atp_c => (-13.8155, -4.60517)
:dhap_c => (-6.23214, -3.23273)
:f6p_c => (-10.4809, -3.3638)
:fdp_c => (-5.30199, -2.30259)
:g3p_c => (-8.5, -8.5)
:g6p_c => (-9.41969, -2.30259)
:glc__D_e => (-13.8155, -2.30259)
:lac__D_c => (-13.8155, -2.30259)
:nad_c => (-13.8155, -4.34281)
:nadh_c => (-11.7753, -2.30259)
:pep_c => (-13.8155, -3.12524)
:pi_c => (-3.06118, -2.30259)
:pyr_c => (-13.8155, -2.30259)This page was generated using Literate.jl.