From 48cbc6ed0540d8b75504882c7e0d8a9242c5d425 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Sat, 9 Dec 2023 19:55:51 +0100 Subject: [PATCH] split pFBA into CT and easy-modeling part --- docs/src/examples/03-qp-problems.jl | 2 +- src/analysis/flux_balance.jl | 62 ++++---- src/analysis/parsimonious_flux_balance.jl | 171 ++++++++++++---------- src/builders/objectives.jl | 8 +- src/solver.jl | 25 ---- 5 files changed, 133 insertions(+), 135 deletions(-) diff --git a/docs/src/examples/03-qp-problems.jl b/docs/src/examples/03-qp-problems.jl index b0f1c5147..19bdc56d5 100644 --- a/docs/src/examples/03-qp-problems.jl +++ b/docs/src/examples/03-qp-problems.jl @@ -22,7 +22,7 @@ model = A.load(J.JSONFBCModel, "e_coli_core.json") # load the model # Use the convenience function to run standard pFBA -vt = X.parsimonious_flux_balance(model, Clarabel.Optimizer) +vt = X.parsimonious_flux_balance(model, Clarabel.Optimizer; modifications = [X.silence]) # Or use the piping functionality diff --git a/src/analysis/flux_balance.jl b/src/analysis/flux_balance.jl index fd3583ebe..745c99896 100644 --- a/src/analysis/flux_balance.jl +++ b/src/analysis/flux_balance.jl @@ -1,34 +1,40 @@ """ $(TYPEDSIGNATURES) -Run flux balance analysis (FBA) on the `model`, optionally specifying -`modifications` to the problem. Basically, FBA solves this optimization -problem: -``` -max cᵀx -s.t. S x = b - xₗ ≤ x ≤ xᵤ -``` -See "Orth, J., Thiele, I. & Palsson, B. What is flux balance analysis?. Nat -Biotechnol 28, 245-248 (2010). https://doi.org/10.1038/nbt.1614" for more -information. - -The `optimizer` must be set to a `JuMP`-compatible optimizer, such as -`GLPK.Optimizer` or `Tulip.Optimizer`. - -Optionally, you may specify one or more modifications to be applied to the model -before the analysis, such as [`set_objective_sense`](@ref), -[`set_optimizer`](@ref), [`set_optimizer_attribute`](@ref), and -[`silence`](@ref). - -Returns a tree with the optimization solution of the same shape as the model -defined by [`fbc_model_constraints`](@ref). - -# Example -``` -model = load_model("e_coli_core.json") -solution = flux_balance(model, GLPK.optimizer) -``` +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, + args...; + modifications = [], + output = constraints, + kwargs..., +) + om = optimization_model(constraints, args...; 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) + +Compute an optimal objective-optimizing solution of the given `model`. + +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) diff --git a/src/analysis/parsimonious_flux_balance.jl b/src/analysis/parsimonious_flux_balance.jl index 26d842a6e..0658b0d73 100644 --- a/src/analysis/parsimonious_flux_balance.jl +++ b/src/analysis/parsimonious_flux_balance.jl @@ -1,92 +1,111 @@ + """ $(TYPEDSIGNATURES) -Run parsimonious flux balance analysis (pFBA) on the `model`. In short, pFBA -runs two consecutive optimization problems. The first is traditional FBA: -``` -max cᵀx = μ -s.t. S x = b - xₗ ≤ x ≤ xᵤ -``` -And the second is a quadratic optimization problem: -``` -min Σᵢ xᵢ² -s.t. S x = b - xₗ ≤ x ≤ xᵤ - μ = μ⁰ -``` -Where the optimal solution of the FBA problem, μ⁰, has been added as an -additional constraint. See "Lewis, Nathan E, Hixson, Kim K, Conrad, Tom M, -Lerman, Joshua A, Charusanti, Pep, Polpitiya, Ashoka D, Adkins, Joshua N, -Schramm, Gunnar, Purvine, Samuel O, Lopez-Ferrer, Daniel, Weitz, Karl K, Eils, -Roland, König, Rainer, Smith, Richard D, Palsson, Bernhard Ø, (2010) Omic data -from evolved E. coli are consistent with computed optimal growth from -genome-scale models. Molecular Systems Biology, 6. 390. doi: -accession:10.1038/msb.2010.47" for more details. - -pFBA gets the model optimum by standard FBA (using -[`flux_balance`](@ref) with `optimizer` and `modifications`), then -finds a minimal total flux through the model that still satisfies the (slightly -relaxed) optimum. This is done using a quadratic problem optimizer. If the -original optimizer does not support quadratic optimization, it can be changed -using the callback in `qp_modifications`, which are applied after the FBA. See -the documentation of [`flux_balance`](@ref) for usage examples of -modifications. - -The optimum relaxation sequence can be specified in `relax` parameter, it -defaults to multiplicative range of `[1.0, 0.999999, ..., 0.99]` of the original -bound. - -Returns an optimized model that contains the pFBA solution (or an unsolved model -if something went wrong). - -# Performance - -This implementation attempts to save time by executing all pFBA steps on a -single instance of the optimization model problem, trading off possible -flexibility. For slightly less performant but much more flexible use, one can -construct parsimonious models directly using -[`with_parsimonious_objective`](@ref). - -# Example -``` -model = load_model("e_coli_core.json") -parsimonious_flux_balance(model, biomass, Gurobi.Optimizer) |> values_vec -``` +Optimize the system of `constraints` to get the optimal `objective` value. Then +try to find a "parsimonious" solution with the same `objective` value, which +optimizes the `parsimonious_objective` (possibly also switching optimization +sense, optimizer, and adding more modifications). + +For efficiency, everything is performed on a single instance of JuMP model. + +A simpler version suitable for direct work with metabolic models is available +in [`parsimonious_flux_balance`](@ref). """ -function parsimonious_flux_balance( - model::C.ConstraintTree, - optimizer; +function parsimonious_optimized_constraints( + constraints::C.ConstraintTreeElem, + args...; + objective::C.Value, modifications = [], - qp_modifications = [], - relax_bounds = [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99], + parsimonious_objective::C.Value, + parsimonious_optimizer = nothing, + parsimonious_sense = J.MIN_SENSE, + parsimonious_modifications = [], + tolerances = [absolute_tolerance_bound(0)], + output = constraints, + kwargs..., ) - # Run FBA - opt_model = flux_balance(model, optimizer; modifications) - J.is_solved(opt_model) || return nothing # FBA failed - # get the objective - Z = J.objective_value(opt_model) - original_objective = J.objective_function(opt_model) + # first solve the optimization problem with the original objective + om = optimization_model(constraints, args...; kwargs...) + for m in modifications + m(om) + end + J.optimize!(om) + is_solved(om) || return nothing + + target_objective_value = J.objective_value(om) - # prepare the model for pFBA - for mod in qp_modifications - mod(model, opt_model) + # switch to parsimonizing the solution w.r.t. to the objective value + isnothing(parsimonious_optimizer) || J.set_optimizer(om, parsimonious_optimizer) + for m in parsimonious_modifications + m(om) end - # add the minimization constraint for total flux - v = opt_model[:x] # fluxes - J.@objective(opt_model, Min, sum(dot(v, v))) + J.@objective(om, J.MIN_SENSE, C.substitute(parsimonious_objective, om[:x])) - for rb in relax_bounds - # lb, ub = objective_bounds(rb)(Z) - J.@constraint(opt_model, pfba_constraint, lb <= original_objective <= ub) + # try all admissible tolerances + for tolerance in tolerances + (lb, ub) = tolerance(target_objective_value) + J.@constraint( + om, + pfba_tolerance_constraint, + lb <= C.substitute(objective, om[:x]) <= ub + ) - J.optimize!(opt_model) - J.is_solved(opt_model) && break + J.optimize!(om) + is_solved(om) && return C.constraint_values(output, J.value.(om[:x])) - J.delete(opt_model, pfba_constraint) - J.unregister(opt_model, :pfba_constraint) + J.delete(om, pfba_tolerance_constraint) + J.unregister(om, :pfba_tolerance_constraint) end + # all tolerances failed + return nothing end + +export parsimonious_optimized_constraints + +""" +$(TYPEDSIGNATURES) + +Compute a parsimonious flux solution for the given `model`. In short, the +objective value of the parsimonious solution should be the same as the one from +[`flux_balance`](@ref), except the squared sum of reaction fluxes is minimized. +If there are multiple possible fluxes that achieve a given objective value, +parsimonious flux thus represents the "minimum energy" one, thus arguably more +realistic. + +Most arguments are forwarded to [`parsimonious_optimized_constraints`](@ref), +with some (objectives) filled in automatically to fit the common processing of +FBC models, and some (`tolerances`) provided with more practical defaults. + +Similarly to the [`flux_balance`](@ref), returns a tree with the optimization +solutions of the shape as given by [`fbc_model_constraints`](@ref). +""" +function parsimonious_flux_balance( + model::A.AbstractFBCModel, + optimizer; + tolerances = relative_tolerance_bound.(1 .- [0, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2]), + kwargs..., +) + constraints = fbc_model_constraints(model) + parsimonious_optimized_constraints( + constraints; + optimizer, + objective = constraints.objective.value, + parsimonious_objective = squared_sum_objective(constraints.fluxes), + tolerances, + kwargs..., + ) +end + +""" +$(TYPEDSIGNATURES) + +Pipe-able variant of [`parsimonious_flux_balance`](@ref). +""" +parsimonious_flux_balance(optimizer; kwargs...) = + model -> parsimonious_flux_balance(model, optimizer; kwargs...) + +export parsimonious_flux_balance diff --git a/src/builders/objectives.jl b/src/builders/objectives.jl index c27ca5538..e6f7a8822 100644 --- a/src/builders/objectives.jl +++ b/src/builders/objectives.jl @@ -13,9 +13,7 @@ $(TYPEDSIGNATURES) TODO """ squared_sum_error_objective(constraints::C.ConstraintTree, target::Dict{Symbol,Float64}) = - C.Constraint( - sum( - (C.value(c) - target[k]) * (C.value(c) - target[k]) for - (k, c) in constraints if haskey(target, k) - ), + sum( + (C.squared(C.value(c) - target[k]) for (k, c) in constraints if haskey(target, k)), + init = zero(C.LinearValue), ) diff --git a/src/solver.jl b/src/solver.jl index 4a6e41aef..b0d24fd52 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -47,28 +47,3 @@ is_solved(opt_model::J.Model) = J.termination_status(opt_model) in [J.MOI.OPTIMAL, J.MOI.LOCALLY_SOLVED] export is_solved - -""" -$(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`. -""" -function optimized_constraints( - constraints::C.ConstraintTreeElem, - args...; - modifications = [], - output = constraints, - kwargs..., -) - om = optimization_model(constraints, args...; 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