Unverified Commit 9f174fdf authored by Miroslav Kratochvil's avatar Miroslav Kratochvil
Browse files

bulldoze pFBA

parent 26fe2b1c
"""
pfba(model::MetabolicModel, optimizer; modifications, qp_solver, qp_solver_attributes)
function parsimonious_flux_balance_analysis(
model::MetabolicModel,
optimizer;
modifications = [],
qp_modifications = [],
relax_bounds=[1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99],
)
Run parsimonious flux balance analysis (pFBA) on the `model`.
Note, the `optimizer` must be set to perform the analysis, any JuMP solver will work.
Optionally, specify problem `modifications` as in [`flux_balance_analysis`](@ref).
Also, `qp_solver` can be set to be different from `optimizer`, where the latter is then the LP optimizer only.
Also note that `qp_solver_attributes` is meant to set the attributes for the `qp_solver`.
Any solver attributes changed in `modifications` will only affect he LP solver.
This function automatically relaxes the constraint that the FBA solution objective matches the pFBA solution.
This is iteratively relaxed like 1.0, 0.999999, 0.9999, etc. of the bound until 0.99 when the function fails and returns nothing.
Return a solved pFBA model.
pFBA gets the model optimum by standard FBA (using
[`flux_balance_analysis`](@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.
Thhe 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 `nothing` if the
optimization failed.
# Example
```
......@@ -17,49 +29,44 @@ optimizer = Gurobi.Optimizer
atts = Dict("OutputFlag" => 0)
model = load_model(StandardModel, "iJO1366.json")
biomass = findfirst(model.reactions, "BIOMASS_Ec_iJO1366_WT_53p95M")
sol = pfba(model, biomass, optimizer; solver_attributes=atts)
sol = pfba(model, biomass, Gurobi.optimizer)
```
"""
function parsimonious_flux_balance_analysis(
model::MetabolicModel,
optimizer;
modifications = [(model, opt_model) -> nothing],
qp_solver = (model, opt_model) -> nothing,
qp_solver_attributes = [(model, opt_model) -> nothing],
modifications = [],
qp_modifications = [],
relax_bounds = [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99],
)
# Run FBA
opt_model = flux_balance_analysis(model, optimizer; modifications = modifications)
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
return nothing # FBA failed
# Run pFBA
λ = objective_value(opt_model)
old_objective = COBREXA.JuMP.objective_function(opt_model)
qp_solver(model, opt_model) # change the solver if specified, otherwise default argument does nothing
if typeof(qp_solver_attributes) <: AbstractVector # many modifications
for mod in qp_solver_attributes
mod(model, opt_model)
end
else # single modification
qp_solver_attributes(model, opt_model)
# get the objective
Z = objective_value(opt_model)
original_objective = COBREXA.JuMP.objective_function(opt_model)
# prepare the model for pFBA
for mod in qp_modifications
mod(model, opt_model)
end
# add the minimization constraint for total flux
v = opt_model[:x] # fluxes
@objective(opt_model, Min, sum(dot(v, v)))
for lbconval in [1.0, 0.999999, 0.99999, 0.9999, 0.999, 0.99] # sequentially relax bound for stability
λmin = min(lbconval * λ, λ * 1.0 / lbconval)
λmax = max(lbconval * λ, λ * 1.0 / lbconval)
@constraint(
opt_model,
pfbacon,
λmin <= old_objective <= λmax # in case of negative constraints
)
for rb in relax_bounds
lb, ub = objective_bounds(rb)(Z)
@_models_log @info "pFBA step relaxed to [$lb,$ub]"
@constraint(opt_model, pfba_constraint, lb <= original_objective <= ub)
optimize!(opt_model)
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] &&
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
break
COBREXA.JuMP.delete(opt_model, pfbacon)
COBREXA.JuMP.delete(opt_model, pfba_constraint)
end
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
......@@ -69,55 +76,31 @@ function parsimonious_flux_balance_analysis(
end
"""
parsimonious_flux_balance_analysis_dict(model::MetabolicModel, optimizer; modifications, qp_solver, qp_solver_attributes)
parsimonious_flux_balance_analysis_vec(args...; kwargs...)
Perform parsimonious flux balance analysis on `model` using `optimizer`.
Returns a vector of fluxes in the same order as the reactions in `model`.
Calls [`parsimonious_flux_balance_analysis`](@ref) internally.
Arguments are forwarded to [`parsimonious_flux_balance_analysis`](@ref) internally.
"""
function parsimonious_flux_balance_analysis_vec(
model::MetabolicModel,
optimizer;
modifications = [(model, opt_model) -> nothing],
qp_solver = (model, opt_model) -> nothing,
qp_solver_attributes = [(model, opt_model) -> nothing],
)
opt_model = parsimonious_flux_balance_analysis(
model,
optimizer;
modifications = modifications,
qp_solver = qp_solver,
qp_solver_attributes = qp_solver_attributes,
)
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
return nothing
function parsimonious_flux_balance_analysis_vec(args...; kwargs...)
opt_model = parsimonious_flux_balance_analysis(args...; kwargs...)
isnothing(opt_model) && return nothing
return value.(opt_model[:x])
end
"""
parsimonious_flux_balance_analysis_dict(model::MetabolicModel, optimizer; modifications, qp_solver, qp_solver_attributes)
parsimonious_flux_balance_analysis_dict(model::MetabolicModel, args...; kwargs...)
Perform parsimonious flux balance analysis on `model` using `optimizer`.
Returns a dictionary mapping reaction `id`s to fluxes.
Calls [`parsimonious_flux_balance_analysis`](@ref) internally.
Returns a dictionary mapping the reaction IDs to fluxes.
Arguments are forwarded to [`parsimonious_flux_balance_analysis`](@ref) internally.
"""
function parsimonious_flux_balance_analysis_dict(
model::MetabolicModel,
optimizer;
modifications = [(model, opt_model) -> nothing],
qp_solver = (model, opt_model) -> nothing,
qp_solver_attributes = [(model, opt_model) -> nothing],
)
opt_model = parsimonious_flux_balance_analysis(
model,
optimizer;
modifications = modifications,
qp_solver = qp_solver,
qp_solver_attributes = qp_solver_attributes,
)
COBREXA.JuMP.termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] ||
return nothing
function parsimonious_flux_balance_analysis_dict(model::MetabolicModel, args...; kwargs...)
opt_fluxes = parsimonious_flux_balance_analysis_vec(model, args...; kwargs...)
isnothing(opt_fluxes) && return nothing
return Dict(zip(reactions(model), value.(opt_model[:x])))
return Dict(zip(reactions(model), opt_fluxes))
end
......@@ -13,8 +13,10 @@
change_constraint("EX_glc__D_e", -12, -12),
change_solver_attribute("IPM_IterationsLimit", 500),
],
qp_solver = change_solver(OSQP.Optimizer),
qp_solver_attributes = change_solver_attribute("verbose", false),
qp_modifications = [
change_solver(OSQP.Optimizer),
change_solver_attribute("verbose", false),
],
)
v = parsimonious_flux_balance_analysis_vec(
model,
......@@ -23,8 +25,10 @@
change_constraint("EX_glc__D_e", -12, -12),
change_solver_attribute("IPM_IterationsLimit", 500),
],
qp_solver = change_solver(OSQP.Optimizer),
qp_solver_attributes = change_solver_attribute("verbose", false),
qp_modifications = [
change_solver(OSQP.Optimizer),
change_solver_attribute("verbose", false),
],
)
# The used optimizer doesn't really converge to the same answer everytime
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment