From eecce631249b4ae470d32411f4ad6dc97295cd41 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Tue, 30 Jan 2024 13:58:34 +0100 Subject: [PATCH 1/9] add frontend-managing functions --- docs/src/reference/analysis.md | 7 +++ src/COBREXA.jl | 1 + src/analysis/frontend.jl | 91 ++++++++++++++++++++++++++++++++++ 3 files changed, 99 insertions(+) create mode 100644 src/analysis/frontend.jl diff --git a/docs/src/reference/analysis.md b/docs/src/reference/analysis.md index edbe01ce..41a5e084 100644 --- a/docs/src/reference/analysis.md +++ b/docs/src/reference/analysis.md @@ -26,3 +26,10 @@ Pages = ["src/analysis/screen.jl", "src/analysis/variability.jl", "src/analysis/ Modules = [COBREXA] Pages = ["src/analysis/sample.jl"] ``` + +## Analysis front-end API helpers + +```@autodocs +Modules = [COBREXA] +Pages = ["src/analysis/frontend.jl"] +``` diff --git a/src/COBREXA.jl b/src/COBREXA.jl index c3b3feb2..17803826 100644 --- a/src/COBREXA.jl +++ b/src/COBREXA.jl @@ -54,6 +54,7 @@ include("worker_data.jl") # generic analysis functions include("analysis/envelope.jl") +include("analysis/frontend.jl") include("analysis/parsimonious.jl") include("analysis/sample.jl") include("analysis/screen.jl") diff --git a/src/analysis/frontend.jl b/src/analysis/frontend.jl new file mode 100644 index 00000000..6e051cdc --- /dev/null +++ b/src/analysis/frontend.jl @@ -0,0 +1,91 @@ + +# Copyright (c) 2021-2024, University of Luxembourg +# Copyright (c) 2021-2024, Heinrich-Heine University Duesseldorf +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" +$(TYPEDSIGNATURES) + +A helper that converts a front-end constraint `builder` function (the output of +which would normally be just passed through [`optimized_values`](@ref)) to +front-end analysis function. +""" +function frontend_optimized_values( + builder, + args...; + objective = identity, + output = identity, + sense, + optimizer, + settings = [], + kwargs..., +) + constraints = builder(args...; kwargs...) + + # arguments need to be kept in sync + optimized_values( + constraints; + objective = objective(constraints), + output = output(constraints), + sense, + settings, + optimizer, + ) +end + +export frontend_optimized_values + +""" +$(TYPEDSIGNATURES) + +A helper that converts a parsimonious-style front-end constraint `builder` +function to front-end analysis function. + +Like [`frontend_optimized_values`](@ref), but internally calls +[`parsimonious_optimized_values`](@ref). +""" +function frontend_parsimonious_optimized_values( + builder, + args...; + objective = identity, + output = identity, + sense = Maximal, + optimizer, + settings = [], + parsimonious_objective, + parsimonious_optimizer = nothing, + parsimonious_sense = Minimal, + parsimonious_settings = [], + tolerances = [absolute_tolerance_bound(0)], + kwargs..., +) + constraints = builder(args...; kwargs...) + + # arguments need to be kept in sync + parsimonious_optimized_values( + constraints; + objective = objective(constraints), + output = output(constraints), + sense, + settings, + optimizer, + parsimonious_objective = parsimonious_objective(constraints), + parsimonious_optimizer, + parsimonious_sense, + parsimonious_settings, + tolerances, + ) +end + +export frontend_parsimonious_optimized_values From 272019c9426b681ca09098b6cdc7d04e2ede0189 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Tue, 30 Jan 2024 13:58:54 +0100 Subject: [PATCH 2/9] add a note about keeping args in sync --- src/analysis/parsimonious.jl | 4 +++- src/analysis/solver.jl | 2 ++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/analysis/parsimonious.jl b/src/analysis/parsimonious.jl index 9618dde0..2fe99abb 100644 --- a/src/analysis/parsimonious.jl +++ b/src/analysis/parsimonious.jl @@ -33,12 +33,14 @@ function parsimonious_optimized_values( settings = [], parsimonious_objective::C.Value, parsimonious_optimizer = nothing, - parsimonious_sense = J.MIN_SENSE, + parsimonious_sense = Minimal, parsimonious_settings = [], tolerances = [absolute_tolerance_bound(0)], output = constraints, kwargs..., ) + # arguments need to be kept in sync with + # frontend_parsimonious_optimized_values # first solve the optimization problem with the original objective om = optimization_model(constraints; objective, kwargs...) diff --git a/src/analysis/solver.jl b/src/analysis/solver.jl index 6798a872..117b6baa 100644 --- a/src/analysis/solver.jl +++ b/src/analysis/solver.jl @@ -31,6 +31,8 @@ function optimized_values( output::C.ConstraintTreeElem = constraints, kwargs..., ) + # arguments need to be kept in sync with frontend_optimized_values + om = optimization_model(constraints; kwargs...) for m in settings m(om) From fde248ceb4cec6287b9ed51073d521fed58ec4f3 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Tue, 30 Jan 2024 13:59:35 +0100 Subject: [PATCH 3/9] use the frontend function to construct FBA from flux_balance_constraints --- docs/src/examples/02-flux-balance-analysis.jl | 2 +- docs/src/examples/02a-optimizer-parameters.jl | 4 ++-- docs/src/examples/02b-model-modifications.jl | 4 ++-- src/analysis/frontend.jl | 6 +++--- src/frontend/balance.jl | 15 ++++++--------- 5 files changed, 14 insertions(+), 17 deletions(-) diff --git a/docs/src/examples/02-flux-balance-analysis.jl b/docs/src/examples/02-flux-balance-analysis.jl index 7527c869..a30a51fb 100644 --- a/docs/src/examples/02-flux-balance-analysis.jl +++ b/docs/src/examples/02-flux-balance-analysis.jl @@ -43,7 +43,7 @@ model = load_model("e_coli_core.json") # is captured in the default behavior of function # [`flux_balance_analysis`](@ref): -solution = flux_balance_analysis(model, GLPK.Optimizer) +solution = flux_balance_analysis(model, optimizer = GLPK.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 index 0bc355c9..f5d27f56 100644 --- a/docs/src/examples/02a-optimizer-parameters.jl +++ b/docs/src/examples/02a-optimizer-parameters.jl @@ -46,7 +46,7 @@ model = load_model("e_coli_core.json") # limit for IPM algorithm may now look as follows: solution = flux_balance_analysis( model, - Tulip.Optimizer; + optimizer = Tulip.Optimizer, settings = [silence, set_optimizer_attribute("IPM_IterationsLimit", 1000)], ) @@ -59,7 +59,7 @@ solution = flux_balance_analysis( solution = flux_balance_analysis( model, - Tulip.Optimizer; + optimizer = Tulip.Optimizer, settings = [set_optimizer_attribute("IPM_IterationsLimit", 2)], ) diff --git a/docs/src/examples/02b-model-modifications.jl b/docs/src/examples/02b-model-modifications.jl index cb281a9b..0ff693cc 100644 --- a/docs/src/examples/02b-model-modifications.jl +++ b/docs/src/examples/02b-model-modifications.jl @@ -78,7 +78,7 @@ model.reactions["CS"].stoichiometry import GLPK -base_solution = flux_balance_analysis(model, GLPK.Optimizer) +base_solution = flux_balance_analysis(model, optimizer = GLPK.Optimizer) base_solution.objective # Now, for example, we can limit the intake of glucose by the model: @@ -91,7 +91,7 @@ model.reactions["EX_glc__D_e"].lower_bound = -5.0 # ...and solve the modified model: # -low_glucose_solution = flux_balance_analysis(model, GLPK.Optimizer) +low_glucose_solution = flux_balance_analysis(model, optimizer = GLPK.Optimizer) low_glucose_solution.objective @test isapprox(low_glucose_solution.objective, 0.41559777, atol = TEST_TOLERANCE) #src diff --git a/src/analysis/frontend.jl b/src/analysis/frontend.jl index 6e051cdc..ff99e41c 100644 --- a/src/analysis/frontend.jl +++ b/src/analysis/frontend.jl @@ -24,9 +24,9 @@ front-end analysis function. function frontend_optimized_values( builder, args...; - objective = identity, + objective, output = identity, - sense, + sense = Maximal, optimizer, settings = [], kwargs..., @@ -39,8 +39,8 @@ function frontend_optimized_values( objective = objective(constraints), output = output(constraints), sense, - settings, optimizer, + settings, ) end diff --git a/src/frontend/balance.jl b/src/frontend/balance.jl index 0f52bff6..1c43ecee 100644 --- a/src/frontend/balance.jl +++ b/src/frontend/balance.jl @@ -24,14 +24,11 @@ Most arguments are forwarded to [`optimized_values`](@ref). Returns a tree with the optimization solution of the same shape as given by [`flux_balance_constraints`](@ref). """ -function flux_balance_analysis(model::A.AbstractFBCModel, optimizer; kwargs...) - constraints = flux_balance_constraints(model) - optimized_values( - constraints; - objective = constraints.objective.value, - optimizer, - kwargs..., - ) -end +flux_balance_analysis(model::A.AbstractFBCModel; kwargs...) = frontend_optimized_values( + flux_balance_constraints, + model; + objective = x -> x.objective.value, + kwargs..., +) export flux_balance_analysis From 0e78757ee7a2136be6b6c6fcb72e964523f0726c Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 31 Jan 2024 08:36:25 +0100 Subject: [PATCH 4/9] fix parameters of mmdf --- docs/src/examples/06-mmdf.jl | 4 ++-- src/frontend/mmdf.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/examples/06-mmdf.jl b/docs/src/examples/06-mmdf.jl index df3e672c..2111d846 100644 --- a/docs/src/examples/06-mmdf.jl +++ b/docs/src/examples/06-mmdf.jl @@ -107,8 +107,8 @@ reference_flux = Dict( # ## Solving the MMDF problem mmdf_solution = max_min_driving_force_analysis( - model, - reaction_standard_gibbs_free_energies; + model; + reaction_standard_gibbs_free_energies, reference_flux, concentration_ratios = Dict( "atp" => ("atp_c", "adp_c", 10.0), diff --git a/src/frontend/mmdf.jl b/src/frontend/mmdf.jl index 27cf4639..28d723ba 100644 --- a/src/frontend/mmdf.jl +++ b/src/frontend/mmdf.jl @@ -55,8 +55,8 @@ Following arguments are set optionally: forces. """ function max_min_driving_force_analysis( - model::A.AbstractFBCModel, - reaction_standard_gibbs_free_energies::Dict{String,Float64}; + 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}(), From 9ba319d3c88401410de57563d26d1ec3fc2011d9 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 31 Jan 2024 08:36:59 +0100 Subject: [PATCH 5/9] split mmdf constraint system and frontend --- src/frontend/mmdf.jl | 74 ++++++++++++++++++++++++++------------------ 1 file changed, 44 insertions(+), 30 deletions(-) diff --git a/src/frontend/mmdf.jl b/src/frontend/mmdf.jl index 28d723ba..759431f7 100644 --- a/src/frontend/mmdf.jl +++ b/src/frontend/mmdf.jl @@ -17,21 +17,22 @@ """ $(TYPEDSIGNATURES) -Perform a max-min driving force analysis using `optimizer` on the `model` with -supplied reaction standard Gibbs energies in `reaction_standard_gibbs_free_energies`. +Create max-min driving force constraint system from `model`, using the supplied +reaction standard Gibbs free energies in `reaction_standard_gibbs_free_energies`. -The method was described by by: Noor, et al., "Pathway thermodynamics highlights -kinetic obstacles in central metabolism.", PLoS computational biology, 2014. +The method is described by: *Noor, et al., "Pathway thermodynamics highlights +kinetic obstacles in central metabolism.", PLoS computational biology, 2014.* `reference_flux` sets the directions of each reaction in `model`. The scale of the values is not important, only the direction is examined (w.r.t. `reference_flux_atol` tolerance). Ideally, the supplied `reference_flux` should be completely free of internal cycles, which enables the thermodynamic consistency. To get the cycle-free flux, you can use -[`loopless_flux_balance_analysis`](@ref) (computationally complex but gives -very good fluxes), [`parsimonious_flux_balance_analysis`](@ref) or -[`linear_parsimonious_flux_balance_analysis`](@ref) (computationally simplest -but the consistency is not guaranteed). +[`loopless_flux_balance_analysis`](@ref) (computationally demanding, but gives +thermodynamically consistent solutions), +[`parsimonious_flux_balance_analysis`](@ref) or +[`linear_parsimonious_flux_balance_analysis`](@ref) (which is computationally +simple, but the consistency is not guaranteed). Internally, [`log_concentration_constraints`](@ref) is used to lay out the base structure of the problem. @@ -54,7 +55,7 @@ Following arguments are set optionally: log-concentrations to obtain the actual Gibbs energies and thus driving forces. """ -function max_min_driving_force_analysis( +function max_min_driving_force_constraints( model::A.AbstractFBCModel; reaction_standard_gibbs_free_energies::Dict{String,Float64}, reference_flux = Dict{String,Float64}(), @@ -69,8 +70,6 @@ function max_min_driving_force_analysis( R = 8.31446261815324e-3, # kJ/K/mol reference_flux_atol = 1e-6, check_ignored_reactions = missing, - settings = [], - optimizer, ) # First let's just check if all the identifiers are okay because we use @@ -152,6 +151,9 @@ function max_min_driving_force_analysis( Set(Symbol.(ignored_metabolites)), ) + # TODO it might be quite viable to only create a system with the + # metabolites and reactions actually required for the MMDF. + constraints = log_concentration_constraints( model, @@ -180,26 +182,38 @@ function max_min_driving_force_analysis( end for (rid, dGr0) in reaction_standard_gibbs_free_energies ) - constraints = - constraints * - :driving_forces^driving_forces * - :min_driving_force_thresholds^C.map(driving_forces) do c - less_or_equal_constraint(constraints.min_driving_force, c) - end * - :concentration_ratio_constraints^C.ConstraintTree( - Symbol(cid) => difference_constraint( - constraints.log_concentrations[Symbol(m1)], - constraints.log_concentrations[Symbol(m2)], - log(ratio), - ) for (cid, (m1, m2, ratio)) in concentration_ratios - ) - - optimized_values( - constraints; - objective = constraints.min_driving_force.value, - optimizer, - settings, + constraints * + :driving_forces^driving_forces * + :min_driving_force_thresholds^C.map(driving_forces) do c + less_or_equal_constraint(constraints.min_driving_force, c) + end * + :concentration_ratio_constraints^C.ConstraintTree( + Symbol(cid) => difference_constraint( + constraints.log_concentrations[Symbol(m1)], + constraints.log_concentrations[Symbol(m2)], + log(ratio), + ) for (cid, (m1, m2, ratio)) in concentration_ratios ) end +export max_min_driving_force_constraints + +""" +$(TYPEDSIGNATURES) + +Perform the max-min driving force analysis on the `model`, returning the solved +constraint system. + +Arguments are forwarded to [`max_min_driving_force_constraints`](@ref) (see the +documentation for the description of the constraint system); solver +configuration arguments are forwarded to [`optimized_values`](@ref). +""" +max_min_driving_force_analysis(model::A.AbstractFBCModel; kwargs...) = + frontend_optimized_values( + max_min_driving_force_constraints, + model; + objective = x -> x.min_driving_force.value, + kwargs..., + ) + export max_min_driving_force_analysis From c79a3d4a0ab9db5a314d4752f2c2d4c42b71d3af Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 31 Jan 2024 09:01:30 +0100 Subject: [PATCH 6/9] split off the enzyme FBA frontend --- src/frontend/enzymes.jl | 129 +++++++++++++++++++++------------------- 1 file changed, 69 insertions(+), 60 deletions(-) diff --git a/src/frontend/enzymes.jl b/src/frontend/enzymes.jl index 5642a52c..11e2b5f0 100644 --- a/src/frontend/enzymes.jl +++ b/src/frontend/enzymes.jl @@ -35,26 +35,24 @@ export Isozyme """ $(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 +Construct a enzyme-constrained flux-balance constraint system. The model is +parameterized by `reaction_isozymes`, which is a mapping of reaction identifiers to [`Isozyme`](@ref) descriptions. -Additionally, one typically wants to supply `gene_product_molar_masses` to -describe the weights of enzyme building material, and `capacity` which limits -the mass of enzymes in the whole model. +Additionally, the computation requires `gene_product_molar_masses` to describe +the weights of enzyme building material, and `capacity`, which limits the mass +of enzymes in the whole model. `capacity` may be a single number, which sets the limit for "all described enzymes". Alternatively, `capacity` may be a vector of identifier-genes-limit triples that make a constraint (identified by the given identifier) that limits the listed genes to the given limit. """ -function enzyme_constrained_flux_balance_analysis( +function enzyme_constrained_flux_balance_constraints( model::A.AbstractFBCModel; reaction_isozymes::Dict{String,Dict{String,Isozyme}}, gene_product_molar_masses::Dict{String,Float64}, capacity::Union{Vector{Tuple{String,Vector{String},Float64}},Float64}, - optimizer, - settings = [], ) constraints = flux_balance_constraints(model) @@ -78,63 +76,74 @@ function enzyme_constrained_flux_balance_analysis( ) # connect all parts with constraints - constraints = - constraints * - :directional_flux_balance^sign_split_constraints( - positive = constraints.fluxes_forward, - negative = constraints.fluxes_reverse, - signed = constraints.fluxes, - ) * - :isozyme_flux_forward_balance^isozyme_flux_constraints( - constraints.isozyme_forward_amounts, - constraints.fluxes_forward, - (rid, isozyme) -> maybemap( - x -> x.kcat_forward, - maybeget(reaction_isozymes, string(rid), string(isozyme)), - ), - ) * - :isozyme_flux_reverse_balance^isozyme_flux_constraints( - constraints.isozyme_reverse_amounts, - constraints.fluxes_reverse, - (rid, isozyme) -> maybemap( - x -> x.kcat_reverse, - maybeget(reaction_isozymes, string(rid), string(isozyme)), - ), - ) * - :gene_product_isozyme_balance^gene_product_isozyme_constraints( - constraints.gene_product_amounts, - (constraints.isozyme_forward_amounts, constraints.isozyme_reverse_amounts), - (rid, isozyme) -> maybemap( - x -> [(Symbol(k), v) for (k, v) in x.gene_product_stoichiometry], - maybeget(reaction_isozymes, string(rid), string(isozyme)), + constraints * + :directional_flux_balance^sign_split_constraints( + positive = constraints.fluxes_forward, + negative = constraints.fluxes_reverse, + signed = constraints.fluxes, + ) * + :isozyme_flux_forward_balance^isozyme_flux_constraints( + constraints.isozyme_forward_amounts, + constraints.fluxes_forward, + (rid, isozyme) -> maybemap( + x -> x.kcat_forward, + maybeget(reaction_isozymes, string(rid), string(isozyme)), + ), + ) * + :isozyme_flux_reverse_balance^isozyme_flux_constraints( + constraints.isozyme_reverse_amounts, + constraints.fluxes_reverse, + (rid, isozyme) -> maybemap( + x -> x.kcat_reverse, + maybeget(reaction_isozymes, string(rid), string(isozyme)), + ), + ) * + :gene_product_isozyme_balance^gene_product_isozyme_constraints( + constraints.gene_product_amounts, + (constraints.isozyme_forward_amounts, constraints.isozyme_reverse_amounts), + (rid, isozyme) -> maybemap( + x -> [(Symbol(k), v) for (k, v) in x.gene_product_stoichiometry], + maybeget(reaction_isozymes, string(rid), string(isozyme)), + ), + ) * + :gene_product_capacity^( + capacity isa Float64 ? + C.Constraint( + value = sum( + gpa.value * gene_product_molar_masses[String(gp)] for + (gp, gpa) in constraints.gene_product_amounts ), - ) * - :gene_product_capacity^( - capacity isa Float64 ? - C.Constraint( + bound = C.Between(0, capacity), + ) : + C.ConstraintTree( + Symbol(id) => C.Constraint( value = sum( - gpa.value * gene_product_molar_masses[String(gp)] for - (gp, gpa) in constraints.gene_product_amounts + constraints.gene_product_amounts[Symbol(gp)].value * + gene_product_molar_masses[gp] for gp in gps ), - bound = C.Between(0, capacity), - ) : - C.ConstraintTree( - Symbol(id) => C.Constraint( - value = sum( - constraints.gene_product_amounts[Symbol(gp)].value * - gene_product_molar_masses[gp] for gp in gps - ), - bound = C.Between(0, limit), - ) for (id, gps, limit) in capacity_limits - ) + bound = C.Between(0, limit), + ) for (id, gps, limit) in capacity_limits ) - - optimized_values( - constraints; - objective = constraints.objective.value, - optimizer, - settings, ) end +export enzyme_constrained_flux_balance_constraints + +""" +$(TYPEDSIGNATURES) + +Perform the enzyme-constrained flux balance analysis on the `model` and return the solved constraint system. + +Arguments are forwarded to +[`enzyme_constrained_flux_balance_constraints`](@ref); solver configuration +arguments are forwarded to [`optimized_values`](@ref). +""" +enzyme_constrained_flux_balance_analysis(model::A.AbstractFBCModel; kwargs...) = + frontend_optimized_values( + enzyme_constrained_flux_balance_constraints, + model; + objective = x -> x.objective.value, + kwargs..., + ) + export enzyme_constrained_flux_balance_analysis From 8642cdfb818b4b8fb5e1538652b8d839a93d3404 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 31 Jan 2024 09:24:47 +0100 Subject: [PATCH 7/9] split off loopless constraints from frontend --- src/frontend/loopless.jl | 65 ++++++++++++++++++++++++---------------- 1 file changed, 40 insertions(+), 25 deletions(-) diff --git a/src/frontend/loopless.jl b/src/frontend/loopless.jl index 7e0494ac..8315b62e 100644 --- a/src/frontend/loopless.jl +++ b/src/frontend/loopless.jl @@ -17,11 +17,12 @@ """ $(TYPEDSIGNATURES) -Perform a flux balance analysis with added quasi-thermodynamic constraints that -ensure that thermodynamically infeasible internal cycles can not occur. The -method is closer described by: Schellenberger, Lewis, and, Palsson. -"Elimination of thermodynamically infeasible loops in steady-state metabolic -models.", Biophysical journal, 2011`. +Construct a flux-balance constraint system from `model` with added +quasi-thermodynamic constraints that ensure that thermodynamically infeasible +internal cycles can not occur. The method is closer described by: +*Schellenberger, Lewis, and, Palsson. "Elimination of thermodynamically +infeasible loops in steady-state metabolic models.", Biophysical journal, +2011`*. The loopless condition comes with a performance overhead: the computation needs to find the null space of the stoichiometry matrix (essentially inverting it); @@ -29,17 +30,18 @@ and the subsequently created optimization problem contains binary variables for each internal reaction, thus requiring a MILP solver and a potentially exponential solving time. +Internally, the system is constructed by combining +[`flux_balance_constraints`](@ref) and [`loopless_constraints`](@ref). + The arguments `driving_force_max_bound` and `driving_force_nonzero_bound` set the bounds (possibly negated ones) on the virtual "driving forces" (G_i in the paper). """ -function loopless_flux_balance_analysis( +function loopless_flux_balance_constraints( model::A.AbstractFBCModel; flux_infinity_bound = 10000.0, driving_force_nonzero_bound = 1.0, driving_force_infinity_bound = 1000.0, - settings = [], - optimizer, ) constraints = flux_balance_constraints(model) @@ -54,24 +56,37 @@ function loopless_flux_balance_analysis( :loopless_directions^C.variables(keys = internal_reactions, bounds = Switch(0, 1)) + :loopless_driving_forces^C.variables(keys = internal_reactions) - constraints *= - :loopless_constraints^loopless_constraints(; - fluxes = constraints.fluxes, - loopless_direction_indicators = constraints.loopless_directions, - loopless_driving_forces = constraints.loopless_driving_forces, - internal_reactions, - internal_nullspace = LinearAlgebra.nullspace(Matrix(stoi[:, internal_mask])), - flux_infinity_bound, - driving_force_nonzero_bound, - driving_force_infinity_bound, - ) - - optimized_values( - constraints; - objective = constraints.objective.value, - optimizer, - settings, + constraints * + :loopless_constraints^loopless_constraints(; + fluxes = constraints.fluxes, + loopless_direction_indicators = constraints.loopless_directions, + loopless_driving_forces = constraints.loopless_driving_forces, + internal_reactions, + internal_nullspace = LinearAlgebra.nullspace(Matrix(stoi[:, internal_mask])), + flux_infinity_bound, + driving_force_nonzero_bound, + driving_force_infinity_bound, ) end +export loopless_flux_balance_constraints + +""" +$(TYPEDSIGNATURES) + +Perform the loopless flux balance analysis on the `model`, returning the solved +constraint system. + +Arguments are forwarded to [`loopless_flux_balance_constraints`](@ref) (see the +documentation for the description of the constraint system); solver +configuration arguments are forwarded to [`optimized_values`](@ref). +""" +loopless_flux_balance_analysis(model::A.AbstractFBCModel; kwargs...) = + frontend_optimized_values( + loopless_flux_balance_constraints, + model; + objective = x -> x.objective.value, + kwargs..., + ) + export loopless_flux_balance_analysis From c253acec2c9b5d9e69cc2e1e5782838fcaca9248 Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 31 Jan 2024 15:48:39 +0100 Subject: [PATCH 8/9] frontend for parsimonious analysis functions --- .../examples/03-parsimonious-flux-balance.jl | 18 +-- ...04-minimization-of-metabolic-adjustment.jl | 12 +- src/frontend/parsimonious.jl | 109 ++++++++++++------ 3 files changed, 89 insertions(+), 50 deletions(-) diff --git a/docs/src/examples/03-parsimonious-flux-balance.jl b/docs/src/examples/03-parsimonious-flux-balance.jl index 9465d5e3..40171b54 100644 --- a/docs/src/examples/03-parsimonious-flux-balance.jl +++ b/docs/src/examples/03-parsimonious-flux-balance.jl @@ -22,16 +22,16 @@ # # TODO pFBA citation -# If it is not already present, download the model and load the package: -import Downloads: download +using COBREXA -!isfile("e_coli_core.json") && - download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json") +download_model( + "http://bigg.ucsd.edu/static/models/e_coli_core.json", + "e_coli_core.json", + "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8", +) # next, load the necessary packages -using COBREXA - import JSONFBCModels import Clarabel # can solve QPs @@ -39,7 +39,11 @@ model = load_model("e_coli_core.json") # load the model # Use the convenience function to run standard pFBA on -vt = parsimonious_flux_balance_analysis(model, Clarabel.Optimizer; settings = [silence]) +vt = parsimonious_flux_balance_analysis( + model; + optimizer = Clarabel.Optimizer, + settings = [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 index 48165fdb..bf73a576 100644 --- a/docs/src/examples/04-minimization-of-metabolic-adjustment.jl +++ b/docs/src/examples/04-minimization-of-metabolic-adjustment.jl @@ -18,12 +18,14 @@ # TODO MOMA citation -import Downloads: download +using COBREXA -!isfile("e_coli_core.json") && - download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json") +download_model( + "http://bigg.ucsd.edu/static/models/e_coli_core.json", + "e_coli_core.json", + "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8", +) -using COBREXA import AbstractFBCModels.CanonicalModel as CM import JSONFBCModels import Clarabel @@ -36,7 +38,7 @@ model = convert(CM.Model, load_model("e_coli_core.json")) reference_fluxes = parsimonious_flux_balance_analysis( model, - Clarabel.Optimizer, + optimizer = Clarabel.Optimizer, settings = [silence], ).fluxes diff --git a/src/frontend/parsimonious.jl b/src/frontend/parsimonious.jl index f94ecf46..538c973f 100644 --- a/src/frontend/parsimonious.jl +++ b/src/frontend/parsimonious.jl @@ -17,57 +17,59 @@ """ $(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_analysis`](@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. The optimized squared distance is present in the result as -`parsimonious_objective`. - -Most arguments are forwarded to [`parsimonious_optimized_values`](@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_analysis`](@ref), returns a tree with the optimization -solutions of the shape as given by [`flux_balance_constraints`](@ref). +A constraint system like from [`flux_balance_constraints`](@ref), but with the +parsimonious objective present under key `parsimonious_objective`. Best used +via [`parsimonious_flux_balance_analysis`](@ref). """ -function parsimonious_flux_balance_analysis( - model::A.AbstractFBCModel, - optimizer; +function parsimonious_flux_balance_constraints(model::A.AbstractFBCModel) + constraints = flux_balance_constraints(model) + + constraints * + :parsimonious_objective^C.Constraint(squared_sum_value(constraints.fluxes)) +end + +export parsimonious_flux_balance_constraints + +""" +$(TYPEDSIGNATURES) + +Compute a parsimonious flux solution for the `model`, using the constraints +given by [`parsimonious_flux_balance_constraints`](@ref). + +In short, the objective value of the parsimonious solution should be the same +as the one from [`flux_balance_analysis`](@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, which is arguably more realistic. + +Solver configuration arguments are forwarded to +[`parsimonious_optimized_values`](@ref). +""" +parsimonious_flux_balance_analysis( + model::A.AbstractFBCModel; tolerances = relative_tolerance_bound.(1 .- [0, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2]), kwargs..., +) = frontend_parsimonious_optimized_values( + parsimonious_flux_balance_constraints, + model; + objective = x -> x.objective.value, + parsimonious_objective = x -> x.parsimonious_objective.value, + tolerances, + kwargs..., ) - constraints = flux_balance_constraints(model) - parsimonious_objective = squared_sum_value(constraints.fluxes) - parsimonious_optimized_values( - constraints * :parsimonious_objective^C.Constraint(parsimonious_objective); - optimizer, - objective = constraints.objective.value, - parsimonious_objective, - tolerances, - kwargs..., - ) -end export parsimonious_flux_balance_analysis """ $(TYPEDSIGNATURES) -Like [`parsimonious_flux_balance_analysis`](@ref), but uses a L1 metric for -solving the parsimonious problem. - -In turn, the solution is often faster, does not require a solver capable of -quadratic objectives, and has many beneficial properties of the usual -parsimonious solutions (such as the general lack of unnecessary loops). On the -other hand, like with plain flux balance analysis there is no strong guarantee -of uniqueness of the solution. +Like [`parsimonious_flux_balance_constraints`](@ref), but uses a L1 metric for +solving the parsimonious problem. The `parsimonious_objective` constraint is +thus linear. """ -function linear_parsimonious_flux_balance_analysis( +function linear_parsimonious_flux_balance_constraints( model::A.AbstractFBCModel, optimizer; - tolerances = relative_tolerance_bound.(1 .- [0, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2]), kwargs..., ) constraints = flux_balance_constraints(model) @@ -94,4 +96,35 @@ function linear_parsimonious_flux_balance_analysis( ) end +export linear_parsimonious_flux_balance_constraints + +""" +$(TYPEDSIGNATURES) + +Like [`parsimonious_flux_balance_analysis`](@ref), but uses the L1-metric +parsimonious system given by +[`linear_parsimonious_flux_balance_constraints`](@ref). + +In turn, the solution is often faster, does not require a solver capable of +quadratic objectives, and has many beneficial properties of the usual +parsimonious solutions (such as the general lack of unnecessary loops). On the +other hand, like with plain flux balance analysis there is no strong guarantee +of uniqueness of the solution. + +Solver configuration arguments are forwarded to +[`parsimonious_optimized_values`](@ref). +""" +linear_parsimonious_flux_balance_analysis( + model::A.AbstractFBCModel; + tolerances = relative_tolerance_bound.(1 .- [0, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2]), + kwargs..., +) = frontend_parsimonious_optimized_values( + linear_parsimonious_flux_balance_constraints, + model; + objective = x -> x.objective.value, + parsimonious_objective = x -> x.parsimonious_objective.value, + tolerances, + kwargs..., +) + export linear_parsimonious_flux_balance_analysis From 73f2852432eaa1d0940d66cbf7a7a9bb14f9b11e Mon Sep 17 00:00:00 2001 From: Mirek Kratochvil Date: Wed, 31 Jan 2024 16:37:44 +0100 Subject: [PATCH 9/9] frontend for MOMA (involves quite some argument juggling) --- .../examples/03-parsimonious-flux-balance.jl | 5 +- src/analysis/frontend.jl | 6 +- src/frontend/moma.jl | 206 +++++++++++------- 3 files changed, 131 insertions(+), 86 deletions(-) diff --git a/docs/src/examples/03-parsimonious-flux-balance.jl b/docs/src/examples/03-parsimonious-flux-balance.jl index 40171b54..8d359788 100644 --- a/docs/src/examples/03-parsimonious-flux-balance.jl +++ b/docs/src/examples/03-parsimonious-flux-balance.jl @@ -16,9 +16,8 @@ # # Parsimonious flux balance analysis -# We will use [`parsimonious_flux_balance_analysis`](@ref) and -# [`minimization_of_metabolic_adjustment`](@ref) to find the optimal flux -# distribution in the *E. coli* "core" model. +# We will use [`parsimonious_flux_balance_analysis`](@ref) to find the optimal +# flux distribution in the *E. coli* "core" model. # # TODO pFBA citation diff --git a/src/analysis/frontend.jl b/src/analysis/frontend.jl index ff99e41c..af4f0d84 100644 --- a/src/analysis/frontend.jl +++ b/src/analysis/frontend.jl @@ -24,6 +24,7 @@ front-end analysis function. function frontend_optimized_values( builder, args...; + builder_kwargs = NamedTuple(), objective, output = identity, sense = Maximal, @@ -31,7 +32,7 @@ function frontend_optimized_values( settings = [], kwargs..., ) - constraints = builder(args...; kwargs...) + constraints = builder(args...; builder_kwargs..., kwargs...) # arguments need to be kept in sync optimized_values( @@ -58,6 +59,7 @@ Like [`frontend_optimized_values`](@ref), but internally calls function frontend_parsimonious_optimized_values( builder, args...; + builder_kwargs = NamedTuple(), objective = identity, output = identity, sense = Maximal, @@ -70,7 +72,7 @@ function frontend_parsimonious_optimized_values( tolerances = [absolute_tolerance_bound(0)], kwargs..., ) - constraints = builder(args...; kwargs...) + constraints = builder(args...; builder_kwargs..., kwargs...) # arguments need to be kept in sync parsimonious_optimized_values( diff --git a/src/frontend/moma.jl b/src/frontend/moma.jl index de90094f..ad887e0e 100644 --- a/src/frontend/moma.jl +++ b/src/frontend/moma.jl @@ -17,9 +17,9 @@ """ $(TYPEDSIGNATURES) -Find a feasible solution of the "minimal metabolic adjustment analysis" (MOMA) +Find a solution of the "minimization of metabolic adjustment" (MOMA) analysis for the `model`, which is the "closest" feasible solution to the given -`reference_fluxes`, in the sense of squared-sum error distance. The minimized +`reference_fluxes`, in the sense of squared-sum distance. The minimized squared distance (the objective) is present in the result tree as `minimal_adjustment_objective`. @@ -33,21 +33,14 @@ objective is constructed via [`squared_sum_error_value`](@ref)). Additional parameters are forwarded to [`optimized_values`](@ref). """ -function minimization_of_metabolic_adjustment( +function metabolic_adjustment_minimization_constraints( model::A.AbstractFBCModel, - reference_fluxes::Dict{Symbol,Float64}, - optimizer; - kwargs..., + reference_fluxes::C.Tree, ) constraints = flux_balance_constraints(model) - objective = - squared_sum_error_value(constraints.fluxes, x -> get(reference_fluxes, x, nothing)) - optimized_values( - constraints * :minimal_adjustment_objective^C.Constraint(objective); - optimizer, - objective, - sense = Minimal, - kwargs..., + constraints * + :minimal_adjustment_objective^C.Constraint( + squared_sum_error_value(constraints.fluxes, reference_fluxes), ) end @@ -55,61 +48,85 @@ end $(TYPEDSIGNATURES) A slightly easier-to-use version of -[`minimization_of_metabolic_adjustment`](@ref) that computes the reference flux -as the optimal solution of the `reference_model`. The reference flux is -calculated using `reference_optimizer` and `reference_modifications`, which -default to the `optimizer` and `settings`. - -Leftover arguments are passed to the overload of -[`minimization_of_metabolic_adjustment`](@ref) that accepts the -reference flux dictionary. +[`metabolic_adjustment_minimization_constraints`](@ref) that computes the +reference flux as the parsimonious optimal solution of the `reference_model`. +The reference flux is calculated using `reference_optimizer` and +`reference_modifications`, which default to the `optimizer` and `settings`. + +Other arguments are forwarded to the internal call of +[`parsimonious_optimized_values`](@ref). + +Returns `nothing` if no feasible solution is found. """ -function minimization_of_metabolic_adjustment( +function metabolic_adjustment_minimization_constraints( model::A.AbstractFBCModel, - reference_model::A.AbstractFBCModel, - optimizer; - reference_optimizer = optimizer, - settings = [], - reference_settings = settings, + reference_model::A.AbstractFBCModel; kwargs..., ) - reference_constraints = flux_balance_constraints(reference_model) - reference_fluxes = optimized_values( + reference_constraints = parsimonious_flux_balance_constraints(reference_model) + reference_fluxes = parsimonious_optimized_values( reference_constraints; - optimizer = reference_optimizer, - settings = reference_settings, + objective = reference_constraints.objective.value, + parsimonious_objective = reference_constraints.parsimonious_objective.value, output = reference_constraints.fluxes, - ) - isnothing(reference_fluxes) && return nothing - minimization_of_metabolic_adjustment( - model, - reference_fluxes, - optimizer; - settings, kwargs..., ) + isnothing(reference_fluxes) && return nothing + metabolic_adjustment_minimization_constraints(model, reference_fluxes) end -export minimization_of_metabolic_adjustment +export metabolic_adjustment_minimization_constraints """ $(TYPEDSIGNATURES) -Like [`minimization_of_metabolic_adjustment`](@ref) but optimizes the L1 norm. +TODO +""" +metabolic_adjustment_minimization_analysis( + model::A.AbstractFBCModel, + args...; + optimizer, + settings, + reference_parsimonious_optimizer = optimizer, + reference_parsimonious_settings = settings, + reference_optimizer = optimizer, + reference_settings = settings, + kwargs..., +) = frontend_optimized_values( + metabolic_adjustment_minimization_constraints, + model, + args...; + builder_kwargs = ( + optimizer = reference_optimizer, + settings = reference_settings, + parsimonious_optimizer = reference_parsimonious_optimizer, + parsimonious_settings = reference_parsimonious_settings, + ), + objective = x -> x.minimal_adjustment_objective.value, + sense = Minimal, + optimizer, + settings, + kwargs..., +) + +export metabolic_adjustment_minimization_analysis + +""" +$(TYPEDSIGNATURES) + +Like [`metabolic_adjustment_minimization_constraints`](@ref) but optimizes the L1 norm. This typically produces a sufficiently good result with less resources, depending on the situation. See documentation of [`linear_parsimonious_flux_balance_analysis`](@ref) for some of the considerations. """ -function linear_minimization_of_metabolic_adjustment( +function linear_metabolic_adjustment_minimization_constraints( model::A.AbstractFBCModel, - reference_fluxes::Dict{Symbol,Float64}, - optimizer; - kwargs..., + reference_fluxes::C.Tree, ) constraints = flux_balance_constraints(model) - difference = C.zip(ct.fluxes, C.Tree(reference_fluxes)) do orig, ref + difference = C.zip(ct.fluxes, reference_fluxes) do orig, ref C.Constraint(orig.value - ref) end @@ -120,50 +137,77 @@ function linear_minimization_of_metabolic_adjustment( # `difference` actually doesn't need to go to the CT, but we include it # anyway to calm the curiosity of good neighbors. - constraints *= :reference_diff^difference - constraints *= - :reference_directional_diff_balance^sign_split_constraints( - positive = constraints.reference_positive_diff, - negative = constraints.reference_negative_diff, - signed = difference, - ) - - objective = - sum_value(constraints.reference_positive_diff, constraints.reference_negative_diff) - - optimized_values( - constraints * :linear_minimal_adjustment_objective^C.Constraint(objective); - optimizer, - objective, - sense = Minimal, - kwargs..., + constraints * + :reference_diff^difference * + :reference_directional_diff_balance^sign_split_constraints( + positive = constraints.reference_positive_diff, + negative = constraints.reference_negative_diff, + signed = difference, + ) * + :linear_minimal_adjustment_objective^C.Constraint( + sum_value(constraints.reference_positive_diff, constraints.reference_negative_diff), ) end -function linear_minimization_of_metabolic_adjustment( +""" +$(TYPEDSIGNATURES) + +Like [`metabolic_adjustment_minimization_constraints`](@ref) but optimizes the L1 norm. +This typically produces a sufficiently good result with less resources, +depending on the situation. See documentation of +[`linear_parsimonious_flux_balance_analysis`](@ref) for some of the +considerations. +""" +function linear_metabolic_adjustment_minimization_constraints( model::A.AbstractFBCModel, - reference_model::A.AbstractFBCModel, - optimizer; - reference_optimizer = optimizer, - settings = [], - reference_settings = settings, + reference_model::A.AbstractFBCModel; kwargs..., ) - reference_constraints = flux_balance_constraints(reference_model) - reference_fluxes = optimized_values( + reference_constraints = parsimonious_flux_balance_constraints(reference_model) + reference_fluxes = parsimonious_optimized_values( reference_constraints; - optimizer = reference_optimizer, - settings = reference_settings, + objective = reference_constraints.objective.value, + parsimonious_objective = reference_constraints.parsimonious_objective.value, output = reference_constraints.fluxes, - ) - isnothing(reference_fluxes) && return nothing - linear_minimization_of_metabolic_adjustment( - model, - reference_fluxes, - optimizer; - settings, kwargs..., ) + isnothing(reference_fluxes) && return nothing + + linear_metabolic_adjustment_minimization_constraints(model, reference_fluxes) end -export linear_minimization_of_metabolic_adjustment +export linear_metabolic_adjustment_minimization_constraints + +""" +$(TYPEDSIGNATURES) + +TODO +""" +linear_metabolic_adjustment_minimization_analysis( + model::A.AbstractFBCModel, + args...; + optimizer, + settings, + reference_parsimonious_optimizer = optimizer, + reference_parsimonious_settings = settings, + reference_optimizer = optimizer, + reference_settings = settings, + kwargs..., +) = frontend_optimized_values( + linear_metabolic_adjustment_minimization_constraints, + model, + args...; + builder_kwargs = ( + optimizer = reference_optimizer, + settings = reference_settings, + parsimonious_optimizer = reference_parsimonious_optimizer, + parsimonious_settings = reference_parsimonious_settings, + ), + objective = x -> x.linear_minimal_adjustment_objective.value, + sense = Minimal, + optimizer, + settings, + kwargs..., +) + +export linear_metabolic_adjustment_minimization_analysis