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.007249566561101
   7.477381918982209
  -5.06437536070867
   0.2234617471361091
  -3.2148950303904567
   2.504309432584187
  21.799492758464645
   4.959985078649209
   1.496983802794076
   ⋮
   3.374635951743139e-7
  29.17582720268202
   9.050785528684483e-9
   4.819754815144906e-8
   9.955685365006965e-9
 -21.799492758464645
  -0.0
  -1.4335949884652254e-9
   3.2148950303904567
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.88405e-8
  "R_EX_lac__D_e" => 2.39305e-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"; 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.26578e-10
  "R_EX_lac__D_e" => 4.37595e-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.23975681387692277, "R_TPI" => 9.75898765499219, "R_SUCOAS" => -0.0032176957414243924, "R_GLNS" => 0.06020056413958599, "R_EX_pi_e" => -0.771072348712683, "R_PPC" => 0.6485114617285411, "R_O2t" => 4.7582203196085796e-17, "R_G6PDH2r" => 0.09755602385269292, "R_TALA" => -0.004979583784785076…), "R_ACONTb" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952399026474, "R_TPI" => 9.75131501024247, "R_SUCOAS" => -5.747880951869093e-11, "R_GLNS" => 0.06218133498161174, "R_EX_pi_e" => -0.7708580440964772, "R_PPC" => 0.6578792999413844, "R_O2t" => 9.41141357513098e-15, "R_G6PDH2r" => 0.12074779842226277, "R_TALA" => 0.0027614296522254166…), "R_TPI" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952434522427, "R_TPI" => 9.48580410486164, "R_SUCOAS" => -8.339944924655711e-11, "R_GLNS" => 0.05358099493381805, "R_EX_pi_e" => -0.770858044539316, "R_PPC" => 0.6004759354094737, "R_O2t" => 5.55270433858668e-15, "R_G6PDH2r" => 0.9172805137167619, "R_TALA" => 0.26827233468830797…), "R_SUCOAS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.29557047715126533, "R_TPI" => 9.791564276282001, "R_SUCOAS" => -0.06949095362007719, "R_GLNS" => 0.053580993819859306, "R_EX_pi_e" => -0.7708580432294316, "R_PPC" => 0.6699668853153778, "R_O2t" => 3.2015419802205517e-15, "R_G6PDH2r" => 9.998532860360709e-10, "R_TALA" => -0.03748783611371247…), "R_GLNS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2406424165846551, "R_TPI" => 9.757081445418361, "R_SUCOAS" => -0.003444294983139651, "R_GLNS" => 0.05358099373941311, "R_EX_pi_e" => -0.770858043243698, "R_PPC" => 0.6515935079655129, "R_O2t" => 2.1428826754852485e-16, "R_G6PDH2r" => 0.10344849359136386, "R_TALA" => -0.0030050052495367886…), "R_EX_pi_e" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22836315673274749, "R_TPI" => 9.789458863462313, "R_SUCOAS" => -7.227521040033789e-11, "R_GLNS" => 0.05412221657627669, "R_EX_pi_e" => -0.7786444930369188, "R_PPC" => 0.6065413498231971, "R_O2t" => 4.510997093485763e-15, "R_G6PDH2r" => 1.4292899902836635e-9, "R_TALA" => -0.037866501224055854…), "R_PPC" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.24854818019576025, "R_TPI" => 9.758144222947491, "R_SUCOAS" => -1.7480083773118704e-13, "R_GLNS" => 0.06071141982240615, "R_EX_pi_e" => -0.770858043206266, "R_PPC" => 0.6004759313481098, "R_O2t" => 1.0686963483993124e-16, "R_G6PDH2r" => 0.10026016103430606, "R_TALA" => -0.004067782766739063…), "R_O2t" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2397547676108426, "R_TPI" => 9.75898409516846, "R_SUCOAS" => -0.003217059911234019, "R_GLNS" => 0.06020007282573227, "R_EX_pi_e" => -0.771072334182146, "R_PPC" => 0.6485122957502956, "R_O2t" => 4.900317536827516e-17, "R_G6PDH2r" => 0.09756671511080944, "R_TALA" => -0.004976019325439087…), "R_G6PDH2r" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.24130975619463457, "R_TPI" => 9.791500043154533, "R_SUCOAS" => -0.0035854003725075737, "R_GLNS" => 0.06121262406038752, "R_EX_pi_e" => -0.7710955979044609, "R_PPC" => 0.654531917733157, "R_O2t" => 2.1547196134586131e-16, "R_G6PDH2r" => 9.284037711554736e-12, "R_TALA" => -0.037499389037928946…), "R_TALA" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.228363187444137, "R_TPI" => 9.789458867641173, "R_SUCOAS" => -8.181399868570223e-9, "R_GLNS" => 0.05412221663542802, "R_EX_pi_e" => -0.7786444672271328, "R_PPC" => 0.6065415024097051, "R_O2t" => 4.108609526215684e-14, "R_G6PDH2r" => 9.826772224935852e-9, "R_TALA" => -0.0378664971692575…)…), Dict("R_EX_fum_e" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.23975681387692277, "R_TPI" => 9.75898765499219, "R_SUCOAS" => -0.0032176957414243924, "R_GLNS" => 0.06020056413958599, "R_EX_pi_e" => -0.771072348712683, "R_PPC" => 0.6485114617285411, "R_O2t" => 4.7582203196085796e-17, "R_G6PDH2r" => 0.09755602385269292, "R_TALA" => -0.004979583784785076…), "R_ACONTb" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.29557047722032287, "R_TPI" => 9.791564276452823, "R_SUCOAS" => -0.0694909508870376, "R_GLNS" => 0.05358099380245838, "R_EX_pi_e" => -0.7708580432013972, "R_PPC" => 0.6699668850561054, "R_O2t" => 5.485806245522459e-15, "R_G6PDH2r" => 5.199605170346847e-10, "R_TALA" => -0.0374878362714516…), "R_TPI" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2417737656014464, "R_TPI" => 9.791564276596779, "R_SUCOAS" => -0.0037135787747932894, "R_GLNS" => 0.061441428632050175, "R_EX_pi_e" => -0.7708580432384502, "R_PPC" => 0.655928635430382, "R_O2t" => 3.5069656080560067e-16, "R_G6PDH2r" => 6.017234894216759e-11, "R_TALA" => -0.037487836426361006…), "R_SUCOAS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.236944501732374, "R_TPI" => 9.757854141412192, "R_SUCOAS" => -4.2138091418429306e-12, "R_GLNS" => 0.060469197785589186, "R_EX_pi_e" => -0.771079467731071, "R_PPC" => 0.6491018624950026, "R_O2t" => 2.0204132234977804e-16, "R_G6PDH2r" => 0.10095078976383551, "R_TALA" => -0.003848341354928834…), "R_GLNS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952365569245, "R_TPI" => 9.791564276059306, "R_SUCOAS" => -5.467877097452668e-11, "R_GLNS" => 0.24468111476804333, "R_EX_pi_e" => -0.7708580432075124, "R_PPC" => 0.6004759321859618, "R_O2t" => 1.2605588558741944e-14, "R_G6PDH2r" => 1.692735141955006e-9, "R_TALA" => -0.037487835880975616…), "R_EX_pi_e" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.24017982275905453, "R_TPI" => 9.758339921865147, "R_SUCOAS" => -0.003337568292850414, "R_GLNS" => 0.060440100765891444, "R_EX_pi_e" => -0.7708580432134876, "R_PPC" => 0.6496831146948908, "R_O2t" => 4.733707741466131e-16, "R_G6PDH2r" => 0.09967306427541843, "R_TALA" => -0.004263481686722996…), "R_PPC" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952347652188, "R_TPI" => 9.791564276534835, "R_SUCOAS" => -1.2664701600525253e-11, "R_GLNS" => 0.05358099376102698, "R_EX_pi_e" => -0.7708580432219333, "R_PPC" => 0.9826761788062314, "R_O2t" => -7.453482237529607e-16, "R_G6PDH2r" => 2.4910685402287665e-10, "R_TALA" => -0.037487836363391924…), "R_O2t" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2397588586943713, "R_TPI" => 9.758991246140967, "R_SUCOAS" => -0.0032183313812690247, "R_GLNS" => 0.060201086707893314, "R_EX_pi_e" => -0.771072363038656, "R_PPC" => 0.6485106398048069, "R_O2t" => 4.625872330505629e-17, "R_G6PDH2r" => 0.09754523878537202, "R_TALA" => -0.004983179503915471…), "R_G6PDH2r" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952341994342, "R_TPI" => 9.485804077752663, "R_SUCOAS" => -6.54325258461258e-14, "R_GLNS" => 0.05358099373347381, "R_EX_pi_e" => -0.7708580431954545, "R_PPC" => 0.6004759313306418, "R_O2t" => 2.4905012766984472e-17, "R_G6PDH2r" => 0.9172805966273293, "R_TALA" => 0.26827236243144115…), "R_TALA" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952341983883, "R_TPI" => 9.485804077753238, "R_SUCOAS" => -6.73908436300807e-14, "R_GLNS" => 0.05358099373340432, "R_EX_pi_e" => -0.7708580431954012, "R_PPC" => 0.6004759313306722, "R_O2t" => 2.2990512232112453e-18, "R_G6PDH2r" => 0.9172805966256342, "R_TALA" => 0.2682723624308808…)…))
