--- name: alterlab-cobrapy description: Build and analyze genome-scale constraint-based metabolic models with COBRApy — flux balance analysis (FBA), flux variability analysis (FVA), gene and reaction knockouts, flux sampling, and SBML model I/O. Use when simulating metabolic networks, predicting growth or knockout phenotypes, or running systems-biology and metabolic-engineering analyses on SBML genome-scale models. Part of the AlterLab Academic Skills suite. license: GPL-2.0 allowed-tools: Read Write Edit Bash(python:*) Bash(uv:*) compatibility: "Self-contained — runs under `uv run python` with the skill's Python package installed; no API key or account required." metadata: skill-author: AlterLab version: "1.0.0" ---

# 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.

## Installation and Requirements

```bash uv pip install 'cobra>=0.29,<0.32' ```

**Requirements:** Python 3.8+; depends on optlang and a solver (GLPK ships by default via swiglpk; CPLEX/Gurobi optional). Optional for plots: matplotlib, seaborn, pandas (pandas is already a hard dependency). Flux sampling and parallel deletions use multiprocessing.

> **Bundled model names.** `load_model` only recognizes three bundled aliases: `"textbook"` (the E. coli core model, `e_coli_core`), `"iJO1366"` (full E. coli genome-scale), and `"salmonella"`. Any other identifier (e.g. `"ecoli"`, `"e_coli_core"`) is not a bundled alias and triggers a remote lookup against BiGG/BioModels, which requires network access and may fail. Use `read_sbml_model("path.xml")` for your own models.

## Core Capabilities

COBRApy provides comprehensive tools organized into several key areas:

### 1. Model Management

Load existing models from repositories or files: ```python from cobra.io import load_model

# Load bundled test models model = load_model("textbook") # E. coli core model model = load_model("iJO1366") # Full E. coli genome-scale model (BiGG) model = load_model("salmonella")

# Load from files 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: ```python from cobra.io import write_sbml_model, save_json_model, save_yaml_model write_sbml_model(model, "output.xml") # Preferred format save_json_model(model, "output.json") # For Escher compatibility save_yaml_model(model, "output.yml") # Human-readable ```

### 2. Model Structure and Components

Access and inspect model components: ```python # Access components model.reactions # DictList of all reactions model.metabolites # DictList of all metabolites model.genes # DictList of all genes

# Get specific items by ID or index reaction = model.reactions.get_by_id("PFK") metabolite = model.metabolites[0]

# Inspect properties print(reaction.reaction) # Stoichiometric equation print(reaction.bounds) # Flux constraints print(reaction.gene_reaction_rule) # GPR logic print(metabolite.formula) # Chemical formula print(metabolite.compartment) # Cellular location ```

### 3. Flux Balance Analysis (FBA)

Perform standard FBA simulation: ```python # Basic optimization solution = model.optimize() print(f"Objective value: {solution.objective_value}") print(f"Status: {solution.status}")

# Access fluxes print(solution.fluxes["PFK"]) print(solution.fluxes.head())

# Fast optimization (objective value only) objective_value = model.slim_optimize()

# Change objective model.objective = "ATPM" solution = model.optimize() ```

Parsimonious FBA (minimize total flux): ```python from cobra.flux_analysis import pfba solution = pfba(model) ```

Geometric FBA (find central solution): ```python from cobra.flux_analysis import geometric_fba solution = geometric_fba(model) ```

### 4. Flux Variability Analysis (FVA)

Determine flux ranges for all reactions: ```python from cobra.flux_analysis import flux_variability_analysis

# Standard FVA fva_result = flux_variability_analysis(model)

# FVA at 90% optimality fva_result = flux_variability_analysis(model, fraction_of_optimum=0.9)

# Loopless FVA (eliminates thermodynamically infeasible loops) fva_result = flux_variability_analysis(model, loopless=True)

# FVA for specific reactions fva_result = flux_variability_analysis( model, reaction_list=["PFK", "FBA", "PGI"] ) ```

### 5. Gene and Reaction Deletion Studies

Perform knockout analyses: ```python from cobra.flux_analysis import ( single_gene_deletion, single_reaction_deletion, double_gene_deletion, double_reaction_deletion )

# Single deletions gene_results = single_gene_deletion(model) reaction_results = single_reaction_deletion(model)

# Double deletions (uses multiprocessing) double_gene_results = double_gene_deletion( model, processes=4 # Number of CPU cores )

# Manual knockout using context manager with model: model.genes.get_by_id("b0008").knock_out() solution = model.optimize() print(f"Growth after knockout: {solution.objective_value}") # Model automatically reverts after context exit ```

### 6. Growth Media and Minimal Media

