LCSB-BioCore · stelmo · Dec 28, 2023 · Dec 28, 2023 · Dec 28, 2023
diff --git a/docs/src/examples/02-flux-balance-analysis.jl b/docs/src/examples/02-flux-balance-analysis.jl
@@ -1,7 +1,7 @@
 
 # # Flux balance analysis (FBA)
 
-# Here we use [`flux_balance`](@ref) and several related functions to
+# Here we use [`flux_balance_analysis`](@ref) and several related functions to
 # find an optimal flux in the *E. coli* "core" model. We will need the model,
 # which we can download using [`download_model`](@ref):
 
@@ -26,9 +26,9 @@ model = load_model("e_coli_core.json")
 # There are many possibilities on how to arrange the metabolic model into the
 # optimization framework and how to actually solve it. The "usual" assumed one
 # is captured in the default behavior of function
-# [`flux_balance`](@ref):
+# [`flux_balance_analysis`](@ref):
 
-solution = flux_balance(model, Tulip.Optimizer)
+solution = flux_balance_analysis(model, Tulip.Optimizer)
 
 @test isapprox(solution.objective, 0.8739, atol = TEST_TOLERANCE) #src
 

diff --git a/docs/src/examples/02a-optimizer-parameters.jl b/docs/src/examples/02a-optimizer-parameters.jl
@@ -3,7 +3,7 @@
 #
 # Many optimizers require fine-tuning to produce best results. You can pass in
 # additional optimizer settings via the `modifications` parameter of
-# [`flux_balance`](@ref). These include e.g.
+# [`flux_balance_analysis`](@ref). These include e.g.
 #
 # - [`set_optimizer_attribute`](@ref) (typically allowing you to tune e.g.
 #   iteration limits, tolerances, or floating-point precision)
@@ -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)],

diff --git a/docs/src/examples/02b-model-modifications.jl b/docs/src/examples/02b-model-modifications.jl
@@ -59,11 +59,11 @@ model.reactions["CS"].stoichiometry
 
 # ## Running FBA on modified models
 #
-# Since the canonical model is completely mutable, you can change it in any way you like and feed it directly into [`flux_balance`](@ref). Let's first find a "original" solution, so that we have a base solution for comparing:
+# Since the canonical model is completely mutable, you can change it in any way you like and feed it directly into [`flux_balance_analysis`](@ref). Let's first find a "original" solution, so that we have a base solution for comparing:
 
 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

diff --git a/docs/src/examples/03-parsimonious-flux-balance.jl b/docs/src/examples/03-parsimonious-flux-balance.jl
@@ -1,7 +1,7 @@
 
 # # Parsimonious flux balance analysis
 
-# We will use [`parsimonious_flux_balance`](@ref) and
+# We will use [`parsimonious_flux_balance_analysis`](@ref) and
 # [`minimize_metabolic_adjustment`](@ref) to find the optimal flux
 # distribution in the *E. coli* "core" model.
 #
@@ -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

diff --git a/docs/src/examples/04-minimization-of-metabolic-adjustment.jl b/docs/src/examples/04-minimization-of-metabolic-adjustment.jl
@@ -1,5 +1,5 @@
 
-# # Minimization of metabolic adjustment
+# # Minimization of metabolic adjustment analysis
 
 # TODO MOMA citation
 
@@ -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
diff --git a/src/COBREXA.jl b/src/COBREXA.jl
@@ -31,6 +31,7 @@ include("types.jl")
 include("io.jl")
 include("solver.jl")
 
+# these functions build or extend constrainttrees of metabolic models
 include("builders/core.jl")
 include("builders/genes.jl")
 include("builders/objectives.jl")
@@ -39,11 +40,15 @@ include("builders/thermodynamic.jl")
 include("builders/loopless.jl")
 include("builders/communities.jl")
 
-include("analysis/modifications.jl")
-include("analysis/flux_balance.jl")
-include("analysis/parsimonious_flux_balance.jl")
-include("analysis/minimal_metabolic_adjustment.jl")
+# these are the one shot analysis functions
+include("frontend/flux_balance_analysis.jl")
+include("frontend/parsimonious_flux_balance.jl")
+include("frontend/minimization_of_metabolic_adjustment_analysis.jl")
+include("frontend/enzyme_constrained_flux_balance_analysis.jl")
+include("frontend/loopless_flux_balance_analysis.jl")
+include("frontend/max_min_driving_force_analysis.jl")
 