fva_maxs["R_EX_ac_e"]["R_EX_ac_e"] # get the maximal acetate exchange flux
8.518549434877636

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.29497e-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.4903e-18, 0.112692)  (-5.53655e-18, 0.112692)
 (4.67779e-13, 0.112782)  (7.93011, 0.105831)
 (-0.0378665, 0.211663)   (2.62444, 0.105831)
 ⋮                        
 (3.0129e-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.0506e-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.59454e-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.4903e-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),
    ],
)
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)
        silence, # and make it quiet.
    ],
)
95-element Vector{Float64}:
  -0.006231402534474541
   6.847122967595623
   8.914435875868813
  -5.738183656748431
   0.2536315671085878
  -3.8887743417982255
   2.977310372500073
  25.642083413072584
   6.199100975412263
   1.877264043235936
   ⋮
  -0.00021970511206870563
  34.495929443314076
  -0.0020752886986875763
  -0.0014018715430539478
  -0.0020099008686134795
 -25.642083259153033
   0.016156995638076037
   0.00472969038077297
   3.888774347978596

Flux balance analysis with molecular crowding (FBAwMC)

Flux balance with molecular crowding incorporates enzyme capacity constraints into the classic flux balance analysis algorithm. Essentially, an extra constraint is added to the optimization problem: ∑ wᵢ × vᵢ ≤ 1 where wᵢ weights each internal flux vᵢ. See Beg, Qasim K., et al. "Intracellular crowding defines the mode and sequence of substrate uptake by Escherichia coli and constrains its metabolic activity." Proceedings of the National Academy of Sciences 104.31 (2007): 12663-12668. for more details.

