# Bio Applied Metabolic Flux

> Run flux balance analysis (FBA/FVA) on genome-scale metabolic models (E. coli core, Recon3D, AGORA2) with COBRApy; simulate single/double gene knockouts and integrate RNA-seq expression via GIMME/iMAT. Use when predicting metabolic fluxes, finding essential genes or drug targets, doing synthetic lethality screens, or building transcriptomics-constrained metabolic models from SBML/JSON GEMs.

- Skill: `pavel-kravchenko/bio-applied-metabolic-flux` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-metabolic-flux`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-metabolic-flux/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

```bash
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.

```python
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.

```python
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.

```python
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.

```python
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