Homogenize names of functions

docs/src/examples/02-flux-balance-analysis.jl

@@ -28,7 +28,7 @@ model = load_model("e_coli_core.json")
 # is captured in the default behavior of function
 # [`flux_balance`](@ref):
 
-solution = flux_balance(model, Tulip.Optimizer)
+solution = flux_balance_analysis(model, Tulip.Optimizer)
 
 @test isapprox(solution.objective, 0.8739, atol = TEST_TOLERANCE) #src
 

docs/src/examples/02a-optimizer-parameters.jl

@@ -29,7 +29,7 @@ model = load_model("e_coli_core.json")
 
 # Running a FBA with a silent optimizer that has slightly increased iteration
 # limit for IPM algorithm may now look as follows:
-solution = flux_balance(
+solution = flux_balance_analysis(
     model,
     Tulip.Optimizer;
     modifications = [silence, set_optimizer_attribute("IPM_IterationsLimit", 1000)],
@@ -42,7 +42,7 @@ solution = flux_balance(
 # will cause it to fail, return no solution, and verbosely describe what
 # happened:
 
-solution = flux_balance(
+solution = flux_balance_analysis(
     model,
     Tulip.Optimizer;
     modifications = [set_optimizer_attribute("IPM_IterationsLimit", 2)],

docs/src/examples/02b-model-modifications.jl

@@ -63,7 +63,7 @@ model.reactions["CS"].stoichiometry
 
 import GLPK
 
-base_solution = flux_balance(model, GLPK.Optimizer)
+base_solution = flux_balance_analysis(model, GLPK.Optimizer)
 base_solution.objective
 
 # Now, for example, we can limit the intake of glucose by the model:
@@ -76,7 +76,7 @@ model.reactions["EX_glc__D_e"].lower_bound = -5.0
 
 # ...and solve the modified model:
 #
-low_glucose_solution = flux_balance(model, GLPK.Optimizer)
+low_glucose_solution = flux_balance_analysis(model, GLPK.Optimizer)
 low_glucose_solution.objective
 
 @test isapprox(low_glucose_solution.objective, 0.41559777, atol = TEST_TOLERANCE) #src

docs/src/examples/02c-constraint-modifications.jl

@@ -36,7 +36,7 @@ model = load_model("e_coli_core.json")
 
 import ConstraintTrees as C
 
-ctmodel = fbc_model_constraints(model)
+ctmodel = build_flux_balance_model(model)
 
 fermentation = ctmodel.fluxes.EX_ac_e.value + ctmodel.fluxes.EX_etoh_e.value
 

docs/src/examples/03-parsimonious-flux-balance.jl

@@ -24,11 +24,12 @@ model = load_model("e_coli_core.json") # load the model
 
 # Use the convenience function to run standard pFBA on
 
-vt = parsimonious_flux_balance(model, Clarabel.Optimizer; modifications = [silence])
+vt =
+    parsimonious_flux_balance_analysis(model, Clarabel.Optimizer; modifications = [silence])
 
 # Or use the piping functionality
 
-model |> parsimonious_flux_balance(Clarabel.Optimizer; modifications = [silence])
+model |> parsimonious_flux_balance_analysis(Clarabel.Optimizer; modifications = [silence])
 
 @test isapprox(vt.objective, 0.87392; atol = TEST_TOLERANCE) #src
 @test sum(x^2 for x in values(vt.fluxes)) < 15000 #src
@@ -38,7 +39,7 @@ model |> parsimonious_flux_balance(Clarabel.Optimizer; modifications = [silence]
 # Alternatively, you can construct your own constraint tree model with
 # the quadratic objective (this approach is much more flexible).
 
-ctmodel = fbc_model_constraints(model)
+ctmodel = build_flux_balance_model(model)
 ctmodel *= :l2objective^squared_sum_objective(ctmodel.fluxes)
 ctmodel.objective.bound = 0.3 # set growth rate # TODO currently breaks
 
@@ -73,7 +74,7 @@ minimize_metabolic_adjustment(ref_sol, Clarabel.Optimizer; modifications = [sile
 # Alternatively, you can construct your own constraint tree model with
 # the quadratic objective (this approach is much more flexible).
 
-ctmodel = fbc_model_constraints(model)
+ctmodel = build_flux_balance_model(model)
 ctmodel *=
     :minoxphospho^squared_sum_error_objective(
         ctmodel.fluxes,

docs/src/examples/04-minimization-of-metabolic-adjustment.jl

@@ -18,6 +18,6 @@ import Clarabel
 # guessing.
 model = convert(CM.Model, load_model("e_coli_core.json"))
 
-reference_fluxes = parsimonious_flux_balance(model, Clarabel.Optimizer).fluxes
+reference_fluxes = parsimonious_flux_balance_analysis(model, Clarabel.Optimizer).fluxes
 
 # TODO MOMA from here

docs/src/examples/05-enzyme-constrained-models.jl

@@ -108,7 +108,7 @@ ec_solution = enzyme_constrained_flux_balance_analysis(
 ### Building a model incrementally
 
 # create basic flux model
-m = fbc_model_constraints(model)
+m = build_flux_balance_model(model)
 
 # create enzyme variables
 m += :enzymes^enzyme_variables(model)

src/COBREXA.jl

@@ -36,11 +36,14 @@ include("builders/objectives.jl")
 include("builders/enzymes.jl")
 include("builders/thermodynamic.jl")
 
-include("analysis/modifications.jl")
+include("analysis/core.jl")
 include("analysis/flux_balance.jl")
 include("analysis/parsimonious_flux_balance.jl")
 include("analysis/minimal_metabolic_adjustment.jl")
+include("analysis/enzyme_constraints.jl")
+include("analysis/max_min_driving_force.jl")
 
 include("misc/bounds.jl")
+include("misc/modifications.jl")
 
 end # module COBREXA

src/analysis/core.jl

@@ -0,0 +1,27 @@
+
+"""
+$(TYPEDSIGNATURES)
+
+Make an JuMP model out of `constraints` using [`optimization_model`](@ref)
+(most arguments are forwarded there), then apply the `modifications`, optimize
+the model, and return either `nothing` if the optimization failed, or `output`
+substituted with the solved values (`output` defaults to `constraints`.
+
+For a "nice" version for simpler finding of metabolic model optima, use
+[`flux_balance`](@ref).
+"""
+function optimized_constraints(
+    constraints::C.ConstraintTreeElem;
+    modifications = [],
+    output = constraints,
+    kwargs...,
+)
+    om = optimization_model(constraints; kwargs...)
+    for m in modifications
+        m(om)
+    end
+    J.optimize!(om)
+    is_solved(om) ? C.constraint_values(output, J.value.(om[:x])) : nothing
+end
+
+export optimized_constraints

src/analysis/enzyme_constraints.jl

@@ -0,0 +1,50 @@
+
+"""
+$(TYPEDSIGNATURES)
+
+Run a basic enzyme constrained flux balance analysis on `model`. The enzyme
+model is parameterized by `reaction_isozymes`, which is a mapping of reaction
+IDs (those used in the fluxes of the model) to named [`Isozyme`](@ref)s.
+Additionally, `gene_molar_masses` and `capacity_limitations` should be supplied.
+The latter is a vector of tuples, where each tuple represents a distinct bound
+as `(bound_id, genes_in_bound, protein_mass_bound)`. Typically, `model` has
+bounded exchange reactions, which are unnecessary in enzyme constrained models.
+Unbound these reactions by listing their IDs in `unconstrain_reactions`, which
+makes them reversible. Optimization `modifications` are directly forwarded to
+[`optimized_constraints`](@ref).
+
+In the event that your model requires more complex build steps, consider
+constructing it manually by using [`add_enzyme_constraints!`](@ref).
+
+Returns a tree with the optimization solution or `nothing` if the model cannot
+be solved.
+"""
+function enzyme_constrained_flux_balance_analysis(
+    model::A.AbstractFBCModel,
+    reaction_isozymes::Dict{String,Dict{String,Isozyme}},
+    gene_molar_masses::Dict{String,Float64},
+    capacity_limitations::Vector{Tuple{String,Vector{String},Float64}};
+    optimizer,
+    unconstrain_reactions = String[],
+    modifications = [],
+)
+    m = build_flux_balance_model(model)
+
+    # create enzyme variables
+    m += :enzymes^enzyme_variables(model)
+
+    m = add_enzyme_constraints!(
+        m,
+        reaction_isozymes,
+        gene_molar_masses,
+        capacity_limitations,
+    )
+
+    for rid in Symbol.(unconstrain_reactions)
+        m.fluxes[rid].bound = (-1000.0, 1000.0)
+    end
+
+    optimized_constraints(m; objective = m.objective.value, optimizer, modifications)
+end
+
+export enzyme_constrained_flux_balance_analysis

src/analysis/flux_balance.jl

@@ -1,29 +1,3 @@
-"""
-$(TYPEDSIGNATURES)
-
-Make an JuMP model out of `constraints` using [`optimization_model`](@ref)
-(most arguments are forwarded there), then apply the `modifications`, optimize
-the model, and return either `nothing` if the optimization failed, or `output`
-substituted with the solved values (`output` defaults to `constraints`.
-
-For a "nice" version for simpler finding of metabolic model optima, use
-[`flux_balance`](@ref).
-"""
-function optimized_constraints(
-    constraints::C.ConstraintTreeElem;
-    modifications = [],
-    output = constraints,
-    kwargs...,
-)
-    om = optimization_model(constraints; kwargs...)
-    for m in modifications
-        m(om)
-    end
-    J.optimize!(om)
-    is_solved(om) ? C.constraint_values(output, J.value.(om[:x])) : nothing
-end
-
-export optimized_constraints
 
 """
 $(TYPEDSIGNATURES)
@@ -35,8 +9,8 @@ Most arguments are forwarded to [`optimized_constraints`](@ref).
 Returns a tree with the optimization solution of the same shape as
 given by [`fbc_model_constraints`](@ref).
 """
-function flux_balance(model::A.AbstractFBCModel, optimizer; kwargs...)
-    constraints = fbc_model_constraints(model)
+function flux_balance_analysis(model::A.AbstractFBCModel, optimizer; kwargs...)
+    constraints = build_flux_balance_model(model)
     optimized_constraints(
         constraints;
         objective = constraints.objective.value,
@@ -50,6 +24,7 @@ $(TYPEDSIGNATURES)
 
 Pipe-able overload of [`flux_balance`](@ref).
 """
-flux_balance(optimizer; modifications = []) = m -> flux_balance(m, optimizer; modifications)
+flux_balance_analysis(optimizer; modifications = []) =
+    m -> flux_balance_analysis(m, optimizer; modifications)
 
-export flux_balance
+export flux_balance_analysis

src/analysis/max_min_driving_force.jl

@@ -0,0 +1,107 @@
+
+"""
+$(TYPEDSIGNATURES)
+
+Perform a max-min driving force analysis using `optimizer` on the `model` with
+supplied ΔG⁰s in `reaction_standard_gibbs_free_energies`, as defined by Noor, et
+al., "Pathway thermodynamics highlights kinetic obstacles in central
+metabolism.", PLoS computational biology, 2014.
+
+Optionally, `reference_flux` can be used to set the directions of each reaction
+in `model` (all reactions are assumed to proceed forward by default). The
+supplied `reference_flux` should be free of internal cycles i.e.
+thermodynamically consistent. This optional input is important if a reaction in
+`model` normally runs in reverse (negative flux). Note, only the signs are
+extracted from this input, so be careful with floating point precision near 0.
+
+The max-min driving force algorithm returns the Gibbs free energy of the
+reactions, the concentrations of metabolites and the actual maximum minimum
+driving force. The optimization problem solved is:
+```
+max min -ΔᵣG
+s.t. ΔᵣG = ΔrG⁰ + R T S' ln(C)
+     ΔᵣG ≤ 0
+     ln(Cₗ) ≤ ln(C) ≤ ln(Cᵤ)
+```
+where `ΔrG` are the Gibbs energies dissipated by the reactions, R is the gas
+constant, T is the temperature, S is the stoichiometry of the model, and C is
+the vector of metabolite concentrations (and their respective lower and upper
+bounds).
+
+In case no feasible solution exists, `nothing` is returned.
+
+Reactions specified in `ignore_reaction_ids` are internally ignored when
+calculating the max-min driving force. Importantly, this should include water
+and proton importers.
+
+Since biochemical thermodynamics are assumed, the `proton_ids` and `water_ids`
+need to be specified so that they can be ignored in the calculations.
+Effectively this assumes an aqueous environment at constant pH is used.
+
+`constant_concentrations` is used to fix the concentrations of certain
+metabolites (such as CO₂) by passing a dictionary mapping metabolite id to its
+constant value. `concentration_ratios` is used to specify additional constraints
+on metabolite pair concentrations (typically, this is done with various
+cofactors, such as the ATP/ADP ratio. For example, you can fix the concentration
+of ATP to be always 5× higher than of ADP by specifying `Dict("atp_ratio" =>
+("atp_c","adp_c", 5.0))` if these metabolites are called `atp_c` and `adp_c` in
+the model. `concentration_lb` and `concentration_ub` set the `Cₗ` and `Cᵤ` in
+the optimization problems (these are overwritten by `constant_concentrations` if
+supplied).
+
+`T` and `R` can be specified in the corresponding units; defaults are K and
+kJ/K/mol. The unit of metabolite concentrations is typically molar. As usual,
+optimizer settings can be changed with `modifications`.
+
+Returns a tree with the optimization solution or `nothing` if the model cannot
+be solved.
+"""
+function max_min_driving_force_analysis(
+    model::A.AbstractFBCModel,
+    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}(),
+    proton_ids,
+    water_ids,
+    concentration_lb = 1e-9, # M
+    concentration_ub = 1e-1, # M
+    T = 298.15, # Kelvin
+    R = 8.31446261815324e-3, # kJ/K/mol
+    ignore_reaction_ids = String[],
+    modifications = [],
+    optimizer,
+)
+    m = build_max_min_driving_force_model(
+        model,
+        reaction_standard_gibbs_free_energies;
+        reference_flux,
+        concentration_lb,
+        concentration_ub,
+        R,
+        T,
+        ignore_reaction_ids,
+        water_ids,
+        proton_ids,
+    )
+
+    for (mid, val) in constant_concentrations
+        m.log_metabolite_concentrations[Symbol(mid)].bound = log(val)
+    end
+
+    m *=
+        :metabolite_ratio_constraints^log_ratio_constraints(
+            concentration_ratios,
+            m.log_metabolite_concentrations,
+        )
+
+
+    optimized_constraints(
+        m;
+        objective = m.max_min_driving_force.value,
+        optimizer,
+        modifications,
+    )
+end
+
+export max_min_driving_force_analysis