Flux variability analysis (FVA)
Here we will use flux_variability_analysis to analyze the E. coli core model.
As usual, if not already present, download the model and load the required packages. We picked the GLPK solver, but others may work as well:
!isfile("e_coli_core.xml") &&
download("http://bigg.ucsd.edu/static/models/e_coli_core.xml", "e_coli_core.xml")
using COBREXA, GLPK
model = load_model("e_coli_core.xml")Metabolic model of type SBMLModel 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
The FVA implementation in flux_variability_analysis returns maximized and minimized reaction fluxes in a 2-column matrix. The bounds parameter function here (constructed with objective_bounds) sets that the objective value is allowed to vary by 1% from the optimum found by FBA on the same model:
flux_variability_analysis(model, GLPK.Optimizer; bounds = objective_bounds(0.99))95×2 Matrix{Float64}:
-0.254237 0.0
-0.254237 0.0
-0.381356 0.0
4.8143 7.56006
4.8143 7.56006
-0.381356 0.0
-1.11022e-15 1.7161
1.18233 6.62661
-0.143008 0.0
-0.221432 0.0
⋮
0.0 2.28813
1.84464e-15 2.28813
3.88085 1000.0
-6.62661 -1.18233
0.0691339 2.81489
-8.88178e-16 6.8644
0.0691339 2.81489
-0.243197 2.50256
6.16973 8.91549(You may also use gamma_bounds.)
Detailed variability analysis with modifications
A dictionary-returning variant in flux_variability_analysis_dict, returns the result in a slightly more structured way. At the same time, we can specify additional modifications to be applied to the model:
min_fluxes, max_fluxes = flux_variability_analysis_dict(
model,
GLPK.Optimizer;
bounds = objective_bounds(0.99),
modifications = [
change_constraint("R_EX_glc__D_e"; lb = -10, ub = -10),
change_constraint("R_EX_o2_e"; lb = 0.0, ub = 0.0),
],
)(Dict("R_EX_fum_e" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.22607952490472283, "R_GLNS" => 0.053580994084844846, "R_SUCOAS" => 6.225897942841894e-15, "R_TPI" => 9.791564275259212, "R_EX_pi_e" => -0.7708580482593426, "R_PPC" => 0.6004759352738883, "R_O2t" => 1.7763568394002587e-15, "R_G6PDH2r" => 2.438809336803311e-15, "R_TALA" => -0.03748783669056421…), "R_ACONTb" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.22607952490472283, "R_GLNS" => 0.053580994084844624, "R_SUCOAS" => 9.383674941524645e-15, "R_TPI" => 9.791564275259212, "R_EX_pi_e" => -0.7708580482593473, "R_PPC" => 0.6004759352738988, "R_O2t" => 2.6645352591003883e-15, "R_G6PDH2r" => 2.1553062906483042e-15, "R_TALA" => -0.03748783669056422…), "R_GLNS" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.22607952490472882, "R_GLNS" => 0.053580994084844895, "R_SUCOAS" => 0.0, "R_TPI" => 9.485804275258849, "R_EX_pi_e" => -0.7708580482593624, "R_PPC" => 0.6004759352738828, "R_O2t" => 0.0, "R_G6PDH2r" => 0.9172800000013737, "R_TALA" => 0.2682721633098879…), "R_SUCOAS" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.29557043399574195, "R_GLNS" => 0.05358099408484489, "R_SUCOAS" => -0.06949090909101316, "R_TPI" => 9.791564275259306, "R_EX_pi_e" => -0.7708580482593625, "R_PPC" => 0.669966844364897, "R_O2t" => 0.0, "R_G6PDH2r" => -9.39717574588439e-17, "R_TALA" => -0.03748783669057004…), "R_TPI" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.22607952490472882, "R_GLNS" => 0.053580994084844895, "R_SUCOAS" => 0.0, "R_TPI" => 9.485804275258849, "R_EX_pi_e" => -0.7708580482593624, "R_PPC" => 0.6004759352738828, "R_O2t" => 0.0, "R_G6PDH2r" => 0.9172800000013737, "R_TALA" => 0.2682721633098879…), "R_EX_pi_e" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.22836315646940947, "R_GLNS" => 0.054122216247316675, "R_SUCOAS" => 3.874791485706844e-15, "R_TPI" => 9.789458863898215, "R_EX_pi_e" => -0.7786444931912458, "R_PPC" => 0.6065413487614906, "R_O2t" => 1.7763568394002587e-15, "R_G6PDH2r" => -4.768316294366138e-15, "R_TALA" => -0.03786650170764047…), "R_PPC" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.2872315249048204, "R_GLNS" => 0.05358099408484489, "R_SUCOAS" => 0.0, "R_TPI" => 9.791564275259306, "R_EX_pi_e" => -0.7708580482593624, "R_PPC" => 0.6004759352738829, "R_O2t" => 0.0, "R_G6PDH2r" => -9.39717574588439e-17, "R_TALA" => -0.03748783669057004…), "R_O2t" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.29254909012221963, "R_GLNS" => 0.05358099408484489, "R_SUCOAS" => 0.0, "R_TPI" => 9.791564275259306, "R_EX_pi_e" => -0.7708580482593625, "R_PPC" => 0.6669455004913749, "R_O2t" => 0.0, "R_G6PDH2r" => -1.879435149176878e-16, "R_TALA" => -0.03748783669057007…), "R_G6PDH2r" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.22607952490472105, "R_GLNS" => 0.05358099408484448, "R_SUCOAS" => 4.9960036108132044e-15, "R_TPI" => 9.791564275259217, "R_EX_pi_e" => -0.7708580482593473, "R_PPC" => 0.7915759352734985, "R_O2t" => -1.7763568394002587e-15, "R_G6PDH2r" => -6.240482714761413e-15, "R_TALA" => -0.03748783669056543…), "R_TALA" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.22836315646942648, "R_GLNS" => 0.05412221624731889, "R_SUCOAS" => 0.0, "R_TPI" => 9.789458863898286, "R_EX_pi_e" => -0.778644493191287, "R_PPC" => 0.6065413487615077, "R_O2t" => 0.0, "R_G6PDH2r" => 0.0, "R_TALA" => -0.03786650170764705…)…), Dict("R_EX_fum_e" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.28723152490482035, "R_GLNS" => 0.05358099408484489, "R_SUCOAS" => 0.0, "R_TPI" => 9.791564275259306, "R_EX_pi_e" => -0.7708580482593624, "R_PPC" => 0.6004759352738829, "R_O2t" => 0.0, "R_G6PDH2r" => 0.0, "R_TALA" => -0.037487836690570014…), "R_ACONTb" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.29557043399574195, "R_GLNS" => 0.05358099408484489, "R_SUCOAS" => -0.06949090909101316, "R_TPI" => 9.791564275259306, "R_EX_pi_e" => -0.7708580482593625, "R_PPC" => 0.669966844364897, "R_O2t" => 0.0, "R_G6PDH2r" => -9.39717574588439e-17, "R_TALA" => -0.03748783669057004…), "R_GLNS" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.22607952490472868, "R_GLNS" => 0.24468099408513094, "R_SUCOAS" => 0.0, "R_TPI" => 9.791564275259308, "R_EX_pi_e" => -0.7708580482593625, "R_PPC" => 0.6004759352738851, "R_O2t" => 0.0, "R_G6PDH2r" => 2.664535259100375e-16, "R_TALA" => -0.0374878366905699…), "R_SUCOAS" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.22607952490472877, "R_GLNS" => 0.05358099408484488, "R_SUCOAS" => 0.0, "R_TPI" => 9.791564275259306, "R_EX_pi_e" => -0.7708580482593625, "R_PPC" => 0.600475935273883, "R_O2t" => 0.0, "R_G6PDH2r" => 0.0, "R_TALA" => -0.03748783669056999…), "R_TPI" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.28977952490482417, "R_GLNS" => 0.05358099408484488, "R_SUCOAS" => 0.0, "R_TPI" => 9.791564275259306, "R_EX_pi_e" => -0.7708580482593625, "R_PPC" => 0.6641759352739797, "R_O2t" => 0.0, "R_G6PDH2r" => 0.0, "R_TALA" => -0.03748783669056999…), "R_EX_pi_e" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.22607952490472874, "R_GLNS" => 0.05358099408484488, "R_SUCOAS" => 0.0, "R_TPI" => 9.791564275259306, "R_EX_pi_e" => -0.7708580482593625, "R_PPC" => 0.6004759352738848, "R_O2t" => 1.387778780781446e-17, "R_G6PDH2r" => 0.0, "R_TALA" => -0.03748783669056999…), "R_PPC" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.22607952490472877, "R_GLNS" => 0.05358099408484488, "R_SUCOAS" => 0.0, "R_TPI" => 9.791564275259306, "R_EX_pi_e" => -0.7708580482593625, "R_PPC" => 0.982675935274456, "R_O2t" => 0.0, "R_G6PDH2r" => 0.0, "R_TALA" => -0.03748783669056999…), "R_O2t" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.2897795249048397, "R_GLNS" => 0.053580994084844895, "R_SUCOAS" => 0.0, "R_TPI" => 9.79156427525938, "R_EX_pi_e" => -0.7708580482593625, "R_PPC" => 0.6641759352739937, "R_O2t" => 0.0, "R_G6PDH2r" => -2.2302160118670142e-13, "R_TALA" => -0.03748783669064436…), "R_G6PDH2r" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.22607952490472882, "R_GLNS" => 0.05358099408484489, "R_SUCOAS" => 0.0, "R_TPI" => 9.485804275258849, "R_EX_pi_e" => -0.7708580482593624, "R_PPC" => 0.6004759352738829, "R_O2t" => 0.0, "R_G6PDH2r" => 0.9172800000013739, "R_TALA" => 0.26827216330988796…), "R_TALA" => Dict("R_EX_fum_e" => 0.0, "R_ACONTb" => 0.22607952490472877, "R_GLNS" => 0.05358099408484488, "R_SUCOAS" => 0.0, "R_TPI" => 9.485804275258849, "R_EX_pi_e" => -0.7708580482593624, "R_PPC" => 0.6004759352738848, "R_O2t" => 0.0, "R_G6PDH2r" => 0.9172800000013748, "R_TALA" => 0.2682721633098883…)…))The dictionaries can be easily used to explore the whole state of the model when certain reactions are maximized or minimized. For example, we can take the maximal acetate exchange flux when the acetate exchange is maximized:
max_fluxes["R_EX_ac_e"]["R_EX_ac_e"]8.51854942518171
We can also check that the modifications really had the desired effect on oxygen consumption:
max_fluxes["R_EX_ac_e"]["R_O2t"]0.0
...and see how much carbon dioxide would produced under at the given metabolic extreme:
max_fluxes["R_EX_ac_e"]["R_EX_co2_e"]-0.5654964103694404
Summarizing the flux variability
A convenience function flux_variability_summary is able to display this information in a nice overview:
flux_variability_summary((min_fluxes, max_fluxes))Biomass Lower bound Upper bound R_BIOMASS_Ecoli_core_w_GAM: 0.2095 0.2095 Exchange R_EX_h_e: 28.2555 30.7398 R_EX_for_e: 16.2978 17.8266 R_EX_glc__D_e: -10.0 -10.0 R_EX_etoh_e: 8.0419 9.0611 R_EX_ac_e: 7.4484 8.5185 R_EX_h2o_e: -7.1446 -6.3802 R_EX_nh4_e: -1.2063 -1.1426 R_EX_co2_e: -0.5655 1.1544 R_EX_pi_e: -0.7786 -0.7709 R_EX_lac__D_e: 0.0 0.5096 R_EX_acald_e: 0.0 0.3058 R_EX_pyr_e: 0.0 0.2548 R_EX_succ_e: -0.0 0.1911 R_EX_akg_e: 0.0 0.0665 R_EX_glu__L_e: -0.0 0.0637 R_EX_o2_e: -0.0 -0.0 R_EX_fum_e: 0.0 0.0 R_EX_mal__L_e: 0.0 0.0 R_EX_gln__L_e: 0.0 0.0 R_EX_fru_e: 0.0 0.0
Retrieving details about FVA output
Parameter ret of flux_variability_analysis can be used to extract specific pieces of information from the individual solved (minimized and maximized) optimization problems. Here we show how to extract the value of biomass "growth" along with the minimized/maximized reaction flux.
First, find the index of biomass reaction in all reactions
biomass_idx = first(indexin(["R_BIOMASS_Ecoli_core_w_GAM"], reactions(model)))13
Now run the FVA:
vs = flux_variability_analysis(
model,
GLPK.Optimizer;
bounds = objective_bounds(0.50), # objective can vary by up to 50% of the optimum
modifications = [
change_constraint("R_EX_glc__D_e"; lb = -10, ub = -10),
change_constraint("R_EX_o2_e"; lb = 0.0, ub = 0.0),
],
ret = optimized_model -> (
COBREXA.JuMP.objective_value(optimized_model),
COBREXA.JuMP.value(optimized_model[:x][biomass_idx]),
),
)95×2 Matrix{Tuple{Float64, Float64}}:
(-18.3915, 0.105831) (-1.77636e-15, 0.105831)
(-15.288, 0.105831) (0.0, 0.105831)
(-9.25179, 0.105831) (0.0, 0.105831)
(0.114182, 0.105831) (3.4504, 0.105831)
(0.114182, 0.105831) (3.4504, 0.105831)
(-9.25179, 0.105831) (0.0, 0.105831)
(7.18869e-15, 0.105831) (9.555, 0.105831)
(3.48454e-15, 0.105831) (3.08393, 0.105831)
(-3.09994, 0.105831) (0.0, 0.105831)
(-18.3915, 0.105831) (9.54792e-15, 0.105831)
⋮
(0.0, 0.105831) (12.74, 0.105831)
(0.0, 0.105831) (12.74, 0.105831)
(0.0, 0.105831) (1000.0, 0.105831)
(-3.08393, 0.105831) (4.44089e-16, 0.105831)
(-0.0378665, 0.211663) (2.62444, 0.105831)
(0.0, 0.207557) (20.9246, 0.105831)
(-0.0378665, 0.211663) (2.62444, 0.105831)
(-0.114277, 0.211663) (2.58623, 0.105831)
(7.25136, 0.105831) (9.89473, 0.105831)This page was generated using Literate.jl.