{
  "skill": "alterlab-cobrapy",
  "evals": [
    {
      "id": "fba-growth-prediction",
      "prompt": "I have the E. coli iJO1366 genome-scale model as an SBML file. I want to run flux balance analysis to predict the growth rate and see which pathways carry flux at the optimum.",
      "expected_output": "Invokes alterlab-cobrapy. Loads the model with read_sbml_model (or load_model('iJO1366')), calls model.optimize(), reports solution.objective_value as the growth rate and solution.status, and inspects the active fluxes (e.g. solution.fluxes where abs > 1e-6).",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "optimize" },
        { "type": "behavior", "value": "Runs FBA via model.optimize() and reports the objective value as growth rate." }
      ]
    },
    {
      "id": "single-gene-deletion-essential",
      "prompt": "Using my SBML metabolic model, I want to do a single gene knockout screen to find the essential genes whose deletion abolishes growth.",
      "expected_output": "Invokes alterlab-cobrapy. Uses single_gene_deletion from cobra.flux_analysis, filters the results where growth falls below a threshold (e.g. < 0.01) to identify essential genes, and may mention the context-manager knock_out() pattern for individual genes.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "single_gene_deletion" },
        { "type": "behavior", "value": "Identifies essential genes by thresholding growth from a deletion screen." }
      ]
    },
    {
      "id": "fva-flux-variability",
      "prompt": "I want to know the feasible flux ranges for my reactions at 90% of optimal growth, and I need to eliminate thermodynamically infeasible loops from the result.",
      "expected_output": "Invokes alterlab-cobrapy. Runs flux_variability_analysis from cobra.flux_analysis with fraction_of_optimum=0.9 and loopless=True to get minimum/maximum flux per reaction without infeasible cycles.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "flux_variability_analysis" },
        { "type": "behavior", "value": "Uses fraction_of_optimum and loopless to bound flux ranges." }
      ]
    },
    {
      "id": "minimal-media-anaerobic",
      "prompt": "I need to compute the minimal medium that still supports growth for my model, and then re-run the simulation under anaerobic conditions by shutting off oxygen uptake.",
      "expected_output": "Invokes alterlab-cobrapy. Uses minimal_medium from cobra.medium to find the minimal component set, and modifies model.medium (reassigning the whole dict) to set EX_o2_e to 0 for anaerobic conditions before re-optimizing.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "minimal_medium" },
        { "type": "behavior", "value": "Computes minimal medium and edits model.medium to model anaerobic growth." }
      ]
    },
    {
      "id": "flux-sampling-distribution",
      "prompt": "I want to sample 1000 points from the feasible flux space of my model to understand the flux distribution through PFK, using parallel processing.",
      "expected_output": "Invokes alterlab-cobrapy. Uses sample from cobra.sampling (OptGP method with processes for parallelism) to draw samples, then inspects the distribution for the target reaction (e.g. a histogram of samples['PFK']).",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "sample" },
        { "type": "behavior", "value": "Uses cobra.sampling.sample (OptGP/ACHR) to characterize the flux distribution." }
      ]
    },
    {
      "id": "near-miss-pydeseq2",
      "prompt": "I have a counts matrix from a bulk RNA-seq experiment and want to find genes that are differentially expressed between knockout and wild-type samples.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-pydeseq2. The user wants statistical differential gene expression on RNA-seq counts, not constraint-based metabolic flux simulation on a genome-scale model.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-pydeseq2" }
      ]
    },
    {
      "id": "near-miss-arboreto",
      "prompt": "From my expression matrix I want to infer the gene regulatory network linking transcription factors to their target genes.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-arboreto. The user wants TF-target regulatory network inference (GRNBoost2/GENIE3) from expression data, not stoichiometric flux analysis of a metabolic network.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-arboreto" }
      ]
    },
    {
      "id": "near-miss-brenda",
      "prompt": "What are the measured kcat and Km values of E. coli phosphofructokinase for fructose-6-phosphate? I need literature kinetic constants, not a simulation.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-brenda. The user wants curated enzyme kinetic parameters (kcat, Km) from the literature, which BRENDA provides; COBRApy runs steady-state constraint-based simulations and does not look up kinetic constants.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-brenda" }
      ]
    }
  ]
}

# COBRApy API Quick Reference

This document provides quick reference for common COBRApy functions, signatures, and usage patterns.

## Model I/O

### Loading Models

```python
from cobra.io import load_model, read_sbml_model, load_json_model, load_yaml_model, load_matlab_model

# Bundled test models (only these three aliases are bundled)
model = load_model("textbook")   # E. coli core metabolism (e_coli_core)
model = load_model("iJO1366")    # Full E. coli genome-scale model
model = load_model("salmonella") # Salmonella LT2
# NOTE: any other id (e.g. "ecoli") is NOT bundled and falls through to a
# remote BiGG/BioModels lookup (needs network access).

# From files
model = read_sbml_model(filename, f_replace={}, **kwargs)
model = load_json_model(filename)
model = load_yaml_model(filename)
model = load_matlab_model(filename, variable_name=None)
```

### Saving Models

```python
from cobra.io import write_sbml_model, save_json_model, save_yaml_model, save_matlab_model

write_sbml_model(model, filename, f_replace={}, **kwargs)
save_json_model(model, filename, pretty=False, **kwargs)
save_yaml_model(model, filename, **kwargs)
save_matlab_model(model, filename, **kwargs)
```

