Bio Applied Metabolic Flux

Metabolic Flux Analysis and Constraint-Based Modeling

When to Use

• Predicting growth rate and reaction fluxes from a genome-scale metabolic model (GEM) given a nutrient medium.
• Finding essential genes / candidate drug targets via single- or double-gene knockout simulation.
• Testing whether a gene deletion is a synthetic lethal in combination with another gene or drug.
• Constraining a generic GEM (e.g. Recon3D) with RNA-seq/microarray expression data (GIMME, iMAT, E-Flux) to get a context/tissue-specific model.
• Interpreting 13C-labeling flux data (13C-MFA) to resolve parallel pathways (e.g. glycolysis vs. pentose phosphate pathway).

Version Compatibility

• cobrapy ≥ 0.29, Python ≥ 3.9
• optlang ≥ 1.8 (solver interface; default GLPK, or install cplex/gurobipy for large models)
• Model formats: SBML (.xml), JSON, or MATLAB .mat from BiGG (http://bigg.ucsd.edu) or VMH (https://www.vmh.life, Recon3D/AGORA2)

Prerequisites

pip install cobra

Concepts: stoichiometric matrix S (reactions x metabolites), flux vector v, steady-state assumption S·v = 0, reaction bounds (lower/upper define reversibility and capacity), exchange reactions (model nutrient uptake/secretion at the system boundary), objective function (typically biomass maximization).

FBA: Load a Model and Maximize Growth

Goal: predict the maximal growth rate and the flux distribution through central metabolism under a defined medium. Approach: load a GEM, set exchange-reaction bounds to define the medium, call model.optimize(), then inspect nonzero fluxes.

import cobra
from cobra.io import load_model

def run_fba(model_name="e_coli_core", glucose_uptake=-10.0):
    """Load a GEM and maximize biomass under a glucose-limited aerobic medium.

    glucose_uptake is negative (uptake = influx, COBRA convention: exchange
    reaction lower bound is the max uptake rate, in mmol/gDW/h).
    """
    model = load_model(model_name)  # bundled models: "e_coli_core", "iJO1366", "textbook"

    # Set the medium: close all exchange reactions, then open the ones you want
    medium = model.medium
    if "EX_glc__D_e" in medium:
        medium["EX_glc__D_e"] = -glucose_uptake  # cobra medium dict stores uptake as positive value
        model.medium = medium

    solution = model.optimize()
    if solution.status != "optimal":
        raise RuntimeError(f"FBA infeasible/unbounded: status={solution.status}")

    print(f"Growth rate (objective): {solution.objective_value:.4f} 1/h")
    active_fluxes = solution.fluxes[solution.fluxes.abs() > 1e-3].sort_values()
    return solution, active_fluxes

if __name__ == "__main__":
    sol, fluxes = run_fba()
    assert sol.status == "optimal"
    assert sol.objective_value > 0, "expected positive growth on glucose minimal medium"
    print(fluxes.tail(10))

Flux Variability Analysis (FVA)

Goal: find the min/max possible flux through each reaction while keeping growth near-optimal, to reveal alternate optimal pathways. Approach: cobra.flux_analysis.flux_variability_analysis with a fraction_of_optimum threshold.

from cobra.flux_analysis import flux_variability_analysis

def run_fva(model, reaction_ids=None, fraction_of_optimum=0.9):
    """Compute flux ranges for reactions while requiring >=90% of max growth."""
    reaction_ids = reaction_ids or [r.id for r in model.reactions]
    fva_result = flux_variability_analysis(
        model, reaction_list=reaction_ids, fraction_of_optimum=fraction_of_optimum
    )
    return fva_result  # DataFrame indexed by reaction id, columns minimum/maximum

Gene Knockout and Drug Target Prediction

Goal: identify essential genes (candidate drug targets) and synthetic-lethal gene pairs. Approach: single_gene_deletion / double_gene_deletion recompute growth with each gene's reactions constrained to zero via gene-protein-reaction (GPR) rules; an essential gene drops growth near zero.

from cobra.flux_analysis import single_gene_deletion, double_gene_deletion

def find_essential_genes(model, growth_cutoff=0.01):
    """Single-gene deletion scan; genes whose knockout growth < cutoff are essential."""
    wt_growth = model.slim_optimize()
    results = single_gene_deletion(model)  # DataFrame: ids, growth, status
    results = results.reset_index(drop=False)
    results["relative_growth"] = results["growth"] / wt_growth
    essential = results[results["relative_growth"] < growth_cutoff]
    return essential.sort_values("relative_growth")

def synthetic_lethal_pairs(model, gene_ids, growth_cutoff=0.01):
    """Double-gene deletion screen restricted to a candidate gene list (faster than all-pairs)."""
    pairs = double_gene_deletion(model, gene_list1=gene_ids)
    pairs["relative_growth"] = pairs["growth"] / model.slim_optimize()
    return pairs[pairs["relative_growth"] < growth_cutoff]

Compare predicted essential genes against DepMap CRISPR-screen dependency scores for orthogonal validation before proposing drug targets.

Transcriptomics-Constrained FBA (GIMME / iMAT)

Goal: turn a generic GEM into a context-specific model using RNA-seq expression. Approach: map expression to reactions via GPR rules, then constrain fluxes; cobrapy ships GIMME, and cobra.flux_analysis supports iMAT-style integration through create_context_specific_model helpers in COBRA add-ons (troppo/cobamp) for iMAT/E-Flux — GIMME is native.

from cobra.flux_analysis import gimme

def context_specific_model(model, expression_by_gene, cutoff_percentile=25):
    """Build a context-specific model with GIMME: penalize flux through
    reactions whose associated genes are lowly expressed.

    expression_by_gene: dict[gene_id -> float] (e.g. TPM values from RNA-seq).
    """
    import pandas as pd
    expr = pd.Series(expression_by_gene)
    cutoff = expr.quantile(cutoff_percentile / 100)
    # Map gene expression to reaction expression via GPR (min of genes in an AND rule)
    reaction_expression = {}
    for rxn in model.reactions:
        genes = [g.id for g in rxn.genes]
        if genes:
            reaction_expression[rxn.id] = min(expr.get(g, 0.0) for g in genes)
    gimme_result = gimme.gimme(model, reaction_expression, cutoff=cutoff)
    return gimme_result

For iMAT / E-Flux specifically, use the troppo package (pip install troppo), which implements both algorithms against a cobrapy model and expression matrix.

13C Metabolic Flux Analysis

Isotope-labeling experiments (feed U-13C glucose), then measure mass isotopologue distributions (MID) by GC-MS/LC-MS on downstream metabolites. Fit intracellular fluxes that best reproduce the observed MIDs using dedicated 13C-MFA tools (INCA, OpenMebius) — this resolves flux splits that FBA alone cannot (e.g. glycolysis vs. pentose phosphate pathway), because FBA only optimizes an objective and cannot distinguish flux distributions with identical growth rate.

Pitfalls

• Unbounded fluxes: always cap exchange-reaction bounds; a model with unconstrained uptake gives infinite/unrealistic growth.
• Infeasible models (status != "optimal"): usually a missing biomass precursor or a broken GPR rule — check model.slim_optimize() and model.summary() before trusting results.
• Gene knockout artifacts: a single-gene knockout can show no growth defect purely because of isozyme redundancy in the GPR (OR rule) — inspect reaction.gene_reaction_rule before concluding a gene is non-essential.
• Alternate optima: FBA gives one optimal flux vector among possibly many; always pair with FVA before claiming a specific pathway is "the" flux route.
• Medium mismatches: forgetting to close default exchange reactions before opening your intended medium silently leaves unintended nutrients available.

See Also

• bio-metabolomics-pathway-mapping
• bio-metabolomics-statistical-analysis
• cobrapy
• bio-machine-learning-biomarker-discovery