Manage growth medium: ```python # View current medium print(model.medium)

# Modify medium (must reassign entire dict) medium = model.medium medium["EX_glc__D_e"] = 10.0 # Set glucose uptake medium["EX_o2_e"] = 0.0 # Anaerobic conditions model.medium = medium

# Calculate minimal media from cobra.medium import minimal_medium

# Minimize total import flux min_medium = minimal_medium(model, minimize_components=False)

# Minimize number of components (uses MILP, slower) min_medium = minimal_medium( model, minimize_components=True, open_exchanges=True ) ```

### 7. Flux Sampling

Sample the feasible flux space: ```python from cobra.sampling import sample

# Sample using OptGP (default, supports parallel processing) samples = sample(model, n=1000, method="optgp", processes=4)

# Sample using ACHR samples = sample(model, n=1000, method="achr")

# Validate samples from cobra.sampling import OptGPSampler sampler = OptGPSampler(model, processes=4) sampler.sample(1000) validation = sampler.validate(sampler.samples) print(validation.value_counts()) # Should be all 'v' for valid ```

### 8. Production Envelopes

Calculate phenotype phase planes: ```python from cobra.flux_analysis import production_envelope

# Standard production envelope envelope = production_envelope( model, reactions=["EX_glc__D_e", "EX_o2_e"], objective="EX_ac_e" # Acetate production )

# With carbon yield envelope = production_envelope( model, reactions=["EX_glc__D_e", "EX_o2_e"], carbon_sources="EX_glc__D_e" )

# Visualize (use matplotlib or pandas plotting) 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: ```python import cobra from cobra.flux_analysis import gapfill

# Build a universal model of candidate reactions to draw from. # In practice, load a curated reaction database (e.g. the BiGG # universal reactions JSON via load_json_model) or assemble reactions # into an empty cobra.Model — there is no load_model("universal"). universal = cobra.Model("universal_reactions") # ...populate `universal` with candidate cobra.Reaction objects...

with model: # Create a gap for demonstration model.remove_reactions([model.reactions.PGI])

# Find the minimal set of reactions from `universal` to restore feasibility solution = gapfill(model, universal, demand_reactions=False) for reaction in solution[0]: print(reaction.id) ```

> **Note**: `gapfill` returns a list of reaction lists (one per iteration when `iterations > 1`), so iterate `solution[0]` rather than printing `solution` directly.

### 10. Model Building

Build models from scratch: ```python from cobra import Model, Reaction, Metabolite

# Create model model = Model("my_model")

# Create metabolites 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")

# Create reaction reaction = Reaction("ATPASE") reaction.name = "ATP hydrolysis" reaction.subsystem = "Energy" reaction.lower_bound = 0.0 reaction.upper_bound = 1000.0

# Add metabolites with stoichiometry reaction.add_metabolites({ atp_c: -1.0, adp_c: 1.0, pi_c: 1.0 })

# Add gene-reaction rule reaction.gene_reaction_rule = "(gene1 and gene2) or gene3"

# Add to model model.add_reactions([reaction])

# Add boundary reactions model.add_boundary(atp_c, type="exchange") model.add_boundary(adp_c, type="demand")

# Set objective model.objective = "ATPASE" ```

## Common Workflows

### Workflow 1: Load Model and Predict Growth

```python from cobra.io import load_model

# Load model model = load_model("iJO1366")

# Run FBA solution = model.optimize() print(f"Growth rate: {solution.objective_value:.3f} /h")

# Show active pathways print(solution.fluxes[solution.fluxes.abs() > 1e-6]) ```

### Workflow 2: Gene Knockout Screen

```python from cobra.io import load_model from cobra.flux_analysis import single_gene_deletion

# Load model model = load_model("iJO1366")

# Perform single gene deletions results = single_gene_deletion(model)

# Find essential genes (growth < threshold) essential_genes = results[results["growth"] < 0.01] print(f"Found {len(essential_genes)} essential genes")

# Find genes with minimal impact neutral_genes = results[results["growth"] > 0.9 * solution.objective_value] ```

### Workflow 3: Media Optimization

```python from cobra.io import load_model from cobra.medium import minimal_medium

# Load model model = load_model("iJO1366")

# Calculate minimal medium for 50% of max growth 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

```python from cobra.io import load_model from cobra.flux_analysis import flux_variability_analysis from cobra.sampling import sample

# Load model model = load_model("iJO1366")

# First check flux ranges at optimality fva = flux_variability_analysis(model, fraction_of_optimum=1.0)

# For reactions with large ranges, sample to understand distribution samples = sample(model, n=1000)

# Analyze specific reaction 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: ```python # Model remains unchanged outside context with model: # Temporarily change objective model.objective = "ATPM"

# Temporarily modify bounds model.reactions.EX_glc__D_e.lower_bound = -5.0

# Temporarily knock out genes model.genes.b0008.knock_out()

# Optimize with changes solution = model.optimize() print(f"Modified growth: {solution.objective_value}")

# All changes automatically reverted 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: ```python # Access by index first_reaction = model.reactions[0]

# Access by ID pfk