## Model Structure

### Core Classes

```python
from cobra import Model, Reaction, Metabolite, Gene

# Create model
model = Model(id_or_model=None, name=None)

# Create metabolite
metabolite = Metabolite(
    id=None,
    formula=None,
    name="",
    charge=None,
    compartment=None
)

# Create reaction
reaction = Reaction(
    id=None,
    name="",
    subsystem="",
    lower_bound=0.0,
    upper_bound=None
)

# Create gene
gene = Gene(id=None, name="", functional=True)
```

### Model Attributes

```python
# Component access (DictList objects)
model.reactions       # DictList of Reaction objects
model.metabolites     # DictList of Metabolite objects
model.genes          # DictList of Gene objects

# Special reaction lists
model.exchanges      # Exchange reactions (external transport)
model.demands        # Demand reactions (metabolite sinks)
model.sinks          # Sink reactions
model.boundary       # All boundary reactions

# Model properties
model.objective      # Current objective (read/write)
model.objective_direction  # "max" or "min"
model.medium         # Growth medium (dict of exchange: bound)
model.solver         # Optimization solver
```

### DictList Methods

```python
# Access by index
item = model.reactions[0]

# Access by ID
item = model.reactions.get_by_id("PFK")

# Query by string (substring match)
items = model.reactions.query("atp")      # Case-insensitive search
items = model.reactions.query(lambda x: x.subsystem == "Glycolysis")

# List comprehension
items = [r for r in model.reactions if r.lower_bound < 0]

# Check membership
"PFK" in model.reactions
```

## Optimization

### Basic Optimization

```python
# Full optimization (returns Solution object)
solution = model.optimize()

# Attributes of Solution
solution.objective_value   # Objective function value
solution.status           # Optimization status ("optimal", "infeasible", etc.)
solution.fluxes          # Pandas Series of reaction fluxes
solution.shadow_prices   # Pandas Series of metabolite shadow prices
solution.reduced_costs   # Pandas Series of reduced costs

# Fast optimization (returns float only)
objective_value = model.slim_optimize()

# Change objective
model.objective = "ATPM"
model.objective = model.reactions.ATPM
model.objective = {model.reactions.ATPM: 1.0}

# Change optimization direction
model.objective_direction = "max"  # or "min"
```

### Solver Configuration

```python
# Check available solvers
from cobra.util.solver import solvers
print(solvers)

# Change solver
model.solver = "glpk"  # or "cplex", "gurobi", etc.

# Solver-specific configuration
model.solver.configuration.timeout = 60  # seconds
model.solver.configuration.verbosity = 1
model.solver.configuration.tolerances.feasibility = 1e-9
```

## Flux Analysis

### Flux Balance Analysis (FBA)

```python
from cobra.flux_analysis import pfba, geometric_fba

# Parsimonious FBA
solution = pfba(model, fraction_of_optimum=1.0, **kwargs)

# Geometric FBA
solution = geometric_fba(model, epsilon=1e-06, max_tries=200)
```

### Flux Variability Analysis (FVA)

```python
from cobra.flux_analysis import flux_variability_analysis

fva_result = flux_variability_analysis(
    model,
    reaction_list=None,        # List of reaction IDs or None for all
    loopless=False,            # Eliminate thermodynamically infeasible loops
    fraction_of_optimum=1.0,   # Optimality fraction (0.0-1.0)
    pfba_factor=None,          # Optional pFBA constraint
    processes=1                # Number of parallel processes
)

# Returns DataFrame with columns: minimum, maximum
```

### Gene and Reaction Deletions

```python
from cobra.flux_analysis import (
    single_gene_deletion,
    single_reaction_deletion,
    double_gene_deletion,
    double_reaction_deletion
)

# Single deletions
results = single_gene_deletion(
    model,
    gene_list=None,     # None for all genes
    processes=1,
    **kwargs
)

results = single_reaction_deletion(
    model,
    reaction_list=None,  # None for all reactions
    processes=1,
    **kwargs
)

# Double deletions
results = double_gene_deletion(
    model,
    gene_list1=None,
    gene_list2=None,
    processes=1,
    **kwargs
)

results = double_reaction_deletion(
    model,
    reaction_list1=None,
    reaction_list2=None,
    processes=1,
    **kwargs
)

# Returns DataFrame with columns: ids, growth, status
# For double deletions, index is MultiIndex of gene/reaction pairs
```

### Flux Sampling

```python
from cobra.sampling import sample, OptGPSampler, ACHRSampler

# Simple interface
samples = sample(
    model,
    n,                  # Number of samples
    method="optgp",     # default "auto"; also "achr" or "chrr" (needs cobra[chrr])
    thinning=100,       # Thinning factor (sample every n iterations)
    processes=1,        # Parallel processes (OptGP only)
    seed=None          # Random seed
)

# Advanced interface with sampler objects
sampler = OptGPSampler(model, processes=4, thinning=100)
sampler = ACHRSampler(model, thinning=100)

# Generate samples
samples = sampler.sample(n)

# Validate samples
validation = sampler.validate(sampler.samples)
# Returns array of 'v' (valid), 'l' (lower bound violation),
# 'u' (upper bound violation), 'e' (equality violation)

# Batch sampling
sampler.batch(n_samples, n_batches)
```

### Production Envelopes