+include("misc/modifications.jl")
 include("misc/bounds.jl")
 include("misc/utils.jl")
 

diff --git a/src/analysis/flux_balance.jl b/src/analysis/flux_balance.jl
diff --git a/src/analysis/loopless_flux_balance.jl b/src/analysis/loopless_flux_balance.jl
diff --git a/src/builders/core.jl b/src/builders/core.jl
@@ -1,6 +1,4 @@
 
-import SparseArrays: sparse
-
 """
 $(TYPEDSIGNATURES)
 

diff --git a/src/builders/enzymes.jl b/src/builders/enzymes.jl
@@ -234,49 +234,3 @@ function add_enzyme_constraints!(
 end
 
 export add_enzyme_constraints!
-
-"""
-$(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.
-
-In the event that your model requires more complex build steps, consider
-constructing it manually by using [`add_enzyme_constraints!`](@ref).
-"""
-function enzyme_constrained_flux_balance_analysis(
-    model::A.AbstractFBCModel,
-    reaction_isozymes::Dict{String,Dict{String,SimpleIsozyme}},
-    gene_molar_masses::Dict{String,Float64},
-    capacity_limitations::Vector{Tuple{String,Vector{String},Float64}};
-    optimizer,
-    unconstrain_reactions = String[],
-    modifications = [],
-)
-    m = fbc_model_constraints(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 = C.Between(-1000.0, 1000.0)
-    end
-
-    optimized_constraints(m; objective = m.objective.value, optimizer, modifications)
-end
-
-export enzyme_constrained_flux_balance_analysis
diff --git a/src/builders/loopless.jl b/src/builders/loopless.jl
@@ -103,61 +103,3 @@ function add_loopless_constraints!(
 end
 
 export add_loopless_constraints!
-
-"""
-$(TYPEDSIGNATURES)
-
-Add quasi-thermodynamic constraints to the model to ensure that no thermodynamically
-infeasible internal cycles can occur. Adds the following constraints to the problem:
-```
--max_flux_bound × (1 - yᵢ) ≤ xᵢ ≤ max_flux_bound × yᵢ
--max_flux_bound × yᵢ + strict_inequality_tolerance × (1 - yᵢ) ≤ Gᵢ
-Gᵢ ≤ -strict_inequality_tolerance × yᵢ + max_flux_bound × (1 - yᵢ)
-Nᵢₙₜ' × G = 0
-yᵢ ∈ {0, 1}
-Gᵢ ∈ ℝ
-i ∈ internal reactions
-Nᵢₙₜ is the nullspace of the internal stoichiometric matrix
-```
-Note, this modification introduces binary variables, so an optimization solver capable of
-handing mixed integer problems needs to be used. The arguments `max_flux_bound` and
-`strict_inequality_tolerance` implement the "big-M" method of indicator constraints.
-
-For more details about the algorithm, see `Schellenberger, Lewis, and, Palsson. "Elimination
-of thermodynamically infeasible loops in steady-state metabolic models.", Biophysical
-journal, 2011`.
-"""
-function loopless_flux_balance_analysis(
-    model;
-    max_flux_bound = 1000.0, # needs to be an order of magnitude bigger, big M method heuristic
-    strict_inequality_tolerance = 1.0, # heuristic from paper
-    modifications = [],
-    optimizer,
-)
-
-    m = fbc_model_constraints(model)
-
-    # find all internal reactions
-    internal_reactions = [
-        (i, Symbol(rid)) for
-        (i, rid) in enumerate(A.reactions(model)) if !is_boundary(model, rid)
-    ]
-    internal_reaction_ids = last.(internal_reactions)
-    internal_reaction_idxs = first.(internal_reactions) # order needs to match the internal reaction ids below
-
-    internal_reaction_stoichiometry_nullspace_columns =
-        eachcol(nullspace(Array(A.stoichiometry(model)[:, internal_reaction_idxs]))) # no sparse nullspace function
-
-    m = add_loopless_constraints!(
-        m,
-        internal_reaction_ids,
-        internal_reaction_stoichiometry_nullspace_columns;
-        max_flux_bound,
-        strict_inequality_tolerance,
-    )
-
-    # solve
-    optimized_constraints(m; objective = m.objective.value, optimizer, modifications)
-end
-
-export loopless_flux_balance_analysis