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.937900590627708
  :EX_h2o_e => -12.25973750482275
      :ACKr => -12.23112377824434
     :ACt2r => -12.231123778242587
       :PDH => -10.669672606178668
  :EX_co2_e => -9.924850036343367
    :ATPS4r => -8.655939529984154
       :PPC => -4.689341153719435
      :FRD7 => -4.073826730909896
 :EX_etoh_e => -3.727116982395617
            ⋮
     :ACALD => 5.6514793889087
      :THD2 => 6.71524350812716
      :CO2t => 9.924850036343452
   :EX_ac_e => 12.231123778242294
      :PTAr => 12.231123778252899
      :H2Ot => 12.259737504822999
  :EX_for_e => 15.937900592298138
       :PFL => 15.937900592299496
    :EX_h_e => 20.94183649468959

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.217421474277627
         :ENO => -2.1841385417330486
        :FRD7 => -2.0818742855029964
  :EX_acald_e => -1.9243624065292773
      :EX_h_e => -1.7636376825707192
         :PDH => -1.3871400070221949
      :GLCpts => -1.1700704237440558
      :SUCCt3 => -1.0914174611276102
   :EX_succ_e => -1.0914174609219507
      :NADH16 => -1.0914174606140925
              ⋮
        :CO2t => 0.4467050816544056
        :H2Ot => 1.046555418943509
         :FUM => 1.0914174606586302
         :MDH => 1.0914174620392267
 :EX_glc__D_e => 1.1700704237440576
       :ACALD => 1.7125013840399266
      :ACALDt => 1.9243624065292788
         :PGM => 2.184138541733045
         :PGK => 2.2174214742776375

This page was generated using Literate.jl.