```python
from cobra.flux_analysis import production_envelope

envelope = production_envelope(
    model,
    reactions,              # List of 1-2 reaction IDs
    objective=None,         # Objective reaction ID (None uses model objective)
    carbon_sources=None,    # Carbon source for yield calculation
    points=20,              # Number of points to calculate
    threshold=0.01          # Minimum objective value threshold
)

# Returns DataFrame with columns:
# - First reaction flux
# - Second reaction flux (if provided)
# - objective_minimum, objective_maximum
# - carbon_yield_minimum, carbon_yield_maximum (if carbon source specified)
# - mass_yield_minimum, mass_yield_maximum
```

### Gapfilling

```python
from cobra.flux_analysis import gapfill

# Basic gapfilling
solution = gapfill(
    model,
    universal=None,         # Universal model with candidate reactions
    lower_bound=0.05,       # Minimum objective flux
    penalties=None,         # Dict of reaction: penalty
    demand_reactions=True,  # Add demand reactions if needed
    exchange_reactions=False,
    iterations=1
)

# Returns a LIST OF LISTS: one inner list of cobra.Reaction objects per
# iteration. Iterate solution[0] for the first solution set.
for reaction in solution[0]:
    print(reaction.id)

# Multiple alternative solutions in one call
solutions = gapfill(model, universal, iterations=5)
for i, reaction_set in enumerate(solutions):
    print(f"Solution {i}: {[r.id for r in reaction_set]}")
```

### Other Analysis Methods

```python
from cobra.flux_analysis import (
    find_blocked_reactions,
    find_essential_genes,
    find_essential_reactions
)

# Blocked reactions (cannot carry flux)
blocked = find_blocked_reactions(
    model,
    reaction_list=None,
    zero_cutoff=1e-9,
    open_exchanges=False
)

# Essential genes/reactions
essential_genes = find_essential_genes(model, threshold=0.01)
essential_reactions = find_essential_reactions(model, threshold=0.01)
```

## Media and Boundary Conditions

### Medium Management

```python
# Get current medium (returns dict)
medium = model.medium

# Set medium (must reassign entire dict)
medium = model.medium
medium["EX_glc__D_e"] = 10.0
medium["EX_o2_e"] = 20.0
model.medium = medium

# Alternative: individual modification
with model:
    model.reactions.EX_glc__D_e.lower_bound = -10.0
```

### Minimal Media

```python
from cobra.medium import minimal_medium

# Parameter order matters if you call positionally:
min_medium = minimal_medium(
    model,
    min_objective_value=0.1,   # Minimum objective (growth) value
    exports=False,             # Allow secretion fluxes in the result
    minimize_components=False, # If True/int, minimize component COUNT via MILP (slower)
    open_exchanges=False,      # Open all exchanges before optimization
)

# Returns a pandas.Series of required import fluxes (or a DataFrame when
# minimize_components is an int requesting multiple solutions; None if infeasible)
```

### Boundary Reactions

```python
# Add boundary reaction
model.add_boundary(
    metabolite,
    type="exchange",    # or "demand", "sink"
    reaction_id=None,   # Auto-generated if None
    lb=None,
    ub=None,
    sbo_term=None
)

# Access boundary reactions
exchanges = model.exchanges     # System boundary
demands = model.demands         # Intracellular removal
sinks = model.sinks            # Intracellular exchange
boundaries = model.boundary    # All boundary reactions
```

## Model Manipulation

### Adding Components

```python
# Add reactions (always plural; pass a list even for one reaction)
model.add_reactions([reaction1, reaction2, ...])
model.add_reactions([reaction])

# Add metabolites
reaction.add_metabolites({
    metabolite1: -1.0,  # Consumed (negative stoichiometry)
    metabolite2: 1.0    # Produced (positive stoichiometry)
})

# Add metabolites to model
model.add_metabolites([metabolite1, metabolite2, ...])

# Genes are created automatically when you set a reaction's
# gene_reaction_rule, so there is rarely a need to add Gene objects
# directly. To register genes explicitly, set the GPR:
reaction.gene_reaction_rule = "gene1 and gene2"  # creates Gene objects in model.genes
```

### Removing Components

```python
# Remove reactions
model.remove_reactions([reaction1, reaction2, ...])
model.remove_reactions(["PFK", "FBA"])

# Remove metabolites (removes from reactions too)
model.remove_metabolites([metabolite1, metabolite2, ...])

# Remove genes (usually via gene_reaction_rule)
model.genes.remove(gene)
```

### Modifying Reactions

```python
# Set bounds
reaction.bounds = (lower, upper)
reaction.lower_bound = 0.0
reaction.upper_bound = 1000.0

# Modify stoichiometry
reaction.add_metabolites({metabolite: 1.0})
reaction.subtract_metabolites({metabolite: 1.0})

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

# Knock out
reaction.knock_out()
gene.knock_out()
```

### Model Copying

```python
# Deep copy (independent model)
model_copy = model.copy()

# Copy specific reactions
new_model = Model("subset")
reactions_to_copy = [model.reactions.PFK, model.reactions.FBA]
new_model.add_reactions(reactions_to_copy)
```

## Context Management

Use context managers for temporary modifications:

```python
# Changes automatically revert after with block
with model:
    model.objective = "ATPM"
    model.reactions.EX_glc__D_e.lower_bound = -5.0
    model.genes.b0008.knock_out()
    solution = model.optimize()

# Model state restored here

# Multiple nested contexts
with model:
    model.objective = "ATPM"
    with model:
        model.genes.b0008.knock_out()
        # Both modifications active
    # Only objective change active

# Context management with reactions
with model:
    model.reactions.PFK.knock_out()
    # Equivalent to: reaction.lower_bound = reaction.upper_bound = 0
```