First load the model

model = load_model("e_coli_core.json")
Metabolic model of type JSONModel

⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢄⠀⠀⠀⠀⠀⠀⠀⠈⠶⠴⡆⠀⠀⠀⠀⠀⠀
⡀⢐⣀⢀⡀⡒⢒⣐⠀⣂⣂⠀⣂⣂⢂⠀⢀⠀⠀⠀⠀⠀⢀⠄⠀⠀⠀⢂⠀⢂⣀⣐⡒⡀⠆⢙⣀⠀⡀⠀
⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠰⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠈⢑⣀⣀⠀⠀
⠀⠀⠃⠀⠃⠀⠀⠀⠘⠀⡇⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠀⠀⡜⠀⡄⣤⢠⠘⠙⢣⡇⠘
⠀⠐⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠐⠁⠉⠀⠀⠀⠀⠀⠘⠄
⠀⢐⠀⠂⠀⠄⠠⠠⠀⠠⠆⠀⠄⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠠⠀⠠⠀⠀⢀⠀⠀⠠⠀⠀⠁
⢀⠐⠀⠨⢀⠁⠈⣈⠀⢁⣁⠀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠁⢀⠀⢊⠉⠀⠀⠀⢀⠀⣀⠀⢀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡈⠀⡀⠆⠀⠆⠀⡀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠆⠀
⠀⠀⠂⠀⡂⠀⠀⠁⠀⠀⠀⠈⠁⠀⠀⠀⠄⠄⢁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀
⠈⠀⠁⠀⠀⢀⡀⠀⠠⠁⠁⠀⠑⠀⠐⠲⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠂⠀⠀⠀⠀⠀⠀⠊⠀⠀⠀⠈
⠄⠠⢠⠀⠰⠀⠠⠀⠤⠦⠄⠈⠀⠀⠀⠠⠀⠁⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠤⠄⠄⠠⠀⠀⠀⠀⠀
⠂⠐⠀⠀⠐⡠⢐⠘⢃⠒⠂⡀⠄⠀⠀⠐⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠒⠀⢀⢀⠀⠀⣀⠀⢀
⠈⠀⠁⠀⡀⠀⠀⠀⠈⠁⠅⠀⠁⠀⢀⠈⠄⠔⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠀⠈
⠣⠁⠀⠀⠀⠀⠀⠀⠀⠀⠁⠀⠀⠀⠈⠀⠁⠁⠀⠈⡀⠀⠀⠀⠀⠀⠐⢣⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⡄⠀⠀⠀⠀⠂⠄⠤⠀⠀⠈⠂⠀⠀⠀⠀⠠⠀⠊⠒⣠⠀⠀⠀⠀⠀⠀⠀⠀⠀
Number of reactions: 95
Number of metabolites: 72