--- name: alterlab-cobrapy description: Build and analyze genome-scale constraint-based metabolic models with COBRApy — flux balance analysis (FBA), flux variability analysis (FVA), gene and reaction knockouts, flux sampling, and SBML model I/O. Use when simulating metabolic networks, predicting growth or knockout phenotypes, or running systems-biology and metabolic-engineering analyses on SBML genome-scale models. Part of the AlterLab Academic Skills suite. license: GPL-2.0 allowed-tools: Read Write Edit Bash(python:*) Bash(uv:*) compatibility: "Self-contained — runs under `uv run python` with the skill's Python package installed; no API key or account required." metadata: skill-author: AlterLab version: "1.0.0" ---

# 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.

## Installation and Requirements

```bash uv pip install 'cobra>=0.29,<0.32' ```

**Requirements:** Python 3.8+; depends on optlang and a solver (GLPK ships by default via swiglpk; CPLEX/Gurobi optional). Optional for plots: matplotlib, seaborn, pandas (pandas is already a hard dependency). Flux sampling and parallel deletions use multiprocessing.

> **Bundled model names.** `load_model` only recognizes three bundled aliases: `"textbook"` (the E. coli core model, `e_coli_core`), `"iJO1366"` (full E. coli genome-scale), and `"salmonella"`. Any other identifier (e.g. `"ecoli"`, `"e_coli_core"`) is not a bundled alias and triggers a remote lookup against BiGG/BioModels, which requires network access and may fail. Use `read_sbml_model("path.xml")` for your own models.

## Core Capabilities

COBRApy provides comprehensive tools organized into several key areas:

### 1. Model Management

Load existing models from repositories or files: ```python from cobra.io import load_model

# Load bundled test models model = load_model("textbook") # E. coli core model model = load_model("iJO1366") # Full E. coli genome-scale model (BiGG) model = load_model("salmonella")

# Load from files 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: ```python from cobra.io import write_sbml_model, save_json_model, save_yaml_model write_sbml_model(model, "output.xml") # Preferred format save_json_model(model, "output.json") # For Escher compatibility save_yaml_model(model, "output.yml") # Human-readable ```

### 2. Model Structure and Components

Access and inspect model components: ```python # Access components model.reactions # DictList of all reactions model.metabolites # DictList of all metabolites model.genes # DictList of all genes

# Get specific items by ID or index reaction = model.reactions.get_by_id("PFK") metabolite = model.metabolites[0]

# Inspect properties print(reaction.reaction) # Stoichiometric equation print(reaction.bounds) # Flux constraints print(reaction.gene_reaction_rule) # GPR logic print(metabolite.formula) # Chemical formula print(metabolite.compartment) # Cellular location ```

### 3. Flux Balance Analysis (FBA)

Perform standard FBA simulation: ```python # Basic optimization solution = model.optimize() print(f"Objective value: {solution.objective_value}") print(f"Status: {solution.status}")

# Access fluxes print(solution.fluxes["PFK"]) print(solution.fluxes.head())

# Fast optimization (objective value only) objective_value = model.slim_optimize()

# Change objective model.objective = "ATPM" solution = model.optimize() ```

Parsimonious FBA (minimize total flux): ```python from cobra.flux_analysis import pfba solution = pfba(model) ```

Geometric FBA (find central solution): ```python from cobra.flux_analysis import geometric_fba solution = geometric_fba(model) ```

### 4. Flux Variability Analysis (FVA)

Determine flux ranges for all reactions: ```python from cobra.flux_analysis import flux_variability_analysis

# Standard FVA fva_result = flux_variability_analysis(model)

# FVA at 90% optimality fva_result = flux_variability_analysis(model, fraction_of_optimum=0.9)

# Loopless FVA (eliminates thermodynamically infeasible loops) fva_result = flux_variability_analysis(model, loopless=True)

# FVA for specific reactions fva_result = flux_variability_analysis( model, reaction_list=["PFK", "FBA", "PGI"] ) ```

### 5. Gene and Reaction Deletion Studies

Perform knockout analyses: ```python from cobra.flux_analysis import ( single_gene_deletion, single_reaction_deletion, double_gene_deletion, double_reaction_deletion )

# Single deletions gene_results = single_gene_deletion(model) reaction_results = single_reaction_deletion(model)

# Double deletions (uses multiprocessing) double_gene_results = double_gene_deletion( model, processes=4 # Number of CPU cores )

# Manual knockout using context manager with model: model.genes.get_by_id("b0008").knock_out() solution = model.optimize() print(f"Growth after knockout: {solution.objective_value}") # Model automatically reverts after context exit ```

### 6. Growth Media and Minimal Media