## Reaction and Metabolite Properties

### Reaction Attributes

```python
reaction.id                      # Unique identifier
reaction.name                    # Human-readable name
reaction.subsystem               # Pathway/subsystem
reaction.bounds                  # (lower_bound, upper_bound)
reaction.lower_bound
reaction.upper_bound
reaction.reversibility          # Boolean (lower_bound < 0)
reaction.gene_reaction_rule     # GPR string
reaction.genes                  # Set of associated Gene objects
reaction.metabolites            # Dict of {metabolite: stoichiometry}

# Methods
reaction.reaction               # Stoichiometric equation string
reaction.build_reaction_string() # Same as above
reaction.check_mass_balance()   # Returns imbalances or empty dict
reaction.get_coefficient(metabolite_id)
reaction.add_metabolites({metabolite: coeff})
reaction.subtract_metabolites({metabolite: coeff})
reaction.knock_out()
```

### Metabolite Attributes

```python
metabolite.id                   # Unique identifier
metabolite.name                 # Human-readable name
metabolite.formula              # Chemical formula
metabolite.charge               # Charge
metabolite.compartment          # Compartment ID
metabolite.reactions            # FrozenSet of associated reactions

# Methods
metabolite.summary()            # Print production/consumption
metabolite.copy()
```

### Gene Attributes

```python
gene.id                         # Unique identifier
gene.name                       # Human-readable name
gene.functional                 # Boolean activity status
gene.reactions                  # FrozenSet of associated reactions

# Methods
gene.knock_out()
```

## Model Validation

### Consistency Checking

```python
from cobra.manipulation import check_mass_balance, check_metabolite_compartment_formula

# Check all reactions for mass balance
unbalanced = {}
for reaction in model.reactions:
    balance = reaction.check_mass_balance()
    if balance:
        unbalanced[reaction.id] = balance

# Check metabolite formulas are valid
check_metabolite_compartment_formula(model)
```

### Model Statistics

```python
# Basic stats
print(f"Reactions: {len(model.reactions)}")
print(f"Metabolites: {len(model.metabolites)}")
print(f"Genes: {len(model.genes)}")

# Advanced stats
print(f"Exchanges: {len(model.exchanges)}")
print(f"Demands: {len(model.demands)}")

# Blocked reactions
from cobra.flux_analysis import find_blocked_reactions
blocked = find_blocked_reactions(model)
print(f"Blocked reactions: {len(blocked)}")

# Essential genes
from cobra.flux_analysis import find_essential_genes
essential = find_essential_genes(model)
print(f"Essential genes: {len(essential)}")
```

## Summary Methods

```python
# Model summary
model.summary()                  # Overall model info

# Metabolite summary
model.metabolites.atp_c.summary()

# Reaction summary
model.reactions.PFK.summary()

# Summary with FVA
model.summary(fva=0.95)         # Include FVA at 95% optimality
```

## Common Patterns

### Batch Analysis Pattern

```python
results = []
for condition in conditions:
    with model:
        # Apply condition
        setup_condition(model, condition)

        # Analyze
        solution = model.optimize()

        # Store result
        results.append({
            "condition": condition,
            "growth": solution.objective_value,
            "status": solution.status
        })

df = pd.DataFrame(results)
```

### Systematic Knockout Pattern

```python
knockout_results = []
for gene in model.genes:
    with model:
        gene.knock_out()

        solution = model.optimize()

        knockout_results.append({
            "gene": gene.id,
            "growth": solution.objective_value if solution.status == "optimal" else 0,
            "status": solution.status
        })

df = pd.DataFrame(knockout_results)
```

### Parameter Scan Pattern

```python
parameter_values = np.linspace(0, 20, 21)
results = []

for value in parameter_values:
    with model:
        model.reactions.EX_glc__D_e.lower_bound = -value

        solution = model.optimize()

        results.append({
            "glucose_uptake": value,
            "growth": solution.objective_value,
            "acetate_secretion": solution.fluxes["EX_ac_e"]
        })

df = pd.DataFrame(results)
```

This quick reference covers the most commonly used COBRApy functions and patterns. For complete API documentation, see https://cobrapy.readthedocs.io/

# COBRApy Comprehensive Workflows

This document provides detailed step-by-step workflows for common COBRApy tasks in metabolic modeling.

## Workflow 1: Complete Knockout Study with Visualization

This workflow demonstrates how to perform a comprehensive gene knockout study and visualize the results.

