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, HiGHS
Because 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.0
10.0
Finding 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.91
Now, 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.241497
Comparing 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.150152
sort(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.94183649468937
Using 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.241497
The 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.272247
Efficient 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.434651933306435
How 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.217421474279126
This page was generated using Literate.jl.