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
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.
import Pkg
Pkg.add("Tulip")
Pkg.add("OSQP")
using Tulip
using 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" => 8.387567703061316, "R_TPI" => 7.532255418870419, "R_SUCOAS" => -2.1456715700928, "R_GLNS" => 11.72801152935013, "R_EX_pi_e" => -0.781841023618242, "R_PPC" => 12.709428540552382, "R_O2t" => 30.207306741855934, "R_G6PDH2r" => 6.3963713000040885, "R_TALA" => 2.094101813498349…), "R_ACONTb" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2260795234341539, "R_TPI" => 0.27629819544371975, "R_SUCOAS" => -2.8375942206305713e-12, "R_GLNS" => 9.663346629147215, "R_EX_pi_e" => -0.7708580431967099, "R_PPC" => 13.616949989618535, "R_O2t" => 30.324546231187313, "R_G6PDH2r" => 28.16720077460885, "R_TALA" => 9.351579088425249…), "R_TPI" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952365224862, "R_TPI" => -7.023645995325829, "R_SUCOAS" => -5.989684359376451e-11, "R_GLNS" => 14.308313357789281, "R_EX_pi_e" => -0.7708580432585497, "R_PPC" => 18.802685220533647, "R_O2t" => 50.4222650020313, "R_G6PDH2r" => 50.44563081566434, "R_TALA" => 16.777722435440715…), "R_SUCOAS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 17.041289795789567, "R_TPI" => 9.791564276526747, "R_SUCOAS" => -16.815210272246723, "R_GLNS" => 19.47631848364111, "R_EX_pi_e" => -0.770858043196644, "R_PPC" => 25.4355271739484, "R_O2t" => 48.49828617064436, "R_G6PDH2r" => 2.4419522059618934e-10, "R_TALA" => -0.03748783636298141…), "R_GLNS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 8.061973806526826, "R_TPI" => 7.523620762504264, "R_SUCOAS" => -2.0664872403347907, "R_GLNS" => 0.0535809937339624, "R_EX_pi_e" => -0.7708580431966973, "R_PPC" => 14.947771300763096, "R_O2t" => 29.036876963512825, "R_G6PDH2r" => 6.399244385839639, "R_TALA" => 2.0955936255021808…), "R_EX_pi_e" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 8.388812649215613, "R_TPI" => 7.533309613876824, "R_SUCOAS" => -2.1466235942183447, "R_GLNS" => 11.688210901810553, "R_EX_pi_e" => -0.7865095839138668, "R_PPC" => 12.645801274331156, "R_O2t" => 30.188166942395195, "R_G6PDH2r" => 6.388553068133392, "R_TALA" => 2.09126869800525…), "R_PPC" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 9.14053048210921, "R_TPI" => 7.815024382280643, "R_SUCOAS" => -1.4014273460315068, "R_GLNS" => 15.195578469146987, "R_EX_pi_e" => -0.7818087293451119, "R_PPC" => 8.00091769193397e-10, "R_O2t" => 31.204423146044512, "R_G6PDH2r" => 5.543458812160874, "R_TALA" => 1.809799221396707…), "R_O2t" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.23347557615742534, "R_TPI" => 9.783623631549453, "R_SUCOAS" => -0.001629549873797233, "R_GLNS" => 0.06270823060221115, "R_EX_pi_e" => -0.7716851675773635, "R_PPC" => 0.6737424037490961, "R_O2t" => 2.0683004221209886e-10, "R_G6PDH2r" => 0.01992471198008615, "R_TALA" => -0.03088648992748138…), "R_G6PDH2r" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 10.251784177340285, "R_TPI" => 9.659615645294526, "R_SUCOAS" => -3.404492835948122, "R_GLNS" => 11.754546921095894, "R_EX_pi_e" => -0.7818183172437788, "R_PPC" => 12.84942993145147, "R_O2t" => 29.407516338347314, "R_G6PDH2r" => 1.2044971310956056e-13, "R_TALA" => -0.038020848928905275…), "R_TALA" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 10.35557912119309, "R_TPI" => 9.65911239712571, "R_SUCOAS" => -4.306503172360712, "R_GLNS" => 11.818606548613598, "R_EX_pi_e" => -0.786509583913964, "R_PPC" => 13.469379343614088, "R_O2t" => 29.691047792422857, "R_G6PDH2r" => 4.749673144081168e-13, "R_TALA" => -0.03824899137239327…)…), Dict("R_EX_fum_e" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 8.387567703061316, "R_TPI" => 7.532255418870419, "R_SUCOAS" => -2.1456715700928, "R_GLNS" => 11.72801152935013, "R_EX_pi_e" => -0.781841023618242, "R_PPC" => 12.709428540552382, "R_O2t" => 30.207306741855934, "R_G6PDH2r" => 6.3963713000040885, "R_TALA" => 2.094101813498349…), "R_ACONTb" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 17.041289795007938, "R_TPI" => 9.791564276433519, "R_SUCOAS" => -5.688082121108331, "R_GLNS" => 19.780355774841556, "R_EX_pi_e" => -0.7708580433275182, "R_PPC" => 16.936579264014647, "R_O2t" => 48.30857236846392, "R_G6PDH2r" => 3.833409490124206e-10, "R_TALA" => -0.03748783632297794…), "R_TPI" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 10.554697110858058, "R_TPI" => 9.79156427517227, "R_SUCOAS" => -4.36079829072734, "R_GLNS" => 12.195418956692743, "R_EX_pi_e" => -0.7708580443890825, "R_PPC" => 13.848078463930056, "R_O2t" => 30.26371446753615, "R_G6PDH2r" => 1.996871968753506e-9, "R_TALA" => -0.03748783583683612…), "R_SUCOAS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 8.057908581115928, "R_TPI" => 7.332654280738366, "R_SUCOAS" => -1.2187147897617456e-10, "R_GLNS" => 11.572637789614943, "R_EX_pi_e" => -0.781840987603548, "R_PPC" => 11.178325328275292, "R_O2t" => 29.841134825372734, "R_G6PDH2r" => 6.989508449861914, "R_TALA" => 2.2918141985356772…), "R_GLNS" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 17.04128979542336, "R_TPI" => 9.791564276164623, "R_SUCOAS" => -16.769668728120884, "R_GLNS" => 127.20986479797551, "R_EX_pi_e" => -0.7708580431963201, "R_PPC" => 0.55493438797057, "R_O2t" => 50.84039512729372, "R_G6PDH2r" => 1.3335673812250126e-9, "R_TALA" => -0.037487835999915696…), "R_EX_pi_e" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 8.477625183335908, "R_TPI" => 7.508226727537374, "R_SUCOAS" => -2.1740673997402995, "R_GLNS" => 11.955467607327458, "R_EX_pi_e" => -0.7708580431986837, "R_PPC" => 12.809131647158342, "R_O2t" => 30.488106628893533, "R_G6PDH2r" => 6.466262473804337, "R_TALA" => 2.117932988156977…), "R_PPC" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 17.041289790111023, "R_TPI" => 9.791564273069904, "R_SUCOAS" => -16.815210261955922, "R_GLNS" => 0.05358099899733933, "R_EX_pi_e" => -0.7708580432275103, "R_PPC" => 127.7567596864386, "R_O2t" => 50.84039511913468, "R_G6PDH2r" => 1.0228880908192078e-8, "R_TALA" => -0.03748783303680075…), "R_O2t" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 13.488565849967824, "R_TPI" => 6.238840330806542, "R_SUCOAS" => -3.772617647963326, "R_GLNS" => 20.283358516037556, "R_EX_pi_e" => -0.7708580432002806, "R_PPC" => 17.883676633529095, "R_O2t" => 50.84039512688941, "R_G6PDH2r" => 10.658171837390052, "R_TALA" => 3.515236109352078…), "R_G6PDH2r" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.2260795280734479, "R_TPI" => -7.023645980797208, "R_SUCOAS" => -1.2263907749526956e-9, "R_GLNS" => 14.973603529019034, "R_EX_pi_e" => -0.7708580432619148, "R_PPC" => 18.825207611043435, "R_O2t" => 50.530923759776464, "R_G6PDH2r" => 50.445630770442904, "R_TALA" => 16.7777224203659…), "R_TALA" => Dict("R_EX_fum_e" => -0.0, "R_ACONTb" => 0.22607952440356835, "R_TPI" => -7.023645988176, "R_SUCOAS" => -2.8144007763514797e-10, "R_GLNS" => 14.338969784762646, "R_EX_pi_e" => -0.7708580432375506, "R_PPC" => 18.868338300612862, "R_O2t" => 50.4236100549293, "R_G6PDH2r" => 50.445630793530995, "R_TALA" => 16.777722428063736…)…))
fva_maxs["R_EX_ac_e"]["R_EX_ac_e"] # get the maximal acetate exchange flux
16.815210271774536
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.