```python
import pandas as pd
import matplotlib.pyplot as plt
from cobra.io import load_model
from cobra.flux_analysis import single_gene_deletion, double_gene_deletion

# Step 1: Load model
model = load_model("textbook")  # E. coli core model (bundled, runs offline)
print(f"Loaded model: {model.id}")
print(f"Model contains {len(model.reactions)} reactions, {len(model.metabolites)} metabolites, {len(model.genes)} genes")

# Step 2: Get baseline growth rate
baseline = model.slim_optimize()
print(f"Baseline growth rate: {baseline:.3f} /h")

# Step 3: Perform single gene deletions
print("Performing single gene deletions...")
single_results = single_gene_deletion(model)

# Step 4: Classify genes by impact
essential_genes = single_results[single_results["growth"] < 0.01]
severely_impaired = single_results[(single_results["growth"] >= 0.01) &
                                   (single_results["growth"] < 0.5 * baseline)]
moderately_impaired = single_results[(single_results["growth"] >= 0.5 * baseline) &
                                     (single_results["growth"] < 0.9 * baseline)]
neutral_genes = single_results[single_results["growth"] >= 0.9 * baseline]

print(f"\nSingle Deletion Results:")
print(f"  Essential genes: {len(essential_genes)}")
print(f"  Severely impaired: {len(severely_impaired)}")
print(f"  Moderately impaired: {len(moderately_impaired)}")
print(f"  Neutral genes: {len(neutral_genes)}")

# Step 5: Visualize distribution
fig, ax = plt.subplots(figsize=(10, 6))
single_results["growth"].hist(bins=50, ax=ax)
ax.axvline(baseline, color='r', linestyle='--', label='Baseline')
ax.set_xlabel("Growth rate (/h)")
ax.set_ylabel("Number of genes")
ax.set_title("Distribution of Growth Rates After Single Gene Deletions")
ax.legend()
plt.tight_layout()
plt.savefig("single_deletion_distribution.png", dpi=300)

# Step 6: Identify gene pairs for double deletions
# Focus on non-essential genes to find synthetic lethals
target_genes = single_results[single_results["growth"] >= 0.5 * baseline].index.tolist()
target_genes = [list(gene)[0] for gene in target_genes[:50]]  # Limit for performance

print(f"\nPerforming double deletions on {len(target_genes)} genes...")
double_results = double_gene_deletion(
    model,
    gene_list1=target_genes,
    processes=4
)

# Step 7: Find synthetic lethal pairs
synthetic_lethals = double_results[
    (double_results["growth"] < 0.01) &
    (single_results.loc[double_results.index.get_level_values(0)]["growth"].values >= 0.5 * baseline) &
    (single_results.loc[double_results.index.get_level_values(1)]["growth"].values >= 0.5 * baseline)
]

print(f"Found {len(synthetic_lethals)} synthetic lethal gene pairs")
print("\nTop 10 synthetic lethal pairs:")
print(synthetic_lethals.head(10))

# Step 8: Export results
single_results.to_csv("single_gene_deletions.csv")
double_results.to_csv("double_gene_deletions.csv")
synthetic_lethals.to_csv("synthetic_lethals.csv")
```

## Workflow 2: Media Design and Optimization

This workflow shows how to systematically design growth media and find minimal media compositions.

```python
from cobra.io import load_model
from cobra.medium import minimal_medium
import pandas as pd

# Step 1: Load model and check current medium
model = load_model("textbook")  # E. coli core model (bundled, runs offline)
current_medium = model.medium
print("Current medium composition:")
for exchange, bound in current_medium.items():
    metabolite_id = exchange.replace("EX_", "").replace("_e", "")
    print(f"  {metabolite_id}: {bound:.2f} mmol/gDW/h")

# Step 2: Get baseline growth
baseline_growth = model.slim_optimize()
print(f"\nBaseline growth rate: {baseline_growth:.3f} /h")

# Step 3: Calculate minimal medium for different growth targets
growth_targets = [0.25, 0.5, 0.75, 1.0]
minimal_media = {}

for fraction in growth_targets:
    target_growth = baseline_growth * fraction
    print(f"\nCalculating minimal medium for {fraction*100:.0f}% growth ({target_growth:.3f} /h)...")

    min_medium = minimal_medium(
        model,
        target_growth,
        minimize_components=True,
        open_exchanges=True
    )

    minimal_media[fraction] = min_medium
    print(f"  Required components: {len(min_medium)}")
    print(f"  Components: {list(min_medium.index)}")

# Step 4: Compare media compositions
media_df = pd.DataFrame(minimal_media).fillna(0)
media_df.to_csv("minimal_media_comparison.csv")

# Step 5: Test aerobic vs anaerobic conditions
print("\n--- Aerobic vs Anaerobic Comparison ---")

# Aerobic
model_aerobic = model.copy()
aerobic_growth = model_aerobic.slim_optimize()
aerobic_medium = minimal_medium(model_aerobic, aerobic_growth * 0.9, minimize_components=True)

# Anaerobic
model_anaerobic = model.copy()
medium_anaerobic = model_anaerobic.medium
medium_anaerobic["EX_o2_e"] = 0.0
model_anaerobic.medium = medium_anaerobic
anaerobic_growth = model_anaerobic.slim_optimize()
anaerobic_medium = minimal_medium(model_anaerobic, anaerobic_growth * 0.9, minimize_components=True)

print(f"Aerobic growth: {aerobic_growth:.3f} /h (requires {len(aerobic_medium)} components)")
print(f"Anaerobic growth: {anaerobic_growth:.3f} /h (requires {len(anaerobic_medium)} components)")

# Step 6: Identify unique requirements
aerobic_only = set(aerobic_medium.index) - set(anaerobic_medium.index)
anaerobic_only = set(anaerobic_medium.index) - set(aerobic_medium.index)
shared = set(aerobic_medium.index) & set(anaerobic_medium.index)

print(f"\nShared components: {len(shared)}")
print(f"Aerobic-only: {aerobic_only}")
print(f"Anaerobic-only: {anaerobic_only}")

# Step 7: Test custom medium
print("\n--- Testing Custom Medium ---")
custom_medium = {
    "EX_glc__D_e": 10.0,  # Glucose
    "EX_o2_e": 20.0,       # Oxygen
    "EX_nh4_e": 5.0,       # Ammonium
    "EX_pi_e": 5.0,        # Phosphate
    "EX_so4_e": 1.0,       # Sulfate
}

with model:
    model.medium = custom_medium
    custom_growth = model.optimize().objective_value
    print(f"Growth on custom medium: {custom_growth:.3f} /h")

    # Check which nutrients are limiting
    for exchange in custom_medium:
        with model:
            # Double the uptake rate
            medium_test = model.medium
            medium_test[exchange] *= 2
            model.medium = medium_test
            test_growth = model.optimize().objective_value
            improvement = (test_growth - custom_growth) / custom_growth * 100
            if improvement > 1:
                print(f"  {exchange}: +{improvement:.1f}% growth when doubled (LIMITING)")
```