Next, simulate the model over a range of substrate uptake rates.

without_crowding = Dict{Float64,Vector{Float64}}()
with_crowding = Dict{Float64,Vector{Float64}}()
glucose_uptakes = collect(-(1.0:0.5:20))

for glc in glucose_uptakes
    no_crowding = flux_balance_analysis_dict( # classic FBA
        model,
        Tulip.Optimizer;
        modifications = [
            change_optimizer_attribute("IPM_IterationsLimit", 1000),
            change_constraint("EX_glc__D_e"; lb = glc),
        ],
    )

    without_crowding[glc] =
        [no_crowding["BIOMASS_Ecoli_core_w_GAM"], no_crowding["EX_ac_e"]]

    crowding = flux_balance_analysis_dict( # FBAwMC
        model,
        Tulip.Optimizer;
        modifications = [
            change_optimizer_attribute("IPM_IterationsLimit", 1000),
            add_crowding_constraint(0.004), # crowding constraint gets added here
            change_constraint("EX_glc__D_e"; lb = glc),
        ],
    )

    with_crowding[glc] = [crowding["BIOMASS_Ecoli_core_w_GAM"], crowding["EX_ac_e"]]
end

Finally, plot the results to compare classic FBA with FBAwMC.

using CairoMakie
fig = Figure();
ax1 = Axis(fig[1, 1]);
lines!(
    ax1,
    -glucose_uptakes,
    [without_crowding[glc][1] for glc in glucose_uptakes],
    label = "no crowding",
    linewidth = 5,
    linestyle = :dash,
)
lines!(
    ax1,
    -glucose_uptakes,
    [with_crowding[glc][1] for glc in glucose_uptakes],
    label = "with crowding",
    linewidth = 5,
    linestyle = :dot,
)
ax1.ylabel = "Growth rate [1/h]"
ax2 = Axis(fig[2, 1])
lines!(
    ax2,
    -glucose_uptakes,
    [without_crowding[glc][2] for glc in glucose_uptakes],
    label = "no crowding",
    linewidth = 5,
    linestyle = :dash,
)
lines!(
    ax2,
    -glucose_uptakes,
    [with_crowding[glc][2] for glc in glucose_uptakes],
    label = "with crowding",
    linewidth = 5,
    linestyle = :dot,
)
fig[1:2, 2] = Legend(fig, ax1, "Constraint")
ax2.xlabel = "Glucose uptake rate [mmol/gDW/h]"
ax2.ylabel = "Acetate production rate\n[mmol/gDW/h]"
fig

Notice that overflow metabolism is modeled by incorporating the crowding constraint.

Tip: code your own modification like [`add_crowding`](@ref)

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.