Finding balance and variability of constraint-based models
Here we will use flux_balance_analysis, flux_variability_analysis, parsimonious_flux_balance_analysis, and minimize_metabolic_adjustment_analysis, along with the modification functions of COBREXA.jl, to analyze a toy model of E. coli.
If it is not already present, download the model.
!isfile("e_coli_core.xml") &&
download("http://bigg.ucsd.edu/static/models/e_coli_core.xml", "e_coli_core.xml")
using COBREXAWhen you are unsure about how a function works, write ? function_name to see the function reference documentation.
model = load_model("e_coli_core.xml")Metabolic model of type SBMLModel sparse([41, 23, 51, 67, 61, 65, 1, 7, 19, 28 … 72, 3, 8, 33, 57, 66, 31, 45, 46, 57], [1, 2, 2, 2, 3, 3, 4, 4, 4, 4 … 93, 94, 94, 94, 94, 94, 95, 95, 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
Optimization solvers in COBREXA
To actually perform any optimization based analysis we need to load an optimizer. Any JuMP.jl-supported optimizers will work. Here, we will use Tulip.jl to optimize linear programs and OSQP.jl to optimize quadratic programs.
We recommend reading the docs of OSQP before using it, since it may give inconsistent results depending on what settings you use. Commercial solvers like Gurobi, Mosek, CPLEX, etc. require less user engagement.
using Tulip, OSQP, GLPKFlux balance analysis (FBA)
Most analysis functions come in several variants that produce different types of output. All of them usually require a model and JuMP.jl-compatible optimizer to work in the model.
In the case of FBA, you may choose from these variants (here using the Tulip optimizer):
vec_soln = flux_balance_analysis_vec(model, Tulip.Optimizer)95-element Vector{Float64}:
-0.0
6.00724956649032
7.477381918907127
-5.064375360152338
0.2234617471432185
-3.214895030387032
2.504309432010867
21.799492758475754
4.959985078874371
1.496983802869297
⋮
3.375438217960911e-7
29.175827202685298
9.054357964341115e-9
4.817965631705414e-8
9.959461594581987e-9
-21.799492758475754
-0.0
-1.4340676616267298e-9
3.214895030387032dict_soln = flux_balance_analysis_dict(model, Tulip.Optimizer)Dict{String, Float64} with 95 entries:
"R_EX_fum_e" => -0.0
"R_ACONTb" => 6.00725
"R_TPI" => 7.47738
"R_SUCOAS" => -5.06438
"R_GLNS" => 0.223462
"R_EX_pi_e" => -3.2149
"R_PPC" => 2.50431
"R_O2t" => 21.7995
"R_G6PDH2r" => 4.95999
"R_TALA" => 1.49698
"R_PPCK" => 5.88317e-8
"R_EX_lac__D_e" => 2.39394e-9
"R_PGL" => 4.95999
"R_H2Ot" => -29.1758
"R_GLNabc" => -0.0
"R_EX_co2_e" => 22.8098
"R_EX_gln__L_e" => -0.0
"R_EX_nh4_e" => -4.76532
"R_MALt2_2" => -0.0
⋮ => ⋮Extending FBA with modifications
Often it is desirable to add a slight modification to the problem before performing analysis, to see e.g. differences of the model behavior caused by the change introduced.
COBREXA.jl supports several modifications by default, which include changing objective sense, optimizer attributes, flux constraints, optimization objective, reaction and gene knockouts, and others. These modifications are applied in the order they are specified. It is up to the user to ensure that the changes are sensible.
dict_soln = flux_balance_analysis_dict(
model,
OSQP.Optimizer;
modifications = [ # modifications are applied in order
# this changes the objective to maximize the biomass production
change_objective("R_BIOMASS_Ecoli_core_w_GAM"),
# this fixes a specific rate of the glucose exchange
change_constraint("R_EX_glc__D_e"; lb = -12, ub = -12),
# this knocks out two genes, i.e. constrains their associated reactions to zero.
knockout(["b0978", "b0734"]), ## the gene IDs are cytochrome oxidase (CYTBD)
# ignore the optimizer specified above and change it to Tulip
change_optimizer(Tulip.Optimizer),
# set a custom attribute of the Tulip optimizer (see Tulip docs for more possibilities)
change_optimizer_attribute("IPM_IterationsLimit", 110),
# explicitly tell the optimizer to maximize the new objective
change_sense(MAX_SENSE),
],
)Dict{String, Float64} with 95 entries:
"R_EX_fum_e" => -0.0
"R_ACONTb" => 7.03277
"R_TPI" => 8.90908
"R_SUCOAS" => -5.8921
"R_GLNS" => 0.270339
"R_EX_pi_e" => -3.88931
"R_PPC" => 3.02966
"R_O2t" => 25.7859
"R_G6PDH2r" => 6.11782
"R_TALA" => 1.85013
"R_PPCK" => 5.26409e-10
"R_EX_lac__D_e" => 4.37341e-12
"R_PGL" => 6.11782
"R_H2Ot" => -34.7096
"R_GLNabc" => -0.0
"R_EX_co2_e" => 27.0082
"R_EX_gln__L_e" => -0.0
"R_EX_nh4_e" => -5.76498
"R_MALt2_2" => -0.0
⋮ => ⋮This solution can be display using flux_summary. Note, this pretty printing only works on flux solutions that are represented as dictionaries.
flux_summary(dict_soln)Biomass R_BIOMASS_Ecoli_core_w_GAM: 1.0573 Import R_EX_o2_e: -25.7859 R_EX_glc__D_e: -12.0 R_EX_nh4_e: -5.765 R_EX_pi_e: -3.8893 Export R_EX_h_e: 21.2085 R_EX_co2_e: 27.0082 R_EX_h2o_e: 34.7096
Flux variability analysis (FVA)
The default FVA in flux_variability_analysis returns maximized and minimized reaction fluxes in a matrix. Here we use the dictionary variant in fluxvariabilityanalysis_dict, to show how to easily access specific fluxes from its results.
fva_mins, fva_maxs = flux_variability_analysis_dict(
model,
Tulip.Optimizer;
bounds = objective_bounds(0.99), # the objective function is allowed to vary by ~1% from the FBA optimum
modifications = [
change_optimizer_attribute("IPM_IterationsLimit", 500),
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.23975681252388875, "R_TPI" => 9.758987670019215, "R_SUCOAS" => -0.003217695420622425, "R_GLNS" => 0.06020057742388561, "R_EX_pi_e" => -0.7710723486299674, "R_PPC" => 0.6485114654496659, "R_O2t" => 4.7636840466449975e-17, "R_G6PDH2r" => 0.09755597883871935, "R_TALA" => -0.004979598785419143…), "R_ACONTb" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952399024103, "R_TPI" => 9.751314994785577, "R_SUCOAS" => -5.747711214668542e-11, "R_GLNS" => 0.062181323693887736, "R_EX_pi_e" => -0.7708580440960587, "R_PPC" => 0.6578792778328494, "R_O2t" => 9.411084212042624e-15, "R_G6PDH2r" => 0.12074784479375521, "R_TALA" => 0.0027614451094113016…), "R_TPI" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2260795243453041, "R_TPI" => 9.485804104877825, "R_SUCOAS" => -8.319304757951892e-11, "R_GLNS" => 0.05358099493359015, "R_EX_pi_e" => -0.7708580445453711, "R_PPC" => 0.6004759354091982, "R_O2t" => 5.359934820339382e-15, "R_G6PDH2r" => 0.9172805136407763, "R_TALA" => 0.2682723346609655…), "R_SUCOAS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.29557047715224044, "R_TPI" => 9.791564276282955, "R_SUCOAS" => -0.06949095362117334, "R_GLNS" => 0.05358099381978238, "R_EX_pi_e" => -0.7708580432292014, "R_PPC" => 0.6699668853159306, "R_O2t" => 3.191211084827664e-15, "R_G6PDH2r" => 9.982382324749958e-10, "R_TALA" => -0.037487836114279326…), "R_GLNS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.24064242596917626, "R_TPI" => 9.75708160792236, "R_SUCOAS" => -0.0034442975170188377, "R_GLNS" => 0.0535809937394168, "R_EX_pi_e" => -0.7708580432437464, "R_PPC" => 0.6515935674935543, "R_O2t" => 2.1424377258038788e-16, "R_G6PDH2r" => 0.10344800607932744, "R_TALA" => -0.003005167753552366…), "R_EX_pi_e" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22836315673272758, "R_TPI" => 9.789458863462334, "R_SUCOAS" => -7.227058044154967e-11, "R_GLNS" => 0.05412221657624835, "R_EX_pi_e" => -0.7786444930369882, "R_PPC" => 0.606541349823176, "R_O2t" => 4.510775788338821e-15, "R_G6PDH2r" => 1.429168215925255e-9, "R_TALA" => -0.03786650122407263…), "R_PPC" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2485481901699732, "R_TPI" => 9.758144252473867, "R_SUCOAS" => -1.75099877155713e-13, "R_GLNS" => 0.0607114201940413, "R_EX_pi_e" => -0.770858043206322, "R_PPC" => 0.6004759313481519, "R_O2t" => 1.0686845246346513e-16, "R_G6PDH2r" => 0.10026007245513789, "R_TALA" => -0.00406781229313153…), "R_O2t" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.23975476417384725, "R_TPI" => 9.758984092418991, "R_SUCOAS" => -0.0032170588900653063, "R_GLNS" => 0.06020006640693463, "R_EX_pi_e" => -0.771072334219399, "R_PPC" => 0.6485122872632033, "R_O2t" => 4.902121164864372e-17, "R_G6PDH2r" => 0.09756672332899578, "R_TALA" => -0.004976016587856012…), "R_G6PDH2r" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.24130980459385273, "R_TPI" => 9.791500043279633, "R_SUCOAS" => -0.0035854151174607306, "R_GLNS" => 0.0612127234286182, "R_EX_pi_e" => -0.7710955974418061, "R_PPC" => 0.6545319530113679, "R_O2t" => 2.155552554181346e-16, "R_G6PDH2r" => 9.283759662954366e-12, "R_TALA" => -0.037499389015432455…), "R_TALA" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22836318744775283, "R_TPI" => 9.789458867641438, "R_SUCOAS" => -8.182361733333499e-9, "R_GLNS" => 0.05412221663518293, "R_EX_pi_e" => -0.7786444672239631, "R_PPC" => 0.6065415024286762, "R_O2t" => 4.108935134143224e-14, "R_G6PDH2r" => 9.828422042686805e-9, "R_TALA" => -0.03786649716850727…)…), Dict("R_EX_fum_e" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.23975681252388875, "R_TPI" => 9.758987670019215, "R_SUCOAS" => -0.003217695420622425, "R_GLNS" => 0.06020057742388561, "R_EX_pi_e" => -0.7710723486299674, "R_PPC" => 0.6485114654496659, "R_O2t" => 4.7636840466449975e-17, "R_G6PDH2r" => 0.09755597883871935, "R_TALA" => -0.004979598785419143…), "R_ACONTb" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.29557047721975777, "R_TPI" => 9.79156427645274, "R_SUCOAS" => -0.06949095088621661, "R_GLNS" => 0.0535809938024699, "R_EX_pi_e" => -0.7708580432015065, "R_PPC" => 0.669966885055497, "R_O2t" => 5.486413767863138e-15, "R_G6PDH2r" => 5.199963828733239e-10, "R_TALA" => -0.03748783627143602…), "R_TPI" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.24177375362265355, "R_TPI" => 9.791564276596777, "R_SUCOAS" => -0.003713574820679652, "R_GLNS" => 0.061441374870044536, "R_EX_pi_e" => -0.7708580432384676, "R_PPC" => 0.6559286179163718, "R_O2t" => 3.505736433264738e-16, "R_G6PDH2r" => 6.015355388027829e-11, "R_TALA" => -0.037487836426367944…), "R_SUCOAS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.23694450651847213, "R_TPI" => 9.757854115958603, "R_SUCOAS" => -4.214391845861122e-12, "R_GLNS" => 0.060469177452495845, "R_EX_pi_e" => -0.7710794678567561, "R_PPC" => 0.6491017899927606, "R_O2t" => 2.0205854785849238e-16, "R_G6PDH2r" => 0.10095086602264308, "R_TALA" => -0.0038483159414384595…), "R_GLNS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952365568096, "R_TPI" => 9.791564276059372, "R_SUCOAS" => -5.4667871707558614e-11, "R_GLNS" => 0.2446811147673207, "R_EX_pi_e" => -0.7708580432075332, "R_PPC" => 0.6004759321858603, "R_O2t" => 1.2603107594234158e-14, "R_G6PDH2r" => 1.6925839005206303e-9, "R_TALA" => -0.037487835881007216…), "R_EX_pi_e" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.24017982501148905, "R_TPI" => 9.758339870110252, "R_SUCOAS" => -0.003337568935484476, "R_GLNS" => 0.06044007540384528, "R_EX_pi_e" => -0.7708580432135304, "R_PPC" => 0.6496831008287483, "R_O2t" => 4.734147830827286e-16, "R_G6PDH2r" => 0.09967321954006271, "R_TALA" => -0.004263429931844625…), "R_PPC" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2260795234766408, "R_TPI" => 9.791564276535304, "R_SUCOAS" => -1.2654557928695605e-11, "R_GLNS" => 0.05358099376104628, "R_EX_pi_e" => -0.7708580432217733, "R_PPC" => 0.9826761788053896, "R_O2t" => -7.447072296632647e-16, "R_G6PDH2r" => 2.491419059931445e-10, "R_TALA" => -0.03748783636349198…), "R_O2t" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.23975885956391524, "R_TPI" => 9.75899124531553, "R_SUCOAS" => -0.0032183316533536104, "R_GLNS" => 0.06020108774390654, "R_EX_pi_e" => -0.7710723630183868, "R_PPC" => 0.6485106412140964, "R_O2t" => 4.6256792172025116e-17, "R_G6PDH2r" => 0.09754524127813421, "R_TALA" => -0.004983178672009124…), "R_G6PDH2r" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952341983936, "R_TPI" => 9.48580407775266, "R_SUCOAS" => -9.208178288767948e-14, "R_GLNS" => 0.05358099373346173, "R_EX_pi_e" => -0.7708580431951753, "R_PPC" => 0.6004759313305169, "R_O2t" => 2.4880501477390976e-17, "R_G6PDH2r" => 0.9172805966276716, "R_TALA" => 0.26827236243157754…), "R_TALA" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2260795234198597, "R_TPI" => 9.485804077755438, "R_SUCOAS" => -3.829553161833575e-14, "R_GLNS" => 0.053580993733406364, "R_EX_pi_e" => -0.7708580431954954, "R_PPC" => 0.6004759313306914, "R_O2t" => 1.978575548982736e-18, "R_G6PDH2r" => 0.9172805966189224, "R_TALA" => 0.2682723624286374…)…))fva_maxs["R_EX_ac_e"]["R_EX_ac_e"] # get the maximal acetate exchange flux8.518549434877523
Another option is to display this information using flux_variability_summary. This pretty printing only works on flux variability analysis results where dictionary keys indicate which flux is optimized and the associated value is a flux dictionary.
flux_variability_summary((fva_mins, fva_maxs))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_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_o2_e: 0.0 0.0 R_EX_fru_e: -0.0 -0.0
More sophisticated variants of flux_variability_analysis can be used to extract specific pieces of information from the solved optimization problems. Here the objective value of the minimized flux and the associated biomass growth rate is returned instead of every flux.
biomass_idx = first(indexin(["R_BIOMASS_Ecoli_core_w_GAM"], reactions(model))) # index of biomass function
vs = flux_variability_analysis(
model,
Tulip.Optimizer;
bounds = objective_bounds(0.50), # biomass can vary up to 50% less than optimum
modifications = [
change_optimizer_attribute("IPM_IterationsLimit", 500),
change_constraint("R_EX_glc__D_e"; lb = -10, ub = -10),
change_constraint("R_EX_o2_e"; lb = 0.0, ub = 0.0),
],
ret = m ->
(COBREXA.JuMP.objective_value(m), COBREXA.JuMP.value(m[:x][biomass_idx])), # m is the model and m[:x] extracts the fluxes from the model
)95×2 Matrix{Tuple{Float64, Float64}}:
(0.0, 0.112692) (-0.0, 0.112692)
(0.114182, 0.105831) (3.4504, 0.105831)
(7.25136, 0.105831) (9.89473, 0.105831)
(-3.08393, 0.105831) (-1.29494e-11, 0.112787)
(0.0270611, 0.105831) (9.58206, 0.105831)
(-0.778644, 0.211663) (-0.389322, 0.105831)
(0.303271, 0.105831) (10.7656, 0.105831)
(-5.4901e-18, 0.112692) (-5.53601e-18, 0.112692)
(4.67779e-13, 0.112782) (7.93011, 0.105831)
(-0.0378665, 0.211663) (2.62444, 0.105831)
⋮
(3.01295e-11, 0.112701) (20.9246, 0.105831)
(-8.5579, 0.105831) (5.65485, 0.105831)
(3.50341e-11, 0.112866) (9.555, 0.105831)
(1.5999e-11, 0.112882) (9.37857, 0.105831)
(6.05062e-11, 0.112888) (9.555, 0.105831)
(0.0, 0.112692) (-0.0, 0.112692)
(0.0, 0.112692) (-0.0, 0.112692)
(-18.3915, 0.105831) (-1.59455e-11, 0.107865)
(0.389322, 0.105831) (0.778644, 0.211663)fva_mins = Dict(rxn => flux for (rxn, flux) in zip(reactions(model), vs[:, 1]))Dict{String, Tuple{Float64, Float64}} with 95 entries:
"R_EX_fum_e" => (0.0, 0.112692)
"R_ACONTb" => (0.114182, 0.105831)
"R_TPI" => (7.25136, 0.105831)
"R_SUCOAS" => (-3.08393, 0.105831)
"R_GLNS" => (0.0270611, 0.105831)
"R_EX_pi_e" => (-0.778644, 0.211663)
"R_PPC" => (0.303271, 0.105831)
"R_O2t" => (-5.4901e-18, 0.112692)
"R_G6PDH2r" => (4.67779e-13, 0.112782)
"R_TALA" => (-0.0378665, 0.211663)
"R_PPCK" => (1.41777e-11, 0.112871)
"R_EX_lac__D_e" => (6.27486e-11, 0.11313)
"R_PGL" => (4.67543e-13, 0.112782)
"R_H2Ot" => (-5.65485, 0.105831)
"R_GLNabc" => (0.0, 0.112692)
"R_EX_co2_e" => (-9.38627, 0.105831)
"R_EX_gln__L_e" => (0.0, 0.112692)
"R_EX_nh4_e" => (-3.76208, 0.105831)
"R_MALt2_2" => (0.0, 0.112692)
⋮ => ⋮Parsimonious flux balance analysis (pFBA)
Parsimonious flux balance analysis (here in parsimonious_flux_balance_analysis finds a unique flux solution that minimizes the squared sum of fluxes of the system subject, while maintaining the same objective value as the flux balance analysis solution. Since we are optimizing a quadratic objective, we also need to switch to a quadratic optimizer. In this case, OSQP will work. We demonstrate it on the dictionary-returning variant of pFBA, parsimonious_flux_balance_analysis_dict:
dict_soln = parsimonious_flux_balance_analysis_dict(
model,
OSQP.Optimizer;
modifications = [
silence, # silence the optimizer (OSQP is very verbose by default)
change_constraint("R_EX_glc__D_e"; lb = -12, ub = -12),
change_optimizer_attribute("polish", true),
],
)Dict{String, Float64} with 95 entries:
"R_EX_fum_e" => -0.0054306
"R_ACONTb" => 6.51108
"R_TPI" => 8.90211
"R_SUCOAS" => -5.41568
"R_GLNS" => 0.250914
"R_EX_pi_e" => -3.85013
"R_PPC" => 2.94799
"R_O2t" => 25.1823
"R_G6PDH2r" => 6.27109
"R_TALA" => 1.90314
"R_PPCK" => -0.00186526
"R_EX_lac__D_e" => -0.00399055
"R_PGL" => 6.27111
"R_H2Ot" => -33.9472
"R_GLNabc" => 0.0126618
"R_EX_co2_e" => 26.4219
"R_EX_gln__L_e" => -0.0126492
"R_EX_nh4_e" => -5.67116
"R_MALt2_2" => 0.00493596
⋮ => ⋮The function also has the expectable second variant that returns a vector of solutions, in parsimonious_flux_balance_analysis_vec. Here, we utilize it to show how to use different optimizers for finding the optimum and for solving the quadratic problem. That may be preferable if the optimizer qualities differ for the differing tasks. pFBA allows you to specify qp_modifications that are applied after the original optimum is found, and before the quadratic part of the problem solving begins.
vec_soln = parsimonious_flux_balance_analysis_vec(
model,
Tulip.Optimizer; # start with Tulip
modifications = [
change_constraint("R_EX_glc__D_e"; lb = -12, ub = -12),
change_optimizer_attribute("IPM_IterationsLimit", 500), # we may change Tulip-specific attributes here
],
qp_modifications = [
change_optimizer(OSQP.Optimizer), # now switch to OSQP (Tulip wouldn't be able to finish the computation)
change_optimizer_attribute("polish", true), # get an accurate solution, see OSQP's documentation
silence, # and make it quiet.
],
)95-element Vector{Float64}:
-0.006231402533094581
6.847122967587614
8.914435875871503
-5.738183656745528
0.25363156710733
-3.8887743417980833
2.977310372502161
25.64208341306474
6.1991009754149715
1.877264043236841
⋮
-0.0002197051126034427
34.49592944330335
-0.0020752886988600542
-0.0014018715428246171
-0.002009900868927077
-25.642083259145206
0.0161569956400718
0.004729690379188246
3.8887743479784547Minimizing metabolic adjustment analysis (MOMA)
MOMA is a technique used to find a flux distribution that is closest to some reference distribution with respect to the Euclidian norm.
reference_fluxes = parsimonious_flux_balance_analysis_dict( # reference distribution
model,
OSQP.Optimizer;
modifications = [silence, change_optimizer_attribute("polish", true)],
)
moma = minimize_metabolic_adjustment_analysis_dict(
model,
reference_fluxes,
OSQP.Optimizer;
modifications = [
silence,
change_optimizer_attribute("polish", true),
change_constraint("R_CYTBD"; lb = 0.0, ub = 0.0), # find flux distribution closest to the CYTBD knockout
],
)Dict{String, Float64} with 95 entries:
"R_EX_fum_e" => -0.000110206
"R_ACONTb" => 1.02916
"R_TPI" => 8.96839
"R_SUCOAS" => -0.000880332
"R_GLNS" => 0.0157708
"R_EX_pi_e" => -0.228594
"R_PPC" => 7.04745
"R_O2t" => -0.000135095
"R_G6PDH2r" => 0.532365
"R_TALA" => 0.166338
"R_PPCK" => -0.000128188
"R_EX_lac__D_e" => 0.000160938
"R_PGL" => 0.532365
"R_H2Ot" => -5.77982
"R_GLNabc" => 2.40923e-5
"R_EX_co2_e" => 5.23659
"R_EX_gln__L_e" => -2.40628e-5
"R_EX_nh4_e" => -1.29924
"R_MALt2_2" => 4.47106e-5
⋮ => ⋮Composing (more complicated) modifications
COBREXA.jl contains a number of modifications that allow the user to analyze non-standard variants of the classic FBA problem. These include thermodynamic (add_loopless_constraints), capacity (add_crowding_constraints), and kinetic/capacity (add_moment_constraints) modifications. The documentation of each modification details what their purpose is. Here, we will demonstrate how these modifications can be composed to generate even more interesting analyses.
Download the json formatted model so that the reaction identifiers correspond to the Escher map.
!isfile("e_coli_core.json") &&
download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json")
model = load_model("e_coli_core.json")Metabolic model of type JSONModel sparse([9, 51, 55, 64, 65, 34, 44, 59, 66, 64 … 20, 22, 23, 25, 16, 17, 34, 44, 57, 59], [1, 1, 1, 1, 1, 2, 2, 2, 2, 3 … 93, 93, 94, 94, 95, 95, 95, 95, 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
First find a flux distribution that is thermodyamically loopless and incorporates enzyme capacity constraints by composing loopless FBA and FBAwMC.
rid_crowding_weight = Dict(# crowding needs a weight for each flux
rid => 0.004 for rid in reactions(model) if
!looks_like_biomass_reaction(rid) && !looks_like_exchange_reaction(rid)
)
loopless_crowding_fluxes = flux_balance_analysis_dict(
model,
GLPK.Optimizer;
modifications = [
add_crowding_constraints(rid_crowding_weight),
add_loopless_constraints(),
],
)
flux_summary(loopless_crowding_fluxes)Biomass BIOMASS_Ecoli_core_w_GAM: 0.5131 Import EX_o2_e: -11.0901 EX_glc__D_e: -7.4944 EX_nh4_e: -2.7976 EX_pi_e: -1.8874 Export EX_ac_e: 5.7251 EX_co2_e: 11.6832 EX_h2o_e: 15.4205 EX_h_e: 16.017
Next find a flux distribution that satisfies kinetic/capacity constraints using the moment algorithm that is closest (using MOMA) to the loopless crowding solution.
ksas = Dict(rid => 1000.0 for rid in reactions(model)) # make up specific activities of the enyzmes
protein_mass_fraction = 0.56
moment_moma = minimize_metabolic_adjustment_analysis_dict(
model,
loopless_crowding_fluxes,
OSQP.Optimizer;
modifications = [
silence,
change_optimizer_attribute("polish", true),
change_optimizer_attribute("max-iter", 10_000),
change_constraint("EX_glc__D_e", lb = -1000),
change_constraint("CYTBD"; lb = 0, ub = 0),
add_moment_constraints(ksas, protein_mass_fraction;),
],
)
flux_summary(moment_moma)Biomass BIOMASS_Ecoli_core_w_GAM: 0.031 Import EX_glc__D_e: -6.5805 EX_nh4_e: -0.1544 EX_pi_e: -0.1139 EX_akg_e: -0.0065 EX_glu__L_e: -0.0065 EX_gln__L_e: -0.0039 EX_lac__D_e: -0.0 EX_for_e: -0.0 Export EX_pyr_e: 0.0 EX_fru_e: 0.0 EX_fum_e: 0.0 EX_mal__L_e: 0.0 EX_acald_e: 1.1165 EX_ac_e: 1.5923 EX_co2_e: 2.8226 EX_h2o_e: 3.351 EX_etoh_e: 4.7745 EX_succ_e: 5.1153 EX_h_e: 12.4099
Finally, plot the results to inspect the flux distributions visually
using CairoMakie, Escher, ColorSchemes
!isfile("e_coli_core_map.json") && download( # download escher map
"http://bigg.ucsd.edu/escher_map_json/e_coli_core.Core%20metabolism",
"e_coli_core_map.json",
)
maxflux = maximum(abs.(values(moment_moma)))
minflux = minimum(abs.(values(moment_moma)))
color_interp(x) = begin # Scale color of reaction edges to fluxes (manually binned)
normed_x = (abs(x) - minflux) / (maxflux - minflux)
if 0 <= normed_x < 0.01
ColorSchemes.RdYlBu_4[4]
elseif 0.01 <= normed_x < 0.25
ColorSchemes.RdYlBu_4[3]
elseif 0.25 <= normed_x < 0.5
ColorSchemes.RdYlBu_4[2]
else
ColorSchemes.RdYlBu_4[1]
end
end
rc = Dict(k => color_interp(v) for (k, v) in moment_moma) # map reaction id to reaction edge color
fig = Figure();
ax = Axis(fig[1, 1]);
escherplot!(ax, "e_coli_core_map.json", reaction_edge_colors = rc)
hidexdecorations!(ax)
hideydecorations!(ax)
fig
Making custom problem modification functions is really simple due to the tight intergration between COBREXA and JuMP. Look at the source code for the implemented modifications in src\analysis\modifications to get a flavour for it.
This page was generated using Literate.jl.