## Workflow 3: Flux Space Exploration with Sampling

This workflow demonstrates comprehensive flux space analysis using FVA and sampling.

```python
from cobra.io import load_model
from cobra.flux_analysis import flux_variability_analysis
from cobra.sampling import sample
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Step 1: Load model
model = load_model("textbook")  # E. coli core model (bundled, runs offline)
baseline = model.slim_optimize()
print(f"Baseline growth: {baseline:.3f} /h")

# Step 2: Perform FVA at optimal growth
print("\nPerforming FVA at optimal growth...")
fva_optimal = flux_variability_analysis(model, fraction_of_optimum=1.0)

# Step 3: Identify reactions with flexibility
fva_optimal["range"] = fva_optimal["maximum"] - fva_optimal["minimum"]
fva_optimal["relative_range"] = fva_optimal["range"] / (fva_optimal["maximum"].abs() + 1e-9)

flexible_reactions = fva_optimal[fva_optimal["range"] > 1.0].sort_values("range", ascending=False)
print(f"\nFound {len(flexible_reactions)} reactions with >1.0 mmol/gDW/h flexibility")
print("\nTop 10 most flexible reactions:")
print(flexible_reactions.head(10)[["minimum", "maximum", "range"]])

# Step 4: Perform FVA at suboptimal growth (90%)
print("\nPerforming FVA at 90% optimal growth...")
fva_suboptimal = flux_variability_analysis(model, fraction_of_optimum=0.9)
fva_suboptimal["range"] = fva_suboptimal["maximum"] - fva_suboptimal["minimum"]

# Step 5: Compare flexibility at different optimality levels
comparison = pd.DataFrame({
    "range_100": fva_optimal["range"],
    "range_90": fva_suboptimal["range"]
})
comparison["range_increase"] = comparison["range_90"] - comparison["range_100"]

print("\nReactions with largest increase in flexibility at suboptimality:")
print(comparison.sort_values("range_increase", ascending=False).head(10))

# Step 6: Perform flux sampling
print("\nPerforming flux sampling (1000 samples)...")
samples = sample(model, n=1000, method="optgp", processes=4)  # or leave method="auto"

# Step 7: Analyze sampling results for key reactions
key_reactions = ["PFK", "FBA", "TPI", "GAPD", "PGK", "PGM", "ENO", "PYK"]
available_key_reactions = [r for r in key_reactions if r in samples.columns]

if available_key_reactions:
    fig, axes = plt.subplots(2, 4, figsize=(16, 8))
    axes = axes.flatten()

    for idx, reaction_id in enumerate(available_key_reactions[:8]):
        ax = axes[idx]
        samples[reaction_id].hist(bins=30, ax=ax, alpha=0.7)

        # Overlay FVA bounds
        fva_min = fva_optimal.loc[reaction_id, "minimum"]
        fva_max = fva_optimal.loc[reaction_id, "maximum"]
        ax.axvline(fva_min, color='r', linestyle='--', label='FVA min')
        ax.axvline(fva_max, color='r', linestyle='--', label='FVA max')

        ax.set_xlabel("Flux (mmol/gDW/h)")
        ax.set_ylabel("Frequency")
        ax.set_title(reaction_id)
        if idx == 0:
            ax.legend()

    plt.tight_layout()
    plt.savefig("flux_distributions.png", dpi=300)

# Step 8: Calculate correlation between reactions
print("\nCalculating flux correlations...")
correlation_matrix = samples[available_key_reactions].corr()

fig, ax = plt.subplots(figsize=(10, 8))
sns.heatmap(correlation_matrix, annot=True, fmt=".2f", cmap="coolwarm",
            center=0, ax=ax, square=True)
ax.set_title("Flux Correlations Between Key Glycolysis Reactions")
plt.tight_layout()
plt.savefig("flux_correlations.png", dpi=300)

# Step 9: Identify reaction modules (highly correlated groups)
print("\nHighly correlated reaction pairs (|r| > 0.9):")
for i in range(len(correlation_matrix)):
    for j in range(i+1, len(correlation_matrix)):
        corr = correlation_matrix.iloc[i, j]
        if abs(corr) > 0.9:
            print(f"  {correlation_matrix.index[i]} <-> {correlation_matrix.columns[j]}: {corr:.3f}")

# Step 10: Export all results
fva_optimal.to_csv("fva_optimal.csv")
fva_suboptimal.to_csv("fva_suboptimal.csv")
samples.to_csv("flux_samples.csv")
correlation_matrix.to_csv("flux_correlations.csv")
```

## Workflow 4: Production Strain Design

