| name | cobrapy |
| description | Constraint-based metabolic modeling (COBRA). FBA, FVA, gene knockouts, flux sampling, SBML models, for systems biology and metabolic engineering analysis. |
| license | GPL-2.0 license |
| allowed-tools | Read Write Edit Bash |
| compatibility | Requires Python 3.9+ (cobra 0.30+ dropped 3.8). Install with uv pip install. GLPK (swiglpk) is the default solver; CPLEX/Gurobi optional. load_model fetches from bundled data, BiGG, or BioModels (network required for remote models). |
| metadata | version: "1.1" skill-author: K-Dense Inc. |
COBRApy - Constraint-Based Reconstruction and Analysis
Overview
COBRApy is a Python library for constraint-based reconstruction and analysis (COBRA) of metabolic models, essential for systems biology research. Work with genome-scale metabolic models, perform computational simulations of cellular metabolism, conduct metabolic engineering analyses, and predict phenotypic behaviors.
Version note: Examples target cobra 0.31.1 on PyPI (import cobra). Docs: cobrapy.readthedocs.io. Repo: opencobra/cobrapy.
When to Use This Skill
Use this skill when:
- Loading, building, or exporting genome-scale metabolic models (SBML, JSON, YAML)
- Running FBA, pFBA, FVA, or flux sampling on COBRA models
- Performing gene or reaction knockout screens and production envelope analysis
- Designing or optimizing growth media and exchange constraints
- Gap-filling infeasible models or validating model consistency
Installation
uv pip install "cobra==0.31.1"
MATLAB model I/O (optional):
uv pip install "cobra[array]==0.31.1"
COBRApy uses optlang for solvers. GLPK installs automatically via swiglpk. For large MILPs/QPs, cobra 0.29+ adds a hybrid solver (HIGHS/OSQP); model.solver = "osqp" now routes through hybrid and may error on plain LPs in a future release—prefer model.solver = "hybrid" when available.
Core Capabilities
COBRApy provides comprehensive tools organized into several key areas:
1. Model Management
Load existing models from repositories or files:
from cobra.io import load_model
model = load_model("textbook")
model = load_model("e_coli_core")
model = load_model("iJO1366")
model = load_model("salmonella")
model = load_model("iML1515")
from cobra.io import read_sbml_model, load_json_model, load_yaml_model
model = read_sbml_model("path/to/model.xml")
model = load_json_model("path/to/model.json")
model = load_yaml_model("path/to/model.yml")
Save models in various formats:
from cobra.io import write_sbml_model, save_json_model, save_yaml_model
write_sbml_model(model, "output.xml")
save_json_model(model, "output.json")
save_yaml_model(model, "output.yml")
2. Model Structure and Components
Access and inspect model components:
model.reactions
model.metabolites
model.genes
reaction = model.reactions.get_by_id("PFK")
metabolite = model.metabolites[0]
print(reaction.reaction)
print(reaction.bounds)
print(reaction.gene_reaction_rule)
print(metabolite.formula)
print(metabolite.compartment)
3. Flux Balance Analysis (FBA)
Perform standard FBA simulation:
solution = model.optimize()
print(f"Objective value: {solution.objective_value}")
print(f"Status: {solution.status}")
print(solution.fluxes["PFK"])
print(solution.fluxes.head())
objective_value = model.slim_optimize()
model.objective = "ATPM"
solution = model.optimize()
Parsimonious FBA (minimize total flux):
from cobra.flux_analysis import pfba
solution = pfba(model)
Geometric FBA (find central solution):
from cobra.flux_analysis import geometric_fba
solution = geometric_fba(model)
4. Flux Variability Analysis (FVA)
Determine flux ranges for all reactions:
from cobra.flux_analysis import flux_variability_analysis
fva_result = flux_variability_analysis(model)
fva_result = flux_variability_analysis(model, fraction_of_optimum=0.9)
fva_result = flux_variability_analysis(model, loopless=True)
fva_result = flux_variability_analysis(
model,
reaction_list=["PFK", "FBA", "PGI"]
)
5. Gene and Reaction Deletion Studies
Perform knockout analyses:
from cobra.flux_analysis import (
single_gene_deletion,
single_reaction_deletion,
double_gene_deletion,
double_reaction_deletion
)
gene_results = single_gene_deletion(model)
reaction_results = single_reaction_deletion(model)
double_gene_results = double_gene_deletion(
model,
processes=4
)
with model:
model.genes.get_by_id("b0008").knock_out()
solution = model.optimize()
print(f"Growth after knockout: {solution.objective_value}")
6. Growth Media and Minimal Media
Manage growth medium:
print(model.medium)
medium = model.medium
medium["EX_glc__D_e"] = 10.0
medium["EX_o2_e"] = 0.0
model.medium = medium
from cobra.medium import minimal_medium
min_medium = minimal_medium(model, minimize_components=False)
min_medium = minimal_medium(
model,
minimize_components=True,
open_exchanges=True
)
7. Flux Sampling
Sample the feasible flux space:
from cobra.sampling import sample
samples = sample(model, n=1000, method="optgp", processes=4)
samples = sample(model, n=1000, method="achr")
from cobra.sampling import OptGPSampler
sampler = OptGPSampler(model, processes=4)
sampler.sample(1000)
validation = sampler.validate(sampler.samples)
print(validation.value_counts())
8. Production Envelopes
Calculate phenotype phase planes:
from cobra.flux_analysis import production_envelope
envelope = production_envelope(
model,
reactions=["EX_glc__D_e", "EX_o2_e"],
objective="EX_ac_e"
)
envelope = production_envelope(
model,
reactions=["EX_glc__D_e", "EX_o2_e"],
carbon_sources="EX_glc__D_e"
)
import matplotlib.pyplot as plt
envelope.plot(x="EX_glc__D_e", y="EX_o2_e", kind="scatter")
plt.show()
9. Gapfilling
Add reactions to make models feasible:
from cobra.flux_analysis import gapfill
from cobra.io import read_sbml_model
universal = read_sbml_model("path/to/universal_reactions.xml")
with model:
model.remove_reactions([model.reactions.PGI])
solution = gapfill(model, universal)
print(f"Reactions to add: {solution}")
10. Model Building
Build models from scratch:
from cobra import Model, Reaction, Metabolite
model = Model("my_model")
atp_c = Metabolite("atp_c", formula="C10H12N5O13P3",
name="ATP", compartment="c")
adp_c = Metabolite("adp_c", formula="C10H12N5O10P2",
name="ADP", compartment="c")
pi_c = Metabolite("pi_c", formula="HO4P",
name="Phosphate", compartment="c")
reaction = Reaction("ATPASE")
reaction.name = "ATP hydrolysis"
reaction.subsystem = "Energy"
reaction.lower_bound = 0.0
reaction.upper_bound = 1000.0
reaction.add_metabolites({
atp_c: -1.0,
adp_c: 1.0,
pi_c: 1.0
})
reaction.gene_reaction_rule = "(gene1 and gene2) or gene3"
model.add_reactions([reaction])
model.add_boundary(atp_c, type="exchange")
model.add_boundary(adp_c, type="demand")
model.objective = "ATPASE"
Common Workflows
Workflow 1: Load Model and Predict Growth
from cobra.io import load_model
model = load_model("textbook")
solution = model.optimize()
print(f"Growth rate: {solution.objective_value:.3f} /h")
print(solution.fluxes[solution.fluxes.abs() > 1e-6])
Workflow 2: Gene Knockout Screen
from cobra.io import load_model
from cobra.flux_analysis import single_gene_deletion
model = load_model("textbook")
baseline = model.slim_optimize()
results = single_gene_deletion(model)
essential_genes = results[results["growth"] < 0.01]
print(f"Found {len(essential_genes)} essential genes")
neutral_genes = results[results["growth"] > 0.9 * baseline]
Workflow 3: Media Optimization
from cobra.io import load_model
from cobra.medium import minimal_medium
model = load_model("textbook")
target_growth = model.slim_optimize() * 0.5
min_medium = minimal_medium(
model,
target_growth,
minimize_components=True
)
print(f"Minimal medium components: {len(min_medium)}")
print(min_medium)
Workflow 4: Flux Uncertainty Analysis
from cobra.io import load_model
from cobra.flux_analysis import flux_variability_analysis
from cobra.sampling import sample
model = load_model("textbook")
fva = flux_variability_analysis(model, fraction_of_optimum=1.0)
samples = sample(model, n=1000)
reaction_id = "PFK"
import matplotlib.pyplot as plt
samples[reaction_id].hist(bins=50)
plt.xlabel(f"Flux through {reaction_id}")
plt.ylabel("Frequency")
plt.show()
Workflow 5: Context Manager for Temporary Changes
Use context managers to make temporary modifications:
with model:
model.objective = "ATPM"
model.reactions.EX_glc__D_e.lower_bound = -5.0
model.genes.b0008.knock_out()
solution = model.optimize()
print(f"Modified growth: {solution.objective_value}")
solution = model.optimize()
print(f"Original growth: {solution.objective_value}")
Key Concepts
DictList Objects
Models use DictList objects for reactions, metabolites, and genes - behaving like both lists and dictionaries:
first_reaction = model.reactions[0]
pfk = model.reactions.get_by_id("PFK")
atp_reactions = model.reactions.query("atp")
Flux Constraints
Reaction bounds define feasible flux ranges:
- Irreversible:
lower_bound = 0, upper_bound > 0
- Reversible:
lower_bound < 0, upper_bound > 0
- Set both bounds simultaneously with
.bounds to avoid inconsistencies
Gene-Reaction Rules (GPR)
Boolean logic linking genes to reactions:
reaction.gene_reaction_rule = "gene1 and gene2"
reaction.gene_reaction_rule = "gene1 or gene2"
reaction.gene_reaction_rule = "(gene1 and gene2) or (gene3 and gene4)"
Exchange Reactions
Special reactions representing metabolite import/export:
- Named with prefix
EX_ by convention
- Positive flux = secretion, negative flux = uptake
- Managed through
model.medium dictionary
Best Practices
- Use context managers for temporary modifications to avoid state management issues
- Validate models before analysis using
model.slim_optimize() to ensure feasibility
- Check solution status after optimization -
optimal indicates successful solve
- Use loopless FVA when thermodynamic feasibility matters
- Set fraction_of_optimum appropriately in FVA to explore suboptimal space
- Parallelize computationally expensive operations (sampling, double deletions) — start with small
n and processes=1 on genome-scale models
- Prefer SBML format for model exchange and long-term storage
- Use slim_optimize() when only objective value needed for performance
- Validate flux samples to ensure numerical stability
- Confirm output paths before writing CSV/PNG files from workflow examples
Troubleshooting
Infeasible solutions: Check medium constraints, reaction bounds, and model consistency
Slow optimization: Try different solvers (GLPK, CPLEX, Gurobi) via model.solver
Unbounded solutions: Verify exchange reactions have appropriate upper bounds
Import errors: Ensure correct file format and valid SBML identifiers
References
For detailed workflows and API patterns, refer to:
references/workflows.md - Comprehensive step-by-step workflow examplesreferences/api_quick_reference.md - Common function signatures and patterns
Official documentation: https://cobrapy.readthedocs.io/en/latest/