Manage growth medium: ```python # View current medium print(model.medium)

# Modify medium (must reassign entire dict) medium = model.medium medium["EX_glc__D_e"] = 10.0 # Set glucose uptake medium["EX_o2_e"] = 0.0 # Anaerobic conditions model.medium = medium

# Calculate minimal media from cobra.medium import minimal_medium

# Minimize total import flux min_medium = minimal_medium(model, minimize_components=False)

# Minimize number of components (uses MILP, slower) min_medium = minimal_medium( model, minimize_components=True, open_exchanges=True ) ```

### 7. Flux Sampling

Sample the feasible flux space: ```python from cobra.sampling import sample

# Sample using OptGP (default, supports parallel processing) samples = sample(model, n=1000, method="optgp", processes=4)

# Sample using ACHR samples = sample(model, n=1000, method="achr")

# Validate samples from cobra.sampling import OptGPSampler sampler = OptGPSampler(model, processes=4) sampler.sample(1000) validation = sampler.validate(sampler.samples) print(validation.value_counts()) # Should be all 'v' for valid ```

### 8. Production Envelopes

Calculate phenotype phase planes: ```python from cobra.flux_analysis import production_envelope

# Standard production envelope envelope = production_envelope( model, reactions=["EX_glc__D_e", "EX_o2_e"], objective="EX_ac_e" # Acetate production )

# With carbon yield envelope = production_envelope( model, reactions=["EX_glc__D_e", "EX_o2_e"], carbon_sources="EX_glc__D_e" )

# Visualize (use matplotlib or pandas plotting) 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: ```python import cobra from cobra.flux_analysis import gapfill

# Build a universal model of candidate reactions to draw from. # In practice, load a curated reaction database (e.g. the BiGG # universal reactions JSON via load_json_model) or assemble reactions # into an empty cobra.Model — there is no load_model("universal"). universal = cobra.Model("universal_reactions") # ...populate `universal` with candidate cobra.Reaction objects...

with model: # Create a gap for demonstration model.remove_reactions([model.reactions.PGI])

# Find the minimal set of reactions from `universal` to restore feasibility solution = gapfill(model, universal, demand_reactions=False) for reaction in solution[0]: print(reaction.id) ```

> **Note**: `gapfill` returns a list of reaction lists (one per iteration when `iterations > 1`), so iterate `solution[0]` rather than printing `solution` directly.

### 10. Model Building

Build models from scratch: ```python from cobra import Model, Reaction, Metabolite

# Create model model = Model("my_model")

# Create metabolites 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")

# Create reaction reaction = Reaction("ATPASE") reaction.name = "ATP hydrolysis" reaction.subsystem = "Energy" reaction.lower_bound = 0.0 reaction.upper_bound = 1000.0

# Add metabolites with stoichiometry reaction.add_metabolites({ atp_c: -1.0, adp_c: 1.0, pi_c: 1.0 })

# Add gene-reaction rule reaction.gene_reaction_rule = "(gene1 and gene2) or gene3"

# Add to model model.add_reactions([reaction])

# Add boundary reactions model.add_boundary(atp_c, type="exchange") model.add_boundary(adp_c, type="demand")

# Set objective model.objective = "ATPASE" ```

## Common Workflows

### Workflow 1: Load Model and Predict Growth

```python from cobra.io import load_model

# Load model model = load_model("iJO1366")

# Run FBA solution = model.optimize() print(f"Growth rate: {solution.objective_value:.3f} /h")

# Show active pathways print(solution.fluxes[solution.fluxes.abs() > 1e-6]) ```

### Workflow 2: Gene Knockout Screen

```python from cobra.io import load_model from cobra.flux_analysis import single_gene_deletion

# Load model model = load_model("iJO1366")

# Perform single gene deletions results = single_gene_deletion(model)

# Find essential genes (growth < threshold) essential_genes = results[results["growth"] < 0.01] print(f"Found {len(essential_genes)} essential genes")

# Find genes with minimal impact neutral_genes = results[results["growth"] > 0.9 * solution.objective_value] ```

### Workflow 3: Media Optimization

```python from cobra.io import load_model from cobra.medium import minimal_medium

# Load model model = load_model("iJO1366")

# Calculate minimal medium for 50% of max growth 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

```python from cobra.io import load_model from cobra.flux_analysis import flux_variability_analysis from cobra.sampling import sample

# Load model model = load_model("iJO1366")

# First check flux ranges at optimality fva = flux_variability_analysis(model, fraction_of_optimum=1.0)

# For reactions with large ranges, sample to understand distribution samples = sample(model, n=1000)

# Analyze specific reaction 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: ```python # Model remains unchanged outside context with model: # Temporarily change objective model.objective = "ATPM"

# Temporarily modify bounds model.reactions.EX_glc__D_e.lower_bound = -5.0

# Temporarily knock out genes model.genes.b0008.knock_out()

# Optimize with changes solution = model.optimize() print(f"Modified growth: {solution.objective_value}")

# All changes automatically reverted 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: ```python # Access by index first_reaction = model.reactions[0]

# Access by ID pfk

