# 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.

!isfile("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"    => -6.7812e-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"        => 3.96203e-9
"R_EX_lac__D_e" => 5.58955e-10
"R_PGL"         => 4.95998
"R_H2Ot"        => -29.1758
"R_GLNabc"      => 2.25337e-12
"R_EX_co2_e"    => 22.8098
"R_EX_gln__L_e" => -2.25168e-12
"R_EX_nh4_e"    => -4.76532
"R_MALt2_2"     => 1.11042e-11
⋮               => ⋮

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(
changed_model,
reference_flux,
Clarabel.Optimizer;
modifications = [silence],
),
)
Biomass
R_BIOMASS_Ecoli_core_w_GAM: 0.06
Import
R_EX_glc__D_e:              -9.1553
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.1126
R_EX_h2o_e:                 5.7937
R_EX_succ_e:                5.9147
R_EX_etoh_e:                7.6345
R_EX_h_e:                   14.9973


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