Minimization of metabolic adjustment (MOMA)

MOMA allows you to find a feasible solution of the model that is closest (in an Euclidean metric) to a reference solution. Often this gives a realistic estimate of the organism behavior that has undergone a radical change (such as a gene knockout) that prevents it from metabolizing optimally, but the rest of the metabolism has not yet adjusted to compensate for the change.

The original description of MOMA is by Segre, Vitkup, and Church, "Analysis of optimality in natural and perturbed metabolic networks", Proceedings of the National Academy of Sciences, 2002.

As always, let's start with downloading a model.

!isfile("e_coli_core.xml") &&
    download("http://bigg.ucsd.edu/static/models/e_coli_core.xml", "e_coli_core.xml")

using COBREXA

model = load_model(StandardModel, "e_coli_core.xml")

Metabolic model of type StandardModel
sparse([8, 10, 21, 43, 50, 51, 8, 9, 6, 12  …  33, 66, 68, 72, 23, 26, 33, 72, 22, 33], [1, 1, 1, 1, 1, 1, 2, 2, 3, 3  …  93, 93, 93, 93, 94, 94, 94, 94, 95, 95], [-1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0  …  1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0], 72, 95)
Number of reactions: 95
Number of metabolites: 72

MOMA analysis requires solution of a quadratic model, we will thus use Clarabel as the main optimizer.

using Clarabel

We will need a reference solution, which represents the original state of the organism before the change.

reference_flux =
    flux_balance_analysis_dict(model, Clarabel.Optimizer; modifications = [silence])

Dict{String, Float64} with 95 entries:
  "R_EX_fum_e"    => -5.00659e-12
  "R_ACONTb"      => 6.00725
  "R_GLNS"        => 0.223462
  "R_SUCOAS"      => -5.06438
  "R_TPI"         => 7.47738
  "R_EX_pi_e"     => -3.2149
  "R_PPC"         => 2.50431
  "R_O2t"         => 21.7995
  "R_G6PDH2r"     => 4.95998
  "R_TALA"        => 1.49698
  "R_PPCK"        => 5.75426e-10
  "R_EX_lac__D_e" => 6.70413e-11
  "R_PGL"         => 4.95998
  "R_H2Ot"        => -29.1758
  "R_GLNabc"      => 2.95563e-12
  "R_EX_co2_e"    => 22.8098
  "R_EX_gln__L_e" => -2.95564e-12
  "R_EX_nh4_e"    => -4.76532
  "R_MALt2_2"     => 5.34064e-12
  ⋮               => ⋮

As the change here, we manually knock out CYTBD reaction:

changed_model = change_bound(model, "R_CYTBD", lower = 0.0, upper = 0.0);

Now, let's find a flux that minimizes the organism's metabolic adjustment for this model:

flux_summary(
    minimize_metabolic_adjustment_analysis_dict(
        changed_model,
        reference_flux,
        Clarabel.Optimizer;
        modifications = [silence],
    ),
)

Biomass
  R_BIOMASS_Ecoli_core_w_GAM: 0.06
Import
  R_EX_glc__D_e:              -9.1554
  R_EX_nh4_e:                 -1.3091
  R_EX_pi_e:                  -0.2209
Export
  R_EX_glu__L_e:              0.9817
  R_EX_acald_e:               1.7139
  R_EX_co2_e:                 5.1127
  R_EX_h2o_e:                 5.7936
  R_EX_succ_e:                5.9147
  R_EX_etoh_e:                7.6345
  R_EX_h_e:                   14.9974

For illustration, you can compare the result to the flux that is found by simple optimization:

flux_summary(
    flux_balance_analysis_dict(
        changed_model,
        Clarabel.Optimizer;
        modifications = [silence],
    ),
)

Biomass
  R_BIOMASS_Ecoli_core_w_GAM: 0.2117
Import
  R_EX_glc__D_e:              -10.0
  R_EX_h2o_e:                 -7.1158
  R_EX_nh4_e:                 -1.1542
  R_EX_pi_e:                  -0.7786
  R_EX_co2_e:                 -0.3782
Export
  R_EX_etoh_e:                8.2795
  R_EX_ac_e:                  8.5036
  R_EX_for_e:                 17.8047
  R_EX_h_e:                   30.5542

This page was generated using Literate.jl.