--- name: alterlab-cobrapy description: Build and analyze genome-scale constraint-based metabolic models with COBRApy — flux balance analysis (FBA), flux variability analysis (FVA), gene and reaction knockouts, flux sampling, and SBML model I/O. Use when simulating metabolic networks, predicting growth or knockout phenotypes, or running systems-biology and metabolic-engineering analyses on SBML genome-scale models. Part of the AlterLab Academic Skills suite. license: GPL-2.0 allowed-tools: Read Write Edit Bash(python:*) Bash(uv:*) compatibility: "Self-contained — runs under `uv run python` with the skill's Python package installed; no API key or account required." metadata: skill-author: AlterLab version: "1.0.0" ---

# 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.

## Installation and Requirements

```bash uv pip install 'cobra>=0.29,<0.32' ```

**Requirements:** Python 3.8+; depends on optlang and a solver (GLPK ships by default via swiglpk; CPLEX/Gurobi optional). Optional for plots: matplotlib, seaborn, pandas (pandas is already a hard dependency). Flux sampling and parallel deletions use multiprocessing.

> **Bundled model names.** `load_model` only recognizes three bundled aliases: `"textbook"` (the E. coli core model, `e_coli_core`), `"iJO1366"` (full E. coli genome-scale), and `"salmonella"`. Any other identifier (e.g. `"ecoli"`, `"e_coli_core"`) is not a bundled alias and triggers a remote lookup against BiGG/BioModels, which requires network access and may fail. Use `read_sbml_model("path.xml")` for your own models.

## Core Capabilities

COBRApy provides comprehensive tools organized into several key areas:

### 1. Model Management

Load existing models from repositories or files: ```python from cobra.io import load_model

# Load bundled test models model = load_model("textbook") # E. coli core model model = load_model("iJO1366") # Full E. coli genome-scale model (BiGG) model = load_model("salmonella")

# Load from files 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: ```python from cobra.io import write_sbml_model, save_json_model, save_yaml_model write_sbml_model(model, "output.xml") # Preferred format save_json_model(model, "output.json") # For Escher compatibility save_yaml_model(model, "output.yml") # Human-readable ```

### 2. Model Structure and Components

Access and inspect model components: ```python # Access components model.reactions # DictList of all reactions model.metabolites # DictList of all metabolites model.genes # DictList of all genes

# Get specific items by ID or index reaction = model.reactions.get_by_id("PFK") metabolite = model.metabolites[0]

# Inspect properties print(reaction.reaction) # Stoichiometric equation print(reaction.bounds) # Flux constraints print(reaction.gene_reaction_rule) # GPR logic print(metabolite.formula) # Chemical formula print(metabolite.compartment) # Cellular location ```

### 3. Flux Balance Analysis (FBA)

Perform standard FBA simulation: ```python # Basic optimization solution = model.optimize() print(f"Objective value: {solution.objective_value}") print(f"Status: {solution.status}")

# Access fluxes print(solution.fluxes["PFK"]) print(solution.fluxes.head())

# Fast optimization (objective value only) objective_value = model.slim_optimize()

# Change objective model.objective = "ATPM" solution = model.optimize() ```

Parsimonious FBA (minimize total flux): ```python from cobra.flux_analysis import pfba solution = pfba(model) ```

Geometric FBA (find central solution): ```python from cobra.flux_analysis import geometric_fba solution = geometric_fba(model) ```

### 4. Flux Variability Analysis (FVA)

Determine flux ranges for all reactions: ```python from cobra.flux_analysis import flux_variability_analysis

# Standard FVA fva_result = flux_variability_analysis(model)

# FVA at 90% optimality fva_result = flux_variability_analysis(model, fraction_of_optimum=0.9)

# Loopless FVA (eliminates thermodynamically infeasible loops) fva_result = flux_variability_analysis(model, loopless=True)

# FVA for specific reactions fva_result = flux_variability_analysis( model, reaction_list=["PFK", "FBA", "PGI"] ) ```

### 5. Gene and Reaction Deletion Studies

Perform knockout analyses: ```python from cobra.flux_analysis import ( single_gene_deletion, single_reaction_deletion, double_gene_deletion, double_reaction_deletion )

# Single deletions gene_results = single_gene_deletion(model) reaction_results = single_reaction_deletion(model)

# Double deletions (uses multiprocessing) double_gene_results = double_gene_deletion( model, processes=4 # Number of CPU cores )

# Manual knockout using context manager with model: model.genes.get_by_id("b0008").knock_out() solution = model.optimize() print(f"Growth after knockout: {solution.objective_value}") # Model automatically reverts after context exit ```

### 6. Growth Media and Minimal Media