This workflow demonstrates how to design a production strain for a target metabolite.

```python
from cobra.io import load_model
from cobra.flux_analysis import (
    production_envelope,
    flux_variability_analysis,
    single_gene_deletion
)
import pandas as pd
import matplotlib.pyplot as plt

# Step 1: Define production target
TARGET_METABOLITE = "EX_ac_e"  # Acetate production
CARBON_SOURCE = "EX_glc__D_e"  # Glucose uptake

# Step 2: Load model
model = load_model("textbook")  # E. coli core model (bundled, runs offline)
print(f"Designing strain for {TARGET_METABOLITE} production")

# Step 3: Calculate baseline production envelope
print("\nCalculating production envelope...")
envelope = production_envelope(
    model,
    reactions=[CARBON_SOURCE, TARGET_METABOLITE],
    carbon_sources=CARBON_SOURCE
)

# Visualize production envelope
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(envelope[CARBON_SOURCE], envelope["mass_yield_maximum"], 'b-', label='Max yield')
ax.plot(envelope[CARBON_SOURCE], envelope["mass_yield_minimum"], 'r-', label='Min yield')
ax.set_xlabel(f"Glucose uptake (mmol/gDW/h)")
ax.set_ylabel(f"Acetate yield")
ax.set_title("Wild-type Production Envelope")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("production_envelope_wildtype.png", dpi=300)

# Step 4: Maximize production while maintaining growth
print("\nOptimizing for production...")

BIOMASS = "BIOMASS_Ecoli_core_w_GAM"  # biomass reaction id in the textbook model

# Record max growth before re-tasking the objective
max_growth = model.slim_optimize()
MIN_GROWTH = 0.1 * max_growth  # require at least 10% of max growth

# Permanently constrain growth and switch the objective to product formation.
# (These are deliberate model changes that persist through the knockout screen
# below; use `with model:` instead if you want them to revert.)
model.reactions.get_by_id(BIOMASS).lower_bound = MIN_GROWTH
model.objective = TARGET_METABOLITE
model.objective_direction = "max"

production_solution = model.optimize()
max_production = production_solution.objective_value
print(f"Maximum production: {max_production:.3f} mmol/gDW/h")
print(f"Growth rate: {production_solution.fluxes[BIOMASS]:.3f} /h")

# Step 5: Identify beneficial gene knockouts
print("\nScreening for beneficial knockouts...")
# The growth constraint and product objective from Step 4 are still in effect.

knockout_results = []
for gene in model.genes:
    with model:
        gene.knock_out()
        try:
            solution = model.optimize()
            if solution.status == "optimal":
                production = solution.objective_value
                growth = solution.fluxes["BIOMASS_Ecoli_core_w_GAM"]

                if production > max_production * 1.05:  # >5% improvement
                    knockout_results.append({
                        "gene": gene.id,
                        "production": production,
                        "growth": growth,
                        "improvement": (production / max_production - 1) * 100
                    })
        except:
            continue

knockout_df = pd.DataFrame(knockout_results)
if len(knockout_df) > 0:
    knockout_df = knockout_df.sort_values("improvement", ascending=False)
    print(f"\nFound {len(knockout_df)} beneficial knockouts:")
    print(knockout_df.head(10))
    knockout_df.to_csv("beneficial_knockouts.csv", index=False)
else:
    print("No beneficial single knockouts found")

# Step 6: Test combination of best knockouts
if len(knockout_df) > 0:
    print("\nTesting knockout combinations...")
    top_genes = knockout_df.head(3)["gene"].tolist()

    with model:
        for gene_id in top_genes:
            model.genes.get_by_id(gene_id).knock_out()

        solution = model.optimize()
        if solution.status == "optimal":
            combined_production = solution.objective_value
            combined_growth = solution.fluxes["BIOMASS_Ecoli_core_w_GAM"]
            combined_improvement = (combined_production / max_production - 1) * 100

            print(f"\nCombined knockout results:")
            print(f"  Genes: {', '.join(top_genes)}")
            print(f"  Production: {combined_production:.3f} mmol/gDW/h")
            print(f"  Growth: {combined_growth:.3f} /h")
            print(f"  Improvement: {combined_improvement:.1f}%")

# Step 7: Analyze flux distribution in production strain
if len(knockout_df) > 0:
    best_gene = knockout_df.iloc[0]["gene"]

    with model:
        model.genes.get_by_id(best_gene).knock_out()
        solution = model.optimize()

        # Get active pathways
        active_fluxes = solution.fluxes[solution.fluxes.abs() > 0.1]
        active_fluxes.to_csv(f"production_strain_fluxes_{best_gene}_knockout.csv")

        print(f"\nActive reactions in production strain: {len(active_fluxes)}")
```

## Workflow 5: Model Validation and Debugging

This workflow shows systematic approaches to validate and debug metabolic models.

```python
from cobra.io import load_model, read_sbml_model
from cobra.flux_analysis import flux_variability_analysis
import pandas as pd

# Step 1: Load model
model = load_model("textbook")  # E. coli core; or read_sbml_model("your_model.xml")
print(f"Model: {model.id}")
print(f"Reactions: {len(model.reactions)}")
print(f"Metabolites: {len(model.metabolites)}")
print(f"Genes: {len(model.genes)}")

# Step 2: Check model feasibility
print("\n--- Feasibility Check ---")
try:
    objective_value = model.slim_optimize()
    print(f"Model is feasible (objective: {objective_value:.3f})")
except:
    print("Model is INFEASIBLE")
    print("Troubleshooting steps:")

    # Check for blocked reactions
    from cobra.flux_analysis import find_blocked_reactions
    blocked = find_blocked_reactions(model)
    print(f"  Blocked reactions: {len(blocked)}")
    if len(blocked) > 0:
        print(f"  First 10 blocked: {list(blocked)[:10]}")

    # Check medium
    print(f"\n  Current medium: {model.medium}")

    # Try opening all exchanges
    for reaction in model.exchanges:
        reaction.lower_bound = -1000

    try:
        objective_value = model.slim_optimize()
        print(f"\n  Model feasible with open exchanges (objective: {objective_value:.3f})")
        print("  Issue: Medium constraints too restrictive")
    except:
        print("\n  Model still infeasible with open exchanges")
        print("  Issue: Structural problem (missing reactions, mass imbalance, etc.)")

# Step 3: Check mass and charge balance
print("\n--- Mass and Charge Balance Check ---")
unbalanced_reactions = []
for reaction in model.reactions:
    try:
        balance = reaction.check_mass_balance()
        if balance:
            unbalanced_reactions.append({
                "reaction": reaction.id,
                "imbalance": balance
            })
    except:
        pass

if unbalanced_reactions:
    print(f"Found {len(unbalanced_reactions)} unbalanced reactions:")
    for item in unbalanced_reactions[:10]:
        print(f"  {item['reaction']}: {item['imbalance']}")
else:
    print("All reactions are mass balanced")

# Step 4: Identify dead-end metabolites
print("\n--- Dead-end Metabolite Check ---")
dead_end_metabolites = []
for metabolite in model.metabolites:
    producing_reactions = [r for r in metabolite.reactions
                          if r.metabolites[metabolite] > 0]
    consuming_reactions = [r for r in metabolite.reactions
                          if r.metabolites[metabolite] < 0]

    if len(producing_reactions) == 0 or len(consuming_reactions) == 0:
        dead_end_metabolites.append({
            "metabolite": metabolite.id,
            "producers": len(producing_reactions),
            "consumers": len(consuming_reactions)
        })

if dead_end_metabolites:
    print(f"Found {len(dead_end_metabolites)} dead-end metabolites:")
    for item in dead_end_metabolites[:10]:
        print(f"  {item['metabolite']}: {item['producers']} producers, {item['consumers']} consumers")
else:
    print("No dead-end metabolites found")

# Step 5: Check for duplicate reactions
print("\n--- Duplicate Reaction Check ---")
reaction_equations = {}
duplicates = []

for reaction in model.reactions:
    equation = reaction.build_reaction_string()
    if equation in reaction_equations:
        duplicates.append({
            "reaction1": reaction_equations[equation],
            "reaction2": reaction.id,
            "equation": equation
        })
    else:
        reaction_equations[equation] = reaction.id

if duplicates:
    print(f"Found {len(duplicates)} duplicate reaction pairs:")
    for item in duplicates[:10]:
        print(f"  {item['reaction1']} == {item['reaction2']}")
else:
    print("No duplicate reactions found")

# Step 6: Identify orphan genes
print("\n--- Orphan Gene Check ---")
orphan_genes = [gene for gene in model.genes if len(gene.reactions) == 0]

if orphan_genes:
    print(f"Found {len(orphan_genes)} orphan genes (not associated with reactions):")
    print(f"  First 10: {[g.id for g in orphan_genes[:10]]}")
else:
    print("No orphan genes found")

# Step 7: Check for thermodynamically infeasible loops
print("\n--- Thermodynamic Loop Check ---")
fva_loopless = flux_variability_analysis(model, loopless=True)
fva_standard = flux_variability_analysis(model)

loop_reactions = []
for reaction_id in fva_standard.index:
    standard_range = fva_standard.loc[reaction_id, "maximum"] - fva_standard.loc[reaction_id, "minimum"]
    loopless_range = fva_loopless.loc[reaction_id, "maximum"] - fva_loopless.loc[reaction_id, "minimum"]

    if standard_range > loopless_range + 0.1:
        loop_reactions.append({
            "reaction": reaction_id,
            "standard_range": standard_range,
            "loopless_range": loopless_range
        })

if loop_reactions:
    print(f"Found {len(loop_reactions)} reactions potentially involved in loops:")
    loop_df = pd.DataFrame(loop_reactions).sort_values("standard_range", ascending=False)
    print(loop_df.head(10))
else:
    print("No thermodynamically infeasible loops detected")

# Step 8: Generate validation report
print("\n--- Generating Validation Report ---")
validation_report = {
    "model_id": model.id,
    "feasible": objective_value if 'objective_value' in locals() else None,
    "n_reactions": len(model.reactions),
    "n_metabolites": len(model.metabolites),
    "n_genes": len(model.genes),
    "n_unbalanced": len(unbalanced_reactions),
    "n_dead_ends": len(dead_end_metabolites),
    "n_duplicates": len(duplicates),
    "n_orphan_genes": len(orphan_genes),
    "n_loop_reactions": len(loop_reactions)
}

validation_df = pd.DataFrame([validation_report])
validation_df.to_csv("model_validation_report.csv", index=False)
print("Validation report saved to model_validation_report.csv")
```

These workflows provide comprehensive templates for common COBRApy tasks. Adapt them as needed for specific research questions and models.