COBRApy - Constraint-Based Reconstruction and Analysis

When to use

Use for loading/building metabolic networks, optimizing their steady-state fluxes, knockout screens, medium design, feasible-space sampling, and gap-filling hypotheses. FBA is a constraint-based prediction; it does not infer kinetic rates or establish experimental growth, thermodynamic feasibility, or flux identifiability.

Setup and reproducibility

Targets cobra 0.32.1 (import cobra), checked against its released source and current documentation.

uv pip install "cobra==0.32.1"
# Optional SciPy support for MATLAB I/O and array operations:
uv pip install "cobra[array]==0.32.1"

Plotting examples additionally require matplotlib, pandas, and seaborn. The cobra[chrr] extra supplies hopsy for the new CHRR sampler; that optional backend was documentation-reviewed, not executed in this refresh. In 0.32.1, sample() defaults to method="auto" (CHRR when hopsy is installed, otherwise OptGP). Choose the method explicitly to avoid environment-dependent changes.

Record model source/version/checksum, cobra and solver versions, objective, medium, bounds, tolerances, and any random seed with each analysis. GLPK handles the examples below; inspect cobra.util.solver.solvers before choosing an optional solver. Use "hybrid" for its HiGHS/OSQP interface when installed; the legacy "osqp" alias is deprecated. QP methods require a suitable backend. Start with processes=1; scripts using multiprocessing need a guarded entry point.

Workflow

1. Load, inspect, and validate a model

from cobra.io import load_model

# These names are bundled: textbook, iJO1366, salmonella.
model = load_model("textbook")  # model.id is e_coli_core; 95 reactions
model.solver = "glpk"
print(model.id, len(model.reactions), len(model.metabolites), len(model.genes))
print(model.reactions.get_by_id("PFK").reaction)
print(model.reactions.PFK.gene_reaction_rule)
print(model.medium)

solution = model.optimize(raise_error=True)
assert solution.status == "optimal"
print(solution.objective_value, solution.fluxes["PFK"])
# error_value=None raises on a failed solve; the default instead returns NaN.
baseline = model.slim_optimize(error_value=None)
assert baseline > 0

e_coli_core is a remote BiGG identifier, not a bundled alias for textbook. The released remote adapters encounter redirects on current BiGG/BioModels URLs; see model I/O for verified download routes and local SBML loading. Do not assume load_model caches: its 0.32.1 implementation accepts cache but does not use it. Save source files for reproducibility.

Check chemical formulas/charges and boundary annotations separately from solver feasibility. Review excluded biomass/pseudo reactions and missing chemistry; an empty imbalance dictionary alone cannot establish a chemically valid model. The validation workflow distinguishes those cases.

2. Compare FBA, pFBA, and FVA under stated constraints

from cobra.flux_analysis import pfba, geometric_fba, flux_variability_analysis

biomass_id = "Biomass_Ecoli_core"  # inspect IDs/objective for each new model
parsimonious = pfba(model)
print(parsimonious.fluxes[biomass_id])
# parsimonious.objective_value is the minimized total flux, not biomass growth.
centered = geometric_fba(model, processes=1)

fva = flux_variability_analysis(
    model, reaction_list=["PFK", "FBA", "PGI"],
    fraction_of_optimum=0.9, processes=1,
)
loopless_fva = flux_variability_analysis(
    model, reaction_list=["PFK", "FBA", "PGI"],
    fraction_of_optimum=0.9, loopless="fastSNP", processes=1,
)
print(fva)  # index: reaction ID; columns: minimum, maximum

FVA extrema are optimized separately and need not be jointly achievable. loopless="fastSNP" computes loopless bounds; "cycleFreeFlux" is an alternative that need not find the optimal bounds. Boolean loopless arguments are deprecated. Loop removal does not impose measured Gibbs energies or metabolite concentrations. Fractional objective thresholds here assume a positive biomass maximization objective; use an explicit constraint for other objective signs/directions.

3. Knockouts and media

from cobra.flux_analysis import single_gene_deletion
from cobra.medium import minimal_medium

results = single_gene_deletion(model, processes=1)

…