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)
valid = results[results.status.eq("optimal") & results.growth.notna()]
low_growth = valid[valid.growth < 0.01 * baseline]  # declared 1% criterion
unresolved = results[~results.index.isin(valid.index)]
print(low_growth[["ids", "growth"]], unresolved[["ids", "status"]])

with model:
    medium = model.medium
    medium["EX_o2_e"] = 0.0
    model.medium = medium  # editing the returned dict alone does not change model
    anaerobic = model.optimize()
    print(anaerobic.status, anaerobic.objective_value)

min_medium = minimal_medium(model, 0.5 * baseline, minimize_components=True)
if min_medium is None:
    raise RuntimeError("No medium found at the requested growth target")
with model:
    model.medium = min_medium.to_dict()
    assert model.slim_optimize(error_value=None) >= 0.5 * baseline - 1e-6

Deletion results contain ids sets, growth, and status; the index is not a pair MultiIndex. Double-deletion results also contain singleton sets (self-pairs). Classify failed/infeasible solves separately from feasible low-growth mutants. A knockout's growth column means the model's objective, so ensure it is biomass.

model.medium values are positive import bounds, not measured concentrations. Exchange flux signs depend on reaction stoichiometry; the textbook's one-reactant exchanges use negative flux for uptake. open_exchanges=False retains the allowed nutrient universe while the optimizer selects a nutrient subset and its import bounds; it does not fix the selected nutrients or amounts. open_exchanges=True expands that universe and can choose unintended carbon sources. Minimal media can be nonunique; validate the returned medium with the intended growth objective.

4. Sample the same feasible region as the FVA comparison

from cobra.sampling import OptGPSampler

with model:
    model.reactions.get_by_id(biomass_id).lower_bound = 0.9 * baseline
    # The growth floor is already installed. fraction=0 adds no stronger optimum.
    sampled_fva = flux_variability_analysis(
        model, reaction_list=["PFK"], fraction_of_optimum=0.0, processes=1,
    )
    sampler = OptGPSampler(model, processes=1, thinning=100, seed=7)
    samples = sampler.sample(200)
    codes = sampler.validate(samples)
    assert (codes == "v").all(), "Inspect bound/equality violations"
    assert (samples[biomass_id] >= 0.9 * baseline - 1e-6).all()
print(samples["PFK"].describe(), sampled_fva)

FVA does not leave its objective-fraction constraint on the model; sampling a fresh/unconstrained model explores a different space. Feasible samples are not a confidence interval for measured biology. Inspect independent chains, autocorrelation and effective sample size before interpreting distributions; 200 samples are a smoke test. Sampling can include internal cycles.

5. Production envelopes and design hypotheses

from cobra.flux_analysis import production_envelope

with model:
    model.reactions.get_by_id(biomass_id).lower_bound = 0.1 * baseline
    envelope = production_envelope(
        model, reactions=["EX_glc__D_e"], objective="EX_ac_e",
        carbon_sources=["EX_glc__D_e"], points=8,
    )
print(envelope[["EX_glc__D_e", "flux_minimum", "flux_maximum"]])

These are acetate flux extrema; carbon_yield_* and mass_yield_* are separate outputs, potentially NaN when inputs or formulas are unsuitable. A multi-reaction grid needs a surface/heatmap, not an arbitrary connected line. When zero flux is already allowed for the affected reactions, a knockout only removes feasible states and cannot improve the global product maximum under otherwise identical constraints. Check that knockout bound replacement does not relax an original forced nonzero flux. The design workflow screens the minimum product flux at a common growth requirement and also reports maximum growth, an explicit model-based coupling hypothesis requiring validation.

6. Build and gap-fill carefully

Use Model, Reaction, Metabolite, model.add_reactions([reaction]), and reaction.gene_reaction_rule to build the network. Set formulas and charges; include water/protons when needed for balanced chemistry. Exchanges belong to external metabolites, demands remove metabolites, and sinks permit both directions. Do not add arbitrary ATP sources to make a model grow.

The API reference includes a balanced toy network and a gap-fill example with a known missing reaction. gapfill returns a list of reaction lists, one per iteration; it does not modify the input model. Candidate additions are hypotheses: check evidence, directionality, chemistry, energy-generating cycles, and growth after adding copied reactions.

Export and troubleshooting

Prefer SBML for exchange, JSON for interoperable tooling, YAML for inspection. Round-trip the file and compare reaction IDs, bounds, objective, and growth. Use context managers for temporary objective/bound/GPR changes. Inspect statuses before reading fluxes; do not turn NaN or every solver failure into zero growth. When debugging infeasibility, test medium changes inside with model: and record which constraints changed. Feasibility after opening all nutrients does not establish biological validity. Keep CSV/PNG output in the task's chosen directory.

• API contracts and current download routes
• Executed local workflows and their interpretation
• Release notes