Finding balance and variability of constraint-based models

Here we use flux_balance_analysis, flux_variability_analysis, and parsimonious_flux_balance_analysis of COBREXA.jl functions 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 COBREXA
Tip: use `?` to get quick help about functions

When 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

⠀⠈⢀⠀⡀⠀⠀⠀⠀⡠⠂⠀⠀⠀⠀⠈⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠀⠀⠀⠀⢀⠐⡀⠀⠀⠀⠀⠄
⠀⠐⠀⠀⠀⠀⠀⠀⡠⠂⠀⠀⠀⠀⢰⠱⣀⠀⡄⢐⠀⠀⢀⠀⠀⠀⡂⠄⠔⠁⠰⠀⠠⠀⣆⠀⠄⢠⢀⠄
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠠⠀⠀⠐⠀⠀⠀⠀⠀⠀⢀⠀⠀⠐⠀⠂⠀⠀⠀⠄⠀⠐⠀⢁⠄⠀⠀⠀⠀⠀
⠀⢀⠀⠐⡈⠀⡀⠀⠂⠀⣀⠀⠑⡈⢀⠀⠀⠀⠀⠀⡀⡠⠀⡀⠰⠁⠈⠂⠁⠀⠠⠀⠀⠂⡂⠀⠂⠂⠀⠀
⠠⠀⠐⠀⠂⠀⠀⢀⠀⠀⠀⠀⠊⠀⡐⠊⠐⠀⠀⠀⠀⠀⠐⠀⠂⠀⠀⠐⠀⠀⠀⠀⠀⠁⠃⠠⠀⠁⠐⠀
⠀⠠⠀⡀⠄⠀⠀⠂⠀⠀⠀⠠⠀⠠⠀⠀⠄⠀⠨⠀⠀⠀⠐⠀⠀⠄⢀⠀⠀⠀⠈⠀⠀⠀⠁⠄⠀⠀⠀⠀
⠀⢐⠐⠀⠄⠀⡂⠀⢐⠀⠀⠀⠀⠂⢀⢀⠐⠂⡀⠈⠀⠀⠀⠂⠀⠈⠀⡀⡐⠀⢄⠀⢀⠀⡆⠀⡀⣀⡀⡐
⠀⠈⠀⠀⠀⠀⠀⠐⢂⠀⢀⠀⠈⠀⠀⠀⠀⠀⠠⠀⠀⠠⠀⠀⠀⠈⠂⠀⠀⠀⠄⠐⠐⠀⠁⠀⠀⠑⠁⠀
⠂⠠⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠠⠈⠀⠀⠀⠀⠀⠁⠀⠀⠠⠐⠀⠁⠈⠀⠁⢀⠀⠀⠀⠀⠀⠀⠀⠀⠌⠀
⠀⠀⠂⢨⠀⡀⠀⠐⠁⠐⠀⠐⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠢⠒⠈⠐⠐⠁⠂⠀⠀⠀⠄⠓⠕⠂⠃⠁⠀⠐
⠠⠀⠨⠀⠁⠤⠄⠀⠁⡄⠀⠂⠠⠄⢈⠌⠠⠄⠀⢀⠀⠀⠀⠄⠨⠀⡤⠀⢀⠀⢀⠠⠀⠁⡔⠨⠀⠈⠄⠀
⠀⢀⢀⣀⠀⡠⡒⢀⢀⣀⠀⢀⣀⡀⢀⠀⢀⠀⡀⠀⡀⠀⠈⣀⠀⢀⣀⠀⡀⠀⢀⠁⢀⣀⣀⡀⡠⡀⡀⣀
⠀⠄⠀⠀⠀⠀⠀⠀⠀⠂⠁⠀⠀⠀⠀⠀⠀⣠⠀⠀⡀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀
⢀⠂⠀⠀⠂⠀⠈⠀⠐⠀⠀⠀⠁⠀⠀⠀⡀⠔⠑⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠢⠀⠀⡀⠂⠈⠀⠀⠀⠄
⠀⠐⠀⠀⡂⠀⠂⠀⠀⠀⠒⠐⠄⠂⠐⠀⠘⡀⠀⠠⡂⠃⠀⠂⠄⠂⠀⠀⠀⠀⡀⠀⡀⠀⡂⠂⠀⠀⢀⠀
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.

Note: OSQP can be sensitive

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

Flux 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.214895030387032
dict_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
  ⋮               => ⋮

Modifications

Often it is desirable to add a slight modififaction 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.

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", -12, -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
  ⋮               => ⋮

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", -10, -10),
        change_constraint("R_EX_o2_e", 0.0, 0.0),
    ],
)
(Dict("R_EX_fum_e" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.23975681246428646, "R_TPI" => 9.758987664410151, "R_SUCOAS" => -0.0032176953711697304, "R_GLNS" => 0.060200571606221165, "R_EX_pi_e" => -0.7710723486676665, "R_PPC" => 0.648511462789616, "R_O2t" => 4.7620907185810674e-17, "R_G6PDH2r" => 0.09755599563531388, "R_TALA" => -0.004979593188388901…), "R_ACONTb" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952399025047, "R_TPI" => 9.751315002264594, "R_SUCOAS" => -5.7476116849648756e-11, "R_GLNS" => 0.062181345572819555, "R_EX_pi_e" => -0.7708580440964182, "R_PPC" => 0.6578793157880826, "R_O2t" => 9.410658208197004e-15, "R_G6PDH2r" => 0.1207478223560313, "R_TALA" => 0.002761437630187624…), "R_TPI" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952434490422, "R_TPI" => 9.485804104911805, "R_SUCOAS" => -8.335641656995056e-11, "R_GLNS" => 0.053580994933664666, "R_EX_pi_e" => -0.7708580445468374, "R_PPC" => 0.6004759354089005, "R_O2t" => 5.507347171694764e-15, "R_G6PDH2r" => 0.9172805135279788, "R_TALA" => 0.2682723346211751…), "R_SUCOAS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2955704771566707, "R_TPI" => 9.791564276285534, "R_SUCOAS" => -0.06949095362639839, "R_GLNS" => 0.05358099381935753, "R_EX_pi_e" => -0.7708580432289044, "R_PPC" => 0.6699668853183516, "R_O2t" => 3.13394692633907e-15, "R_G6PDH2r" => 9.897941427456658e-10, "R_TALA" => -0.03748783611687239…), "R_GLNS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2406424227465452, "R_TPI" => 9.757081490883413, "R_SUCOAS" => -0.0034442967413112913, "R_GLNS" => 0.05358099373941643, "R_EX_pi_e" => -0.7708580432437482, "R_PPC" => 0.6515935275838877, "R_O2t" => 2.1426954824772113e-16, "R_G6PDH2r" => 0.10344835719616492, "R_TALA" => -0.0030050507146060132…), "R_EX_pi_e" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22836315673274257, "R_TPI" => 9.78945886346248, "R_SUCOAS" => -7.227000020845847e-11, "R_GLNS" => 0.05412221657627899, "R_EX_pi_e" => -0.7786444930367454, "R_PPC" => 0.6065413498231926, "R_O2t" => 4.51046301140303e-15, "R_G6PDH2r" => 1.4292033457498604e-9, "R_TALA" => -0.03786650122403871…), "R_PPC" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.24854826185014645, "R_TPI" => 9.758144304106896, "R_SUCOAS" => -1.7527365950313536e-13, "R_GLNS" => 0.060711429288066124, "R_EX_pi_e" => -0.7708580432063145, "R_PPC" => 0.6004759313481499, "R_O2t" => 1.0676969754863362e-16, "R_G6PDH2r" => 0.10025991755604913, "R_TALA" => -0.004067863926160701…), "R_O2t" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2397547670699771, "R_TPI" => 9.758984098152247, "R_SUCOAS" => -0.003217059763568135, "R_GLNS" => 0.06020007476474629, "R_EX_pi_e" => -0.7710723341748742, "R_PPC" => 0.6485122957709831, "R_O2t" => 4.9009916436786396e-17, "R_G6PDH2r" => 0.09756670616534588, "R_TALA" => -0.004976022306907555…), "R_G6PDH2r" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2413098384035294, "R_TPI" => 9.791500043356482, "R_SUCOAS" => -0.003585425307283908, "R_GLNS" => 0.06121279515312429, "R_EX_pi_e" => -0.7710955971575917, "R_PPC" => 0.6545319804571721, "R_O2t" => 2.1562429801722262e-16, "R_G6PDH2r" => 9.283567622367387e-12, "R_TALA" => -0.03749938900160984…), "R_TALA" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22836318744798303, "R_TPI" => 9.789458867641589, "R_SUCOAS" => -8.182422356774574e-9, "R_GLNS" => 0.054122216635179424, "R_EX_pi_e" => -0.7786444672236877, "R_PPC" => 0.6065415024301886, "R_O2t" => 4.108918888876329e-14, "R_G6PDH2r" => 9.828550567579093e-9, "R_TALA" => -0.03786649716846829…)…), Dict("R_EX_fum_e" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.23975681246428646, "R_TPI" => 9.758987664410151, "R_SUCOAS" => -0.0032176953711697304, "R_GLNS" => 0.060200571606221165, "R_EX_pi_e" => -0.7710723486676665, "R_PPC" => 0.648511462789616, "R_O2t" => 4.7620907185810674e-17, "R_G6PDH2r" => 0.09755599563531388, "R_TALA" => -0.004979593188388901…), "R_ACONTb" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.29557047721986357, "R_TPI" => 9.791564276452807, "R_SUCOAS" => -0.06949095088677924, "R_GLNS" => 0.053580993802476364, "R_EX_pi_e" => -0.7708580432014581, "R_PPC" => 0.669966885055626, "R_O2t" => 5.4867723002985896e-15, "R_G6PDH2r" => 5.200538446278462e-10, "R_TALA" => -0.0374878362714213…), "R_TPI" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.24177376951630916, "R_TPI" => 9.791564276596803, "R_SUCOAS" => -0.0037135796370060203, "R_GLNS" => 0.06144138338428863, "R_EX_pi_e" => -0.7708580432384102, "R_PPC" => 0.655928636691913, "R_O2t" => 3.503241648282044e-16, "R_G6PDH2r" => 6.012119002464294e-11, "R_TALA" => -0.0374878364263768…), "R_SUCOAS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.23694450009958637, "R_TPI" => 9.7578541338508, "R_SUCOAS" => -4.213751498134161e-12, "R_GLNS" => 0.06046918564757093, "R_EX_pi_e" => -0.7710794677351324, "R_PPC" => 0.6491018708187021, "R_O2t" => 2.0203965226776368e-16, "R_G6PDH2r" => 0.10095081244472187, "R_TALA" => -0.0038483337948307357…), "R_GLNS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952365569617, "R_TPI" => 9.791564276059317, "R_SUCOAS" => -5.4673258263777484e-11, "R_GLNS" => 0.24468111476697532, "R_EX_pi_e" => -0.7708580432075716, "R_PPC" => 0.6004759321859329, "R_O2t" => 1.2604324269135335e-14, "R_G6PDH2r" => 1.6926999040525598e-9, "R_TALA" => -0.037487835880989515…), "R_EX_pi_e" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.24017982085639003, "R_TPI" => 9.758339902976445, "R_SUCOAS" => -0.003337567801679041, "R_GLNS" => 0.06044010139704514, "R_EX_pi_e" => -0.7708580432135336, "R_PPC" => 0.6496831096146779, "R_O2t" => 4.733798443891441e-16, "R_G6PDH2r" => 0.09967312094148192, "R_TALA" => -0.004263462798038352…), "R_PPC" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952347654484, "R_TPI" => 9.791564276535054, "R_SUCOAS" => -1.2665682284481008e-11, "R_GLNS" => 0.05358099376103458, "R_EX_pi_e" => -0.7708580432220014, "R_PPC" => 0.9826761788036762, "R_O2t" => -7.451157181351547e-16, "R_G6PDH2r" => 2.4914584091711995e-10, "R_TALA" => -0.03748783636345861…), "R_O2t" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.23975885844590922, "R_TPI" => 9.758991247095977, "R_SUCOAS" => -0.0032183313066518893, "R_GLNS" => 0.060201087253185896, "R_EX_pi_e" => -0.7710723630386553, "R_PPC" => 0.6485106398373328, "R_O2t" => 4.6259878518155594e-17, "R_G6PDH2r" => 0.09754523592034138, "R_TALA" => -0.0049831804589255705…), "R_G6PDH2r" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952341995713, "R_TPI" => 9.485804077754599, "R_SUCOAS" => -7.70632409774389e-14, "R_GLNS" => 0.053580993733477245, "R_EX_pi_e" => -0.7708580431955058, "R_PPC" => 0.6004759313306891, "R_O2t" => 2.491455505240063e-17, "R_G6PDH2r" => 0.9172805966214608, "R_TALA" => 0.2682723624294833…), "R_TALA" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2260795234203813, "R_TPI" => 9.485804077764492, "R_SUCOAS" => -2.2070646706247915e-14, "R_GLNS" => 0.05358099373347766, "R_EX_pi_e" => -0.7708580431970373, "R_PPC" => 0.6004759313313932, "R_O2t" => 2.2011253936851523e-18, "R_G6PDH2r" => 0.9172805965899173, "R_TALA" => 0.2682723624188306…)…))
fva_maxs["R_EX_ac_e"]["R_EX_ac_e"] # get the maximal acetate exchange flux
8.518549434877563

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 = [
        change_optimizer_attribute("verbose", false), # silence the optimizer (OSQP is very verbose by default)
        change_constraint("R_EX_glc__D_e", -12, -12),
    ],
)
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", -12, -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("verbose", false), # and silence it.
    ],
)
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.8887743479784547

This page was generated using Literate.jl.