Manage growth medium: ```python # View current medium print(model.medium)

# Modify medium (must reassign entire dict) medium = model.medium medium["EX_glc__D_e"] = 10.0 # Set glucose uptake medium["EX_o2_e"] = 0.0 # Anaerobic conditions model.medium = medium

# Calculate minimal media from cobra.medium import minimal_medium

# Minimize total import flux min_medium = minimal_medium(model, minimize_components=False)

# Minimize number of components (uses MILP, slower) min_medium = minimal_medium( model, minimize_components=True, open_exchanges=True ) ```

### 7. Flux Sampling

Sample the feasible flux space: ```python from cobra.sampling import sample

# Sample using OptGP (default, supports parallel processing) samples = sample(model, n=1000, method="optgp", processes=4)

# Sample using ACHR samples = sample(model, n=1000, method="achr")

# Validate samples from cobra.sampling import OptGPSampler sampler = OptGPSampler(model, processes=4) sampler.sample(1000) validation = sampler.validate(sampler.samples) print(validation.value_counts()) # Should be all 'v' for valid ```

### 8. Production Envelopes

Calculate phenotype phase planes: ```python from cobra.flux_analysis import production_envelope

# Standard production envelope envelope = production_envelope( model, reactions=["EX_glc__D_e", "EX_o2_e"], objective="EX_ac_e" # Acetate production )

# With carbon yield envelope = production_envelope( model, reactions=["EX_glc__D_e", "EX_o2_e"], carbon_sources="EX_glc__D_e" )

# Visualize (use matplotlib or pandas plotting) 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: ```python import cobra from cobra.flux_analysis import gapfill

# Build a universal model of candidate reactions to draw from. # In practice, load a curated reaction database (e.g. the BiGG # universal reactions JSON via load_json_model) or assemble reactions # into an empty cobra.Model — there is no load_model("universal"). universal = cobra.Model("universal_reactions") # ...populate `universal` with candidate cobra.Reaction objects...

with model: # Create a gap for demonstration model.remove_reactions([model.reactions.PGI])

# Find the minimal set of reactions from `universal` to restore feasibility solution = gapfill(model, universal, demand_reactions=False) for reaction in solution[0]: print(reaction.id) ```

> **Note**: `gapfill` returns a list of reaction lists (one per iteration when `iterations > 1`), so iterate `solution[0]` rather than printing `solution` directly.

### 10. Model Building

Build models from scratch: ```python from cobra import Model, Reaction, Metabolite

# Create model model = Model("my_model")

# Create metabolites 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")

# Create reaction reaction = Reaction("ATPASE") reaction.name = "ATP hydrolysis" reaction.subsystem = "Energy" reaction.lower_bound = 0.0 reaction.upper_bound = 1000.0

# Add metabolites with stoichiometry reaction.add_metabolites({ atp_c: -1.0, adp_c: 1.0, pi_c: 1.0 })

# Add gene-reaction rule reaction.gene_reaction_rule = "(gene1 and gene2) or gene3"

# Add to model model.add_reactions([reaction])

# Add boundary reactions model.add_boundary(atp_c, type="exchange") model.add_boundary(adp_c, type="demand")

# Set objective model.objective = "ATPASE" ```

## Common Workflows

### Workflow 1: Load Model and Predict Growth

```python from cobra.io import load_model

# Load model model = load_model("iJO1366")

# Run FBA solution = model.optimize() print(f"Growth rate: {solution.objective_value:.3f} /h")

# Show active pathways print(solution.fluxes[solution.fluxes.abs() > 1e-6]) ```

### Workflow 2: Gene Knockout Screen

```python from cobra.io import load_model from cobra.flux_analysis import single_gene_deletion

# Load model model = load_model("iJO1366")

# Perform single gene deletions results = single_gene_deletion(model)

# Find essential genes (growth < threshold) essential_genes = results[results["growth"] < 0.01] print(f"Found {len(essential_genes)} essential genes")

# Find genes with minimal impact neutral_genes = results[results["growth"] > 0.9 * solution.objective_value] ```

### Workflow 3: Media Optimization

```python from cobra.io import load_model from cobra.medium import minimal_medium

# Load model model = load_model("iJO1366")

# Calculate minimal medium for 50% of max growth 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

```python from cobra.io import load_model from cobra.flux_analysis import flux_variability_analysis from cobra.sampling import sample

# Load model model = load_model("iJO1366")

# First check flux ranges at optimality fva = flux_variability_analysis(model, fraction_of_optimum=1.0)

# For reactions with large ranges, sample to understand distribution samples = sample(model, n=1000)

# Analyze specific reaction 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: ```python # Model remains unchanged outside context with model: # Temporarily change objective model.objective = "ATPM"

# Temporarily modify bounds model.reactions.EX_glc__D_e.lower_bound = -5.0

# Temporarily knock out genes model.genes.b0008.knock_out()

# Optimize with changes solution = model.optimize() print(f"Modified growth: {solution.objective_value}")

# All changes automatically reverted 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: ```python # Access by index first_reaction = model.reactions[0]

# Access by ID pfk

--- name: alterlab-cobrapy description: Build and analyze genome-scale constraint-based metabolic models with COBRApy — flux balance analysis (FBA), flux variability analysis (FVA), gene and reaction knockouts, flux sampling, and SBML model I/O. Use when simulating metabolic networks, predicting growth or knockout phenotypes, or running systems-biology and metabolic-engineering analyses on SBML genome-scale models. Part of the AlterLab Academic Skills suite. license: GPL-2.0 allowed-tools: Read Write Edit Bash(python:*) Bash(uv:*) compatibility: "Self-contained — runs under `uv run python` with the skill's Python package installed; no API key or account required." metadata: skill-author: AlterLab version: "1.0.0" ---

# 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.

## Installation and Requirements

```bash uv pip install 'cobra>=0.29,<0.32' ```

**Requirements:** Python 3.8+; depends on optlang and a solver (GLPK ships by default via swiglpk; CPLEX/Gurobi optional). Optional for plots: matplotlib, seaborn, pandas (pandas is already a hard dependency). Flux sampling and parallel deletions use multiprocessing.

> **Bundled model names.** `load_model` only recognizes three bundled aliases: `"textbook"` (the E. coli core model, `e_coli_core`), `"iJO1366"` (full E. coli genome-scale), and `"salmonella"`. Any other identifier (e.g. `"ecoli"`, `"e_coli_core"`) is not a bundled alias and triggers a remote lookup against BiGG/BioModels, which requires network access and may fail. Use `read_sbml_model("path.xml")` for your own models.

## Core Capabilities

COBRApy provides comprehensive tools organized into several key areas:

### 1. Model Management

Load existing models from repositories or files: ```python from cobra.io import load_model

# Load bundled test models model = load_model("textbook") # E. coli core model model = load_model("iJO1366") # Full E. coli genome-scale model (BiGG) model = load_model("salmonella")

# Load from files 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: ```python from cobra.io import write_sbml_model, save_json_model, save_yaml_model write_sbml_model(model, "output.xml") # Preferred format save_json_model(model, "output.json") # For Escher compatibility save_yaml_model(model, "output.yml") # Human-readable ```

### 2. Model Structure and Components

Access and inspect model components: ```python # Access components model.reactions # DictList of all reactions model.metabolites # DictList of all metabolites model.genes # DictList of all genes

# Get specific items by ID or index reaction = model.reactions.get_by_id("PFK") metabolite = model.metabolites[0]

# Inspect properties print(reaction.reaction) # Stoichiometric equation print(reaction.bounds) # Flux constraints print(reaction.gene_reaction_rule) # GPR logic print(metabolite.formula) # Chemical formula print(metabolite.compartment) # Cellular location ```

### 3. Flux Balance Analysis (FBA)

Perform standard FBA simulation: ```python # Basic optimization solution = model.optimize() print(f"Objective value: {solution.objective_value}") print(f"Status: {solution.status}")

# Access fluxes print(solution.fluxes["PFK"]) print(solution.fluxes.head())

# Fast optimization (objective value only) objective_value = model.slim_optimize()

# Change objective model.objective = "ATPM" solution = model.optimize() ```

Parsimonious FBA (minimize total flux): ```python from cobra.flux_analysis import pfba solution = pfba(model) ```

Geometric FBA (find central solution): ```python from cobra.flux_analysis import geometric_fba solution = geometric_fba(model) ```

### 4. Flux Variability Analysis (FVA)

Determine flux ranges for all reactions: ```python from cobra.flux_analysis import flux_variability_analysis

# Standard FVA fva_result = flux_variability_analysis(model)

# FVA at 90% optimality fva_result = flux_variability_analysis(model, fraction_of_optimum=0.9)

# Loopless FVA (eliminates thermodynamically infeasible loops) fva_result = flux_variability_analysis(model, loopless=True)

# FVA for specific reactions fva_result = flux_variability_analysis( model, reaction_list=["PFK", "FBA", "PGI"] ) ```

### 5. Gene and Reaction Deletion Studies

Perform knockout analyses: ```python from cobra.flux_analysis import ( single_gene_deletion, single_reaction_deletion, double_gene_deletion, double_reaction_deletion )

# Single deletions gene_results = single_gene_deletion(model) reaction_results = single_reaction_deletion(model)

# Double deletions (uses multiprocessing) double_gene_results = double_gene_deletion( model, processes=4 # Number of CPU cores )

# Manual knockout using context manager with model: model.genes.get_by_id("b0008").knock_out() solution = model.optimize() print(f"Growth after knockout: {solution.objective_value}") # Model automatically reverts after context exit ```

### 6. Growth Media and Minimal Media

Manage growth medium: ```python # View current medium print(model.medium)

# Modify medium (must reassign entire dict) medium = model.medium medium["EX_glc__D_e"] = 10.0 # Set glucose uptake medium["EX_o2_e"] = 0.0 # Anaerobic conditions model.medium = medium

# Calculate minimal media from cobra.medium import minimal_medium

# Minimize total import flux min_medium = minimal_medium(model, minimize_components=False)

# Minimize number of components (uses MILP, slower) min_medium = minimal_medium( model, minimize_components=True, open_exchanges=True ) ```

### 7. Flux Sampling

Sample the feasible flux space: ```python from cobra.sampling import sample

# Sample using OptGP (default, supports parallel processing) samples = sample(model, n=1000, method="optgp", processes=4)

# Sample using ACHR samples = sample(model, n=1000, method="achr")

# Validate samples from cobra.sampling import OptGPSampler sampler = OptGPSampler(model, processes=4) sampler.sample(1000) validation = sampler.validate(sampler.samples) print(validation.value_counts()) # Should be all 'v' for valid ```

### 8. Production Envelopes

Calculate phenotype phase planes: ```python from cobra.flux_analysis import production_envelope

# Standard production envelope envelope = production_envelope( model, reactions=["EX_glc__D_e", "EX_o2_e"], objective="EX_ac_e" # Acetate production )

# With carbon yield envelope = production_envelope( model, reactions=["EX_glc__D_e", "EX_o2_e"], carbon_sources="EX_glc__D_e" )

# Visualize (use matplotlib or pandas plotting) 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: ```python import cobra from cobra.flux_analysis import gapfill

# Build a universal model of candidate reactions to draw from. # In practice, load a curated reaction database (e.g. the BiGG # universal reactions JSON via load_json_model) or assemble reactions # into an empty cobra.Model — there is no load_model("universal"). universal = cobra.Model("universal_reactions") # ...populate `universal` with candidate cobra.Reaction objects...

with model: # Create a gap for demonstration model.remove_reactions([model.reactions.PGI])

# Find the minimal set of reactions from `universal` to restore feasibility solution = gapfill(model, universal, demand_reactions=False) for reaction in solution[0]: print(reaction.id) ```

> **Note**: `gapfill` returns a list of reaction lists (one per iteration when `iterations > 1`), so iterate `solution[0]` rather than printing `solution` directly.

### 10. Model Building

Build models from scratch: ```python from cobra import Model, Reaction, Metabolite

# Create model model = Model("my_model")

# Create metabolites 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")

# Create reaction reaction = Reaction("ATPASE") reaction.name = "ATP hydrolysis" reaction.subsystem = "Energy" reaction.lower_bound = 0.0 reaction.upper_bound = 1000.0

# Add metabolites with stoichiometry reaction.add_metabolites({ atp_c: -1.0, adp_c: 1.0, pi_c: 1.0 })

# Add gene-reaction rule reaction.gene_reaction_rule = "(gene1 and gene2) or gene3"

# Add to model model.add_reactions([reaction])

# Add boundary reactions model.add_boundary(atp_c, type="exchange") model.add_boundary(adp_c, type="demand")

# Set objective model.objective = "ATPASE" ```

## Common Workflows

### Workflow 1: Load Model and Predict Growth

```python from cobra.io import load_model

# Load model model = load_model("iJO1366")

# Run FBA solution = model.optimize() print(f"Growth rate: {solution.objective_value:.3f} /h")

# Show active pathways print(solution.fluxes[solution.fluxes.abs() > 1e-6]) ```

### Workflow 2: Gene Knockout Screen

```python from cobra.io import load_model from cobra.flux_analysis import single_gene_deletion

# Load model model = load_model("iJO1366")

# Perform single gene deletions results = single_gene_deletion(model)

# Find essential genes (growth < threshold) essential_genes = results[results["growth"] < 0.01] print(f"Found {len(essential_genes)} essential genes")

# Find genes with minimal impact neutral_genes = results[results["growth"] > 0.9 * solution.objective_value] ```

### Workflow 3: Media Optimization

```python from cobra.io import load_model from cobra.medium import minimal_medium

# Load model model = load_model("iJO1366")

# Calculate minimal medium for 50% of max growth 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

```python from cobra.io import load_model from cobra.flux_analysis import flux_variability_analysis from cobra.sampling import sample

# Load model model = load_model("iJO1366")

# First check flux ranges at optimality fva = flux_variability_analysis(model, fraction_of_optimum=1.0)

# For reactions with large ranges, sample to understand distribution samples = sample(model, n=1000)

# Analyze specific reaction 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: ```python # Model remains unchanged outside context with model: # Temporarily change objective model.objective = "ATPM"

# Temporarily modify bounds model.reactions.EX_glc__D_e.lower_bound = -5.0

# Temporarily knock out genes model.genes.b0008.knock_out()

# Optimize with changes solution = model.optimize() print(f"Modified growth: {solution.objective_value}")

# All changes automatically reverted 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: ```python # Access by index first_reaction = model.reactions[0]

# Access by ID pfk