Minimization of metabolic adjustment analysis
Minimization of metabolic adjustment analysis (MOMA) finds a flux solution that is closest to some reference solution. This may correspond to realistic adjustment of living organisms to various perturbations, such as gene knockout or environmental stress.
To demonstrate, let's use the E. coli model.
using COBREXA
download_model(
"http://bigg.ucsd.edu/static/models/e_coli_core.json",
"e_coli_core.json",
"7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
)"e_coli_core.json"We shall use both quadratic and linear solvers – the "closest to some reference solution" typically refers to Euclidean ("L2") distance which requires a QP solver, but Manhattan ("L1") distance is also demonstrated below.
import Clarabel, HiGHSBecause we will have to perform some perturbation, we import the model in canonical Julia structures:
import JSONFBCModels
import AbstractFBCModels.CanonicalModel as CM
ecoli = load_model("e_coli_core.json", CM.Model)AbstractFBCModels.CanonicalModel.Model(
reactions = Dict{String, AbstractFBCModels.CanonicalModel.Reaction}("ACALD" =…
metabolites = Dict{String, AbstractFBCModels.CanonicalModel.Metabolite}("glu_…
genes = Dict{String, AbstractFBCModels.CanonicalModel.Gene}("b4301" => Abstra…
couplings = Dict{String, AbstractFBCModels.CanonicalModel.Coupling}(),
)
This will be a good reaction for perturbing:
ecoli.reactions["CYTBD"]AbstractFBCModels.CanonicalModel.Reaction(
name = "Cytochrome oxidase bd (ubiquinol-8: 2 protons)",
lower_bound = 0.0,
upper_bound = 1000.0,
stoichiometry = Dict("h2o_c" => 1.0, "o2_c" => -0.5, "h_e" => 2.0, "q8_c" => …
objective_coefficient = 0.0,
gene_association_dnf = [["b0733", "b0734"], ["b0978", "b0979"]],
annotations = Dict("metanetx.reaction" => ["MNXR97031"], "sbo" => ["SBO:00001…
notes = Dict("original_bigg_ids" => ["CYTBD"]),
)
To do the perturbation, we create a model of a strain which has mild issues with running the CYTBD reaction. We use deepcopy to completely avoid any reference sharing issues.
limited_ecoli = deepcopy(ecoli)
limited_ecoli.reactions["CYTBD"].upper_bound = 10.010.0Finding parsimonious solutions
Becuase we are interested in realistic flux distributions, we have to use an analysis method which gives one – in this case, the parsimonious FBA will do just right. For later comparison, we first get the optimal parsimonious flux distribution in the perturbed model:
solution = parsimonious_flux_balance_analysis(limited_ecoli, optimizer = Clarabel.Optimizer)ConstraintTrees.Tree{Float64} with 5 elements:
:coupling => ConstraintTrees.Tree{Float64}(#= 0 elements =#)
:flux_stoichiometry => ConstraintTrees.Tree{Float64}(#= 72 elements =#)
:fluxes => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
:objective => 0.391648
:parsimonious_objective => 5087.91Now, how much is the flux going to differ if we assume the bacterium did only minimal adjustment from the previous state with unlimited CYTBD?
moma_solution = metabolic_adjustment_minimization_analysis(
limited_ecoli, # the model to be examined
ecoli; # the model that gives the reference flux
optimizer = Clarabel.Optimizer,
)ConstraintTrees.Tree{Float64} with 5 elements:
:coupling => ConstraintTrees.Tree{Float64}(#= 0 elements …
:flux_stoichiometry => ConstraintTrees.Tree{Float64}(#= 72 elements…
:fluxes => ConstraintTrees.Tree{Float64}(#= 95 elements…
:minimal_adjustment_objective => 4937.56
:objective => 0.241497Comparing the results
The difference between the naive and minimally-adjusting solutions can be extracted using the constraint tree functionality:
import ConstraintTrees as C
difference = C.zip(-, solution, moma_solution, Float64)ConstraintTrees.Tree{Float64} with 4 elements:
:coupling => ConstraintTrees.Tree{Float64}(#= 0 elements =#)
:flux_stoichiometry => ConstraintTrees.Tree{Float64}(#= 72 elements =#)
:fluxes => ConstraintTrees.Tree{Float64}(#= 95 elements =#)
:objective => 0.150152sort(collect(difference.fluxes), by = last)95-element Vector{Pair{Symbol, Union{Float64, ConstraintTrees.Tree{Float64}}}}:
:FORt => -15.937900590627788
:EX_h2o_e => -12.259737504824756
:ACKr => -12.231123778244362
:ACt2r => -12.231123778242607
:PDH => -10.669672606180624
:EX_co2_e => -9.924850036344495
:ATPS4r => -8.655939529982444
:PPC => -4.689341153717746
:FRD7 => -4.073826731055466
:EX_etoh_e => -3.7271169823948425
⋮
:ACALD => 5.6514793889089265
:THD2 => 6.71524350812718
:CO2t => 9.924850036344582
:EX_ac_e => 12.231123778242313
:PTAr => 12.23112377825292
:H2Ot => 12.259737504825004
:EX_for_e => 15.937900592298162
:PFL => 15.937900592299522
:EX_h_e => 20.94183649468937Using a custom reference flux
In certain situations, one might want to examine how the model would adjust from a known reaction flux. We can supply it manually as the second argument (instead of the reference model).
ref = parsimonious_flux_balance_analysis(ecoli, optimizer = Clarabel.Optimizer)
ref_closest_solution = metabolic_adjustment_minimization_analysis(
limited_ecoli,
ref.fluxes;
optimizer = Clarabel.Optimizer,
)ConstraintTrees.Tree{Float64} with 5 elements:
:coupling => ConstraintTrees.Tree{Float64}(#= 0 elements …
:flux_stoichiometry => ConstraintTrees.Tree{Float64}(#= 72 elements…
:fluxes => ConstraintTrees.Tree{Float64}(#= 95 elements…
:minimal_adjustment_objective => 4937.56
:objective => 0.241497The flux may even be partial (which is common with measured fluxes):
measured_fluxes =
C.Tree{Float64}(:EX_ac_e => 5.0, :EX_o2_e => -2.0, :BIOMASS_Ecoli_core_w_GAM => 0.7)
solution_close_to_measurement = metabolic_adjustment_minimization_analysis(
limited_ecoli,
measured_fluxes;
optimizer = Clarabel.Optimizer,
)ConstraintTrees.Tree{Float64} with 5 elements:
:coupling => ConstraintTrees.Tree{Float64}(#= 0 elements …
:flux_stoichiometry => ConstraintTrees.Tree{Float64}(#= 72 elements…
:fluxes => ConstraintTrees.Tree{Float64}(#= 95 elements…
:minimal_adjustment_objective => 0.183188
:objective => 0.272247Efficient linear-metric MOMA
The linear version of MOMA avoids having to use the quadratic optimizer in the process, giving more optimizer choices and (typically) much better performance. Linear MOMA has the same interface as the quadratic one:
linear_moma_solution = linear_metabolic_adjustment_minimization_analysis(
limited_ecoli,
ecoli;
optimizer = HiGHS.Optimizer,
)
sort(collect(linear_moma_solution.fluxes), by = last)95-element Vector{Pair{Symbol, Union{Float64, ConstraintTrees.Tree{Float64}}}}:
:PGK => -15.434651933306435
:PGM => -15.106655803899905
:H2Ot => -11.54994054300062
:CO2t => -8.778386666898754
:EX_glc__D_e => -8.322012190986193
:ACALD => -5.755385239519781
:ALCD2x => -5.755385239519781
:ETOHt2r => -5.755385239519781
:EX_o2_e => -5.0
:EX_nh4_e => -3.664419001773459
⋮
:EX_co2_e => 8.778386666898754
:PDH => 9.282532599166617
:CYTBD => 10.0
:EX_h2o_e => 11.54994054300062
:NADH16 => 11.99204468413819
:ATPS4r => 12.020451381843891
:EX_h_e => 13.32001812519156
:ENO => 15.106655803899905
:GAPD => 15.434651933306435How much does the flux distribution differ from the L2 solution?
sort(
collect(C.zip(-, linear_moma_solution.fluxes, moma_solution.fluxes, Float64)),
by = last,
)95-element Vector{Pair{Symbol, Union{Float64, ConstraintTrees.Tree{Float64}}}}:
:GAPD => -2.217421474279119
:ENO => -2.184138541734182
:FRD7 => -2.0818742854994783
:EX_acald_e => -1.924362406530278
:EX_h_e => -1.7636376825709608
:PDH => -1.3871400070241506
:GLCpts => -1.1700704237445851
:SUCCt3 => -1.0914174611248049
:EX_succ_e => -1.0914174609191516
:NADH16 => -1.0914174606113036
⋮
:CO2t => 0.44670508165553535
:H2Ot => 1.046555418945502
:FUM => 1.0914174606558396
:MDH => 1.0914174620363917
:EX_glc__D_e => 1.170070423744587
:ACALD => 1.7125013840401522
:ACALDt => 1.9243624065302796
:PGM => 2.184138541734182
:PGK => 2.217421474279126This page was generated using Literate.jl.