# Metabolic Vulnerability Index (MVI)

**Description:** Compute a composite FBA-based gene essentiality score to prioritize antibiotic drug targets in *E. coli* and *M. tuberculosis*. The MVI combines three orthogonal signals derived from genome-scale metabolic models to rank genes by their vulnerability as drug target candidates.

**Formula:** `MVI(g) = 0.5·ΔGrowth(g) + 0.3·FluxCentrality(g) + 0.2·PathwayEssentiality(g)`

Where:
- `ΔGrowth(g)` = normalised growth reduction on gene knockout (0–1)
- `FluxCentrality(g)` = fraction of total absolute flux carried by gene's reactions
- `PathwayEssentiality(g)` = fraction of gene's reactions whose individual deletion is lethal

**Models used:**
- *E. coli* iJO1366: 1,367 genes, 2,583 reactions, wildtype growth = 0.9824 h⁻¹
- *M. tuberculosis* iEK1008: 1,008 genes, 1,226 reactions, wildtype growth = 0.0582 h⁻¹

---

## Dependencies

```
pip install cobra==0.29.0 numpy==1.26.4 scikit-learn==1.4.2 scipy==1.13.0
```

---

## Setup

Create directory structure:

```bash
mkdir -p scripts data/models artifacts
```

---

## Step 1: Write gold-standard data files

Write to `data/keio_essential_genes.csv`:

```
gene_id,gene_name,source
b0015,dnaJ,Keio2006
b0016,dnaK,Keio2006
b0025,ycaC,Keio2006
b0052,cysB,Keio2006
b0057,araC,Keio2006
b0063,leuD,Keio2006
b0064,leuC,Keio2006
b0071,leuB,Keio2006
b0072,leuA,Keio2006
b0116,accD,Keio2006
b0118,pdhR,Keio2006
b0119,aceE,Keio2006
b0120,aceF,Keio2006
b0121,lpd,Keio2006
b0185,folA,Keio2006
b0186,apaH,Keio2006
b0215,ftsL,Keio2006
b0252,murB,Keio2006
b0254,birA,Keio2006
b0317,ispH,Keio2006
b0318,lptD,Keio2006
b0384,ftsB,Keio2006
b0395,mrdB,Keio2006
b0396,mrdA,Keio2006
b0399,nadD,Keio2006
b0453,coaE,Keio2006
b0474,holB,Keio2006
b0485,lpxC,Keio2006
b0535,murJ,Keio2006
b0544,ribF,Keio2006
b0555,ispA,Keio2006
b0590,ispB,Keio2006
b0618,accA,Keio2006
b0619,accB,Keio2006
b0680,coaA,Keio2006
b0726,fldA,Keio2006
b0727,fpr,Keio2006
b0728,ispG,Keio2006
b0756,dapB,Keio2006
b0758,folC,Keio2006
b0773,hemA,Keio2006
b0812,can,Keio2006
b0870,fabI,Keio2006
b0907,mraY,Keio2006
b0908,murD,Keio2006
b0909,murG,Keio2006
b0910,murC,Keio2006
b0911,ddlB,Keio2006
b0912,ftsQ,Keio2006
b0913,ftsA,Keio2006
b0914,ftsZ,Keio2006
b0928,ribB,Keio2006
b0929,ribA,Keio2006
b0931,coaBC,Keio2006
b1000,yceI,Keio2006
b1001,rpsT,Keio2006
b1002,era,Keio2006
b1003,recO,Keio2006
b1004,rpsB,Keio2006
b1005,tsf,Keio2006
b1006,pyrH,Keio2006
b1007,frr,Keio2006
b1080,murI,Keio2006
b1100,lptE,Keio2006
b1101,lptA,Keio2006
b1130,ftsH,Keio2006
b1161,suhB,Keio2006
b1208,eno,Keio2006
b1209,pyrG,Keio2006
b1210,spoT,Keio2006
b1261,lpxA,Keio2006
b1262,lpxB,Keio2006
b1263,lpxK,Keio2006
b1264,rnhB,Keio2006
b1265,kdtA,Keio2006
b1278,accC,Keio2006
b1389,rplS,Keio2006
b1390,trmD,Keio2006
b1391,rimM,Keio2006
b1392,rpsP,Keio2006
b1393,ffh,Keio2006
b1397,ispE,Keio2006
b1440,rpsA,Keio2006
b1474,kdsB,Keio2006
b1488,lpxH,Keio2006
b1548,rpoA,Keio2006
b1549,rplQ,Keio2006
b1552,rpsD,Keio2006
b1553,rpsK,Keio2006
b1554,rpoB,Keio2006
b1555,rpoC,Keio2006
b1583,ispD,Keio2006
b1584,ispF,Keio2006
b1643,lptB,Keio2006
b1660,rpsO,Keio2006
b1661,pnp,Keio2006
b1716,ffh,Keio2006
b1718,ftsY,Keio2006
b1719,ispU,Keio2006
b1720,clsA,Keio2006
b1723,hemD,Keio2006
b1724,hemC,Keio2006
b1766,murE,Keio2006
b1767,murF,Keio2006
b1768,mraW,Keio2006
b1787,dxs,Keio2006
b1819,nadE,Keio2006
b1820,holA,Keio2006
b1838,fabD,Keio2006
b1839,fabG,Keio2006
b1840,acpP,Keio2006
b1841,fabF,Keio2006
b1900,hemB,Keio2006
b1901,lolA,Keio2006
b1902,msbA,Keio2006
b1917,rho,Keio2006
b2013,yeiA,Keio2006
b2024,rplU,Keio2006
b2025,rpmA,Keio2006
b2026,rplB,Keio2006
b2027,rplW,Keio2006
b2028,rplD,Keio2006
b2029,rplC,Keio2006
b2030,rpsJ,Keio2006
b2033,fusA,Keio2006
b2034,rpsG,Keio2006
b2035,rpsL,Keio2006
b2036,rplJ,Keio2006
b2037,rplA,Keio2006
b2039,rplK,Keio2006
b2040,rplY,Keio2006
b2041,rpsI,Keio2006
b2042,rplM,Keio2006
b2043,rpsQ,Keio2006
b2044,rpmD,Keio2006
b2045,rplX,Keio2006
b2046,rplE,Keio2006
b2047,rpsN,Keio2006
b2048,rpsH,Keio2006
b2049,rplF,Keio2006
b2050,rplR,Keio2006
b2051,rpsE,Keio2006
b2052,rpmC,Keio2006
b2053,rplO,Keio2006
b2054,rpsM,Keio2006
b2055,rpsK,Keio2006
b2056,rpoA,Keio2006
b2058,rplL,Keio2006
b2084,gyrB,Keio2006
b2085,gyrA,Keio2006
b2235,glmU,Keio2006
b2303,dapD,Keio2006
b2304,map,Keio2006
b2305,rpsB,Keio2006
b2392,glmM,Keio2006
b2393,ftsW,Keio2006
b2409,nadC,Keio2006
b2472,hisS,Keio2006
b2530,rplI,Keio2006
b2606,glmS,Keio2006
b2607,pgsA,Keio2006
b2608,yibB,Keio2006
b2655,lolB,Keio2006
b2699,rpsU,Keio2006
b2700,dnaG,Keio2006
b2701,rpoD,Keio2006
b2703,ftsE,Keio2006
b2704,ftsX,Keio2006
b2722,bamA,Keio2006
b2741,pgk,Keio2006
b2742,fba,Keio2006
b2743,tpiA,Keio2006
b2755,grpE,Keio2006
b2784,rpsF,Keio2006
b2785,priB,Keio2006
b2786,rpsR,Keio2006
b2787,rplI,Keio2006
b2819,pyrI,Keio2006
b2820,pyrB,Keio2006
b2821,pyrL,Keio2006
b2874,dnaA,Keio2006
b2875,dnaN,Keio2006
b2903,hemE,Keio2006
b3008,fabH,Keio2006
b3038,bamD,Keio2006
b3071,pssA,Keio2006
b3096,rplT,Keio2006
b3097,rpmI,Keio2006
b3165,rpsC,Keio2006
b3166,rplP,Keio2006
b3167,rpmC,Keio2006
b3168,rpsS,Keio2006
b3169,rplV,Keio2006
b3170,rpsC,Keio2006
b3175,rpsQ,Keio2006
b3176,rplN,Keio2006
b3177,rplX,Keio2006
b3178,rplE,Keio2006
b3195,coaD,Keio2006
b3196,rpmB,Keio2006
b3197,rpmG,Keio2006
b3231,asd,Keio2006
b3261,dapA,Keio2006
b3308,ileS,Keio2006
b3309,lolD,Keio2006
b3310,lolC,Keio2006
b3313,secY,Keio2006
b3314,rpmJ,Keio2006
b3315,rpsI,Keio2006
b3316,rpmJ,Keio2006
b3317,rplM,Keio2006
b3339,dapE,Keio2006
b3389,murA,Keio2006
b3390,yhdE,Keio2006
b3500,ispH,Keio2006
b3553,dnaX,Keio2006
b3554,holE,Keio2006
b3639,hemH,Keio2006
b3640,hemG,Keio2006
b3693,valS,Keio2006
b3720,psd,Keio2006
b3731,fabB,Keio2006
b3732,plsB,Keio2006
b3769,dnaB,Keio2006
b3784,nusA,Keio2006
b3785,infB,Keio2006
b3786,rbfA,Keio2006
b3787,truB,Keio2006
b3788,rpsO,Keio2006
b3789,pnp,Keio2006
b3813,lptC,Keio2006
b3865,trxB,Keio2006
b3931,rpsB,Keio2006
b3932,rpsL,Keio2006
b3953,mukF,Keio2006
b3954,mukE,Keio2006
b3955,mukB,Keio2006
b3960,rne,Keio2006
b3996,rplT,Keio2006
b4009,pdxJ,Keio2006
b4010,ispH,Keio2006
b4054,leuS,Keio2006
b4105,yrfF,Keio2006
b4143,dnaE,Keio2006
b4161,argS,Keio2006
b4172,rpmH,Keio2006
b4173,rnpA,Keio2006
b4174,yjgA,Keio2006
b4175,rpsU,Keio2006
b4176,dnaG,Keio2006
b4177,rpoD,Keio2006
b4232,fabA,Keio2006
b4238,pheS,Keio2006
b4239,pheT,Keio2006
b4274,csdA,Keio2006
b4275,iscU,Keio2006
b4276,iscA,Keio2006
b4277,hscB,Keio2006
b4278,hscA,Keio2006
b4279,fdx,Keio2006
b4282,metK,Keio2006
b4319,dxr,Keio2006
b4321,ispC,Keio2006
b4326,rplB,Keio2006
b4327,rplW,Keio2006
b4328,rplD,Keio2006
b4329,rplC,Keio2006
b4330,rpsJ,Keio2006
b4374,bamB,Keio2006
b4375,ssb,Keio2006
b4376,ubiA,Keio2006
b4382,ligA,Keio2006
b4390,lpxD,Keio2006
b4391,fabZ,Keio2006
b4392,lpxA,Keio2006
b4393,lpxB,Keio2006
b4394,murB,Keio2006
b4395,ubiB,Keio2006
b4399,nrdA,Keio2006
b4400,nrdB,Keio2006
b4401,fldB,Keio2006
```

Write to `data/tb_drug_targets.csv`:

```
gene_id,gene_name,drug,source
Rv1484,inhA,Isoniazid,DrugBank
Rv2243,fabD,Thiolactomycin,DrugBank
Rv2245,kasB,Thiolactomycin,DrugBank
Rv2246,kasA,Thiolactomycin,DrugBank
Rv0447c,murA,Fosfomycin,DrugBank
Rv2152c,murB,Lipo-murB_inhibitors,DrugBank
Rv0573c,dxr,Fosmidomycin,DrugBank
Rv3280,accD5,Soraphen,DrugBank
Rv3281,accD4,Soraphen,DrugBank
Rv2244,acpM,Thiolactomycin,DrugBank
Rv0635,acpS,Pantetheine_analogs,DrugBank
Rv3594,gyrA,Fluoroquinolones,WHO
Rv3593,gyrB,Fluoroquinolones,WHO
Rv0667,rpoB,Rifampicin,WHO
Rv0668,rpoC,Rifampicin,WHO
Rv3457c,rpoA,Rifampicin,WHO
Rv1304,rpoZ,Rifampicin,WHO
Rv0006,gyrA,Moxifloxacin,WHO
Rv0007,gyrB,Moxifloxacin,WHO
Rv1908c,katG,Isoniazid_activation,WHO
Rv3546,panC,Pantothenate_analogs,DrugBank
Rv3548c,panB,Pantothenate_analogs,DrugBank
Rv3593,coaBC,Pantothenate_analogs,DrugBank
Rv1416,coaD,Pantothenate_analogs,DrugBank
Rv2977c,coaE,Pantothenate_analogs,DrugBank
Rv0555,pncA,Pyrazinamide_activation,WHO
Rv1694,tlyA,Capreomycin,WHO
Rv1258c,tap,Ethambutol_efflux,WHO
Rv3795,embB,Ethambutol,WHO
Rv3793,embC,Ethambutol,WHO
Rv3794,embA,Ethambutol,WHO
Rv3423c,rpsL,Streptomycin,WHO
Rv0682,rpsL,Streptomycin,WHO
Rv3919c,gid,Streptomycin,WHO
```

---

## Step 2: Write pipeline scripts

Write to `scripts/download_models.py`:

```python
#!/usr/bin/env python3
"""Download BiGG metabolic models for MVI pipeline.

Usage:
    python scripts/download_models.py --organisms ecoli mtb --outdir data/models/
"""
import argparse
import json
import os
import urllib.request

BIGG_BASE = "http://bigg.ucsd.edu/static/models"

ORGANISM_MAP = {
    "ecoli": "iJO1366",
    "mtb": "iEK1008",
}


def download_model(model_id: str, outdir: str) -> str:
    """Download model JSON from BiGG; skip if already exists. Returns path."""
    os.makedirs(outdir, exist_ok=True)
    outpath = os.path.join(outdir, f"{model_id}.json")
    if os.path.exists(outpath):
        print(f"[skip] {outpath} already exists")
        return outpath
    url = f"{BIGG_BASE}/{model_id}.json"
    print(f"[download] {url} -> {outpath}")
    urllib.request.urlretrieve(url, outpath)
    # Validate JSON
    with open(outpath) as f:
        data = json.load(f)
    assert "reactions" in data, f"Invalid BiGG JSON: missing 'reactions' key in {outpath}"
    print(f"[ok] {model_id}: {len(data['reactions'])} reactions, {len(data['genes'])} genes")
    return outpath


def main():
    parser = argparse.ArgumentParser(description="Download BiGG metabolic models")
    parser.add_argument("--organisms", nargs="+", default=["ecoli", "mtb"],
                        choices=list(ORGANISM_MAP.keys()),
                        help="Organisms to download")
    parser.add_argument("--outdir", default="data/models",
                        help="Output directory for model JSON files")
    args = parser.parse_args()

    for org in args.organisms:
        model_id = ORGANISM_MAP[org]
        download_model(model_id, args.outdir)

    print("[done] All models ready.")


if __name__ == "__main__":
    main()
```

Write to `scripts/run_fba_baseline.py`:

```python
#!/usr/bin/env python3
"""Run FBA baseline: wildtype growth rate for a metabolic model."""
import argparse
import json
import os
import cobra
import cobra.io


def run_baseline(model_path: str, output_path: str) -> None:
    if not os.path.exists(model_path):
        raise FileNotFoundError(f"Model file not found: {model_path}")

    model = cobra.io.load_json_model(model_path)
    solution = model.optimize()

    if solution.status == "optimal":
        flux_dist = {rxn_id: float(flux) for rxn_id, flux in solution.fluxes.items()}
        growth_rate = float(solution.objective_value)
    else:
        print(f"WARNING: FBA status is '{solution.status}' — growth rate set to 0.0", flush=True)
        flux_dist = {}
        growth_rate = 0.0

    # Total absolute flux (sum of |flux| across all reactions)
    total_absolute_flux = float(sum(abs(v) for v in flux_dist.values()))

    # Model ID from filename
    model_id = os.path.splitext(os.path.basename(model_path))[0]

    result = {
        "model_id": model_id,
        "model_path": model_path,
        "growth_rate": growth_rate,
        "total_absolute_flux": total_absolute_flux,
        "n_genes": len(model.genes),
        "n_reactions": len(model.reactions),
        "flux_distribution": flux_dist,
        "solver": model.solver.name,
        "status": solution.status,
    }

    out_dir = os.path.dirname(output_path)
    if out_dir:
        os.makedirs(out_dir, exist_ok=True)
    with open(output_path, "w") as f:
        json.dump(result, f, indent=2)

    print(f"Model: {model_id}")
    print(f"Growth rate: {growth_rate:.4f}")
    print(f"Genes: {result['n_genes']}, Reactions: {result['n_reactions']}")
    print(f"Saved to: {output_path}")


def main():
    parser = argparse.ArgumentParser(description="Run FBA baseline for a metabolic model")
    parser.add_argument("--model", required=True, help="Path to BiGG model JSON")
    parser.add_argument("--output", required=True, help="Output JSON path")
    args = parser.parse_args()
    run_baseline(args.model, args.output)


if __name__ == "__main__":
    main()
```

Write to `scripts/run_knockout_screen.py`:

```python
#!/usr/bin/env python3
"""Run single-gene knockout screen for a metabolic model."""
import argparse
import json
import os
import cobra
import cobra.io
import cobra.flux_analysis


def run_knockout_screen(model_path: str, baseline_path: str, output_path: str) -> None:
    if not os.path.exists(model_path):
        raise FileNotFoundError(f"Model file not found: {model_path}")
    if not os.path.exists(baseline_path):
        raise FileNotFoundError(f"Baseline file not found: {baseline_path}")

    with open(baseline_path) as f:
        baseline = json.load(f)
    wt_growth = baseline["growth_rate"]

    model = cobra.io.load_json_model(model_path)
    model_id = os.path.splitext(os.path.basename(model_path))[0]

    print(f"Running single-gene knockout screen for {model_id} ({len(model.genes)} genes)...")
    deletion_results = cobra.flux_analysis.single_gene_deletion(model)

    # Build lookup: gene_id -> ko_growth
    # cobra may return integer index with 'ids' column containing frozensets,
    # or frozenset index directly depending on version.
    if "ids" in deletion_results.columns:
        # Integer-indexed with 'ids' column of frozensets
        ko_map = {}
        for _, row in deletion_results.iterrows():
            ids = row["ids"]
            growth = row["growth"]
            if not (growth == growth):  # NaN
                growth = 0.0
            growth = max(0.0, float(growth))
            for gid in ids:
                ko_map[gid] = growth
    else:
        # frozenset-indexed
        ko_map = {}
        for idx, row in deletion_results.iterrows():
            growth = row["growth"]
            if not (growth == growth):
                growth = 0.0
            growth = max(0.0, float(growth))
            for gid in idx:
                ko_map[gid] = growth

    genes = []
    for gene in model.genes:
        gene_id = gene.id
        ko_growth = ko_map.get(gene_id, 0.0)
        if ko_growth != ko_growth:  # NaN check
            ko_growth = 0.0

        if wt_growth > 0:
            delta_growth = max(0.0, min(1.0, 1.0 - (ko_growth / wt_growth)))
        else:
            delta_growth = 0.0

        status = "essential" if ko_growth < 0.01 * wt_growth else "nonessential"
        genes.append({
            "gene_id": gene_id,
            "ko_growth": ko_growth,
            "delta_growth": delta_growth,
            "status": status,
        })

    n_essential = sum(1 for g in genes if g["status"] == "essential")
    result = {
        "model_id": model_id,
        "wt_growth": wt_growth,
        "n_genes_screened": len(genes),
        "n_essential": n_essential,
        "genes": genes,
    }

    out_dir = os.path.dirname(output_path)
    if out_dir:
        os.makedirs(out_dir, exist_ok=True)
    with open(output_path, "w") as f:
        json.dump(result, f, indent=2)

    print(f"Screened: {len(genes)} genes, Essential: {n_essential}")
    print(f"Saved to: {output_path}")


def main():
    parser = argparse.ArgumentParser(description="Run single-gene knockout screen")
    parser.add_argument("--model", required=True, help="Path to BiGG model JSON")
    parser.add_argument("--baseline", required=True, help="Path to baseline JSON")
    parser.add_argument("--output", required=True, help="Output JSON path")
    args = parser.parse_args()
    run_knockout_screen(args.model, args.baseline, args.output)


if __name__ == "__main__":
    main()
```

Write to `scripts/compute_mvi.py`:

```python
#!/usr/bin/env python3
"""Compute Metabolic Vulnerability Index (MVI) for each gene."""
import argparse
import json
import math
import os
import cobra
import cobra.io
import cobra.flux_analysis


def compute_flux_centrality(model, flux_distribution: dict, total_absolute_flux: float) -> dict:
    """FluxCentrality(g) = sum(|flux[r]| for r in gene.reactions) / total_absolute_flux"""
    centrality = {}
    for gene in model.genes:
        if total_absolute_flux > 0:
            gene_flux = sum(
                abs(flux_distribution.get(rxn.id, 0.0))
                for rxn in gene.reactions
            )
            centrality[gene.id] = gene_flux / total_absolute_flux
        else:
            centrality[gene.id] = 0.0
    return centrality


def compute_pathway_essentiality(model, wt_growth: float) -> dict:
    """PathwayEssentiality(g) = fraction of gene's reactions where single-reaction KO kills growth."""
    print("Computing pathway essentiality (single-reaction deletions)...")
    rxn_deletion_results = cobra.flux_analysis.single_reaction_deletion(model)

    # Build reaction essentiality map
    rxn_essential = {}
    threshold = 0.01 * wt_growth

    def _get_growth(row, df_columns):
        if "growth" in df_columns:
            return row["growth"]
        # fallback: first non-id, non-status column
        for col in df_columns:
            if col not in ("ids", "status"):
                return row[col]
        return 0.0

    # Handle both old (frozenset-indexed) and new (ids column) cobra formats
    if "ids" in rxn_deletion_results.columns:
        for _, row in rxn_deletion_results.iterrows():
            rxn_ids = row["ids"]
            growth = _get_growth(row, rxn_deletion_results.columns)
            if math.isnan(growth):
                growth = 0.0
            for rxn_id in rxn_ids:
                rxn_essential[rxn_id] = growth < threshold
    else:
        for idx, row in rxn_deletion_results.iterrows():
            growth = _get_growth(row, rxn_deletion_results.columns)
            if math.isnan(growth):
                growth = 0.0
            for rxn_id in idx:
                rxn_essential[rxn_id] = growth < threshold

    # For each gene, fraction of its reactions that are essential
    essentiality = {}
    for gene in model.genes:
        rxns = list(gene.reactions)
        if not rxns:
            essentiality[gene.id] = 0.0
        else:
            n_essential = sum(1 for rxn in rxns if rxn_essential.get(rxn.id, False))
            essentiality[gene.id] = n_essential / len(rxns)
    return essentiality


def compute_mvi(
    knockouts_path: str,
    model_path: str,
    baseline_path: str,
    weights: tuple,
    output_path: str,
) -> None:
    if not os.path.exists(knockouts_path):
        raise FileNotFoundError(f"Knockouts file not found: {knockouts_path}")
    if not os.path.exists(model_path):
        raise FileNotFoundError(f"Model file not found: {model_path}")
    if not os.path.exists(baseline_path):
        raise FileNotFoundError(f"Baseline file not found: {baseline_path}")

    w1, w2, w3 = weights

    with open(knockouts_path) as f:
        knockouts = json.load(f)
    with open(baseline_path) as f:
        baseline = json.load(f)

    wt_growth = baseline["growth_rate"]
    total_absolute_flux = baseline["total_absolute_flux"]
    flux_distribution = baseline["flux_distribution"]
    delta_growth_map = {g["gene_id"]: g["delta_growth"] for g in knockouts["genes"]}


    model = cobra.io.load_json_model(model_path)
    model_id = os.path.splitext(os.path.basename(model_path))[0]

    flux_centrality = compute_flux_centrality(model, flux_distribution, total_absolute_flux)
    pathway_essentiality = compute_pathway_essentiality(model, wt_growth)

    genes = []
    for gene in model.genes:
        gene_id = gene.id
        if gene_id not in delta_growth_map:
            print(f"WARNING: gene {gene_id} not in knockouts, using delta_growth=0.0")
        dg = delta_growth_map.get(gene_id, 0.0)
        fc = flux_centrality.get(gene_id, 0.0)
        pe = pathway_essentiality.get(gene_id, 0.0)
        mvi = w1 * dg + w2 * fc + w3 * pe
        genes.append({
            "gene_id": gene_id,
            "gene_name": gene.name if gene.name else gene_id,
            "mvi": round(mvi, 6),
            "delta_growth": round(dg, 6),
            "flux_centrality": round(fc, 6),
            "pathway_essentiality": round(pe, 6),
        })

    genes.sort(key=lambda x: x["mvi"], reverse=True)
    for i, g in enumerate(genes):
        g["rank"] = i + 1

    result = {
        "model_id": model_id,
        "weights": {"w1_delta_growth": w1, "w2_flux_centrality": w2, "w3_pathway_essentiality": w3},
        "n_genes": len(genes),
        "top_20": [g["gene_id"] for g in genes[:20]],
        "genes": genes,
    }

    out_dir = os.path.dirname(output_path)
    if out_dir:
        os.makedirs(out_dir, exist_ok=True)
    with open(output_path, "w") as f:
        json.dump(result, f, indent=2)

    print(f"MVI computed for {len(genes)} genes in {model_id}")
    if genes:
        print(f"Top gene: {genes[0]['gene_id']} (MVI={genes[0]['mvi']:.4f})")
    print(f"Saved to: {output_path}")


def main():
    parser = argparse.ArgumentParser(description="Compute Metabolic Vulnerability Index")
    parser.add_argument("--knockouts", required=True, help="Path to knockouts JSON")
    parser.add_argument("--model", required=True, help="Path to BiGG model JSON")
    parser.add_argument("--baseline", required=True, help="Path to baseline JSON")
    parser.add_argument("--weights", nargs=3, type=float, default=[0.5, 0.3, 0.2],
                        metavar=("W1", "W2", "W3"), help="MVI weights (default: 0.5 0.3 0.2)")
    parser.add_argument("--output", required=True, help="Output JSON path")
    args = parser.parse_args()
    compute_mvi(args.knockouts, args.model, args.baseline, tuple(args.weights), args.output)


if __name__ == "__main__":
    main()
```

Write to `scripts/validate_mvi.py`:

```python
#!/usr/bin/env python3
"""Validate MVI rankings against gold-standard essential gene lists."""
import argparse
import csv
import json
import os
import math

import numpy as np
from sklearn.metrics import roc_auc_score


def load_gold_standard(path: str) -> set:
    """Load gene IDs from CSV gold standard. Reads 'gene_id' column."""
    gene_ids = set()
    with open(path, newline="") as f:
        reader = csv.DictReader(f)
        for row in reader:
            gene_ids.add(row["gene_id"].strip())
    return gene_ids


def precision_at_n(ranked_genes: list, gold_set: set, n: int) -> float:
    top_n = [g["gene_id"] for g in ranked_genes[:n]]
    hits = sum(1 for gid in top_n if gid in gold_set)
    return hits / n if n > 0 else 0.0


def recall_at_n(ranked_genes: list, gold_set: set, n: int) -> float:
    top_n = [g["gene_id"] for g in ranked_genes[:n]]
    hits = sum(1 for gid in top_n if gid in gold_set)
    return hits / len(gold_set) if gold_set else 0.0


def validate(mvi_path: str, gold_standard_path: str, output_path: str) -> None:
    if not os.path.exists(mvi_path):
        raise FileNotFoundError(f"MVI file not found: {mvi_path}")
    if not os.path.exists(gold_standard_path):
        raise FileNotFoundError(f"Gold standard file not found: {gold_standard_path}")

    with open(mvi_path) as f:
        mvi_data = json.load(f)

    gold_set = load_gold_standard(gold_standard_path)
    ranked_genes = mvi_data["genes"]  # already sorted descending by MVI

    # AUC-ROC: MVI score as predictor, gold standard as binary label
    y_score = np.array([g["mvi"] for g in ranked_genes])
    y_true = np.array([1 if g["gene_id"] in gold_set else 0 for g in ranked_genes])

    if y_true.sum() == 0 or y_true.sum() == len(y_true):
        auc_roc = float("nan")
        print("WARNING: Gold standard has no positive/negative examples — AUC undefined")
    else:
        auc_roc = float(roc_auc_score(y_true, y_score))

    result = {
        "model_id": mvi_data["model_id"],
        "gold_standard": os.path.basename(gold_standard_path),
        "n_genes_ranked": len(ranked_genes),
        "n_gold_standard": len(gold_set),
        "auc_roc": round(auc_roc, 4) if not math.isnan(auc_roc) else None,
        "precision_at_10": round(precision_at_n(ranked_genes, gold_set, 10), 4),
        "precision_at_20": round(precision_at_n(ranked_genes, gold_set, 20), 4),
        "precision_at_50": round(precision_at_n(ranked_genes, gold_set, 50), 4),
        "recall_at_10": round(recall_at_n(ranked_genes, gold_set, 10), 4),
        "recall_at_20": round(recall_at_n(ranked_genes, gold_set, 20), 4),
        "recall_at_50": round(recall_at_n(ranked_genes, gold_set, 50), 4),
        "novel_top20": [
            g["gene_id"] for g in ranked_genes[:20]
            if g["gene_id"] not in gold_set
        ],
    }

    out_dir = os.path.dirname(output_path)
    if out_dir:
        os.makedirs(out_dir, exist_ok=True)
    with open(output_path, "w") as f:
        json.dump(result, f, indent=2)

    print(f"Model: {result['model_id']}")
    print(f"AUC-ROC: {result['auc_roc']}")
    print(f"Precision@10/20/50: {result['precision_at_10']} / {result['precision_at_20']} / {result['precision_at_50']}")
    print(f"Novel top-20 candidates: {len(result['novel_top20'])}")
    print(f"Saved to: {output_path}")


def main():
    parser = argparse.ArgumentParser(description="Validate MVI rankings against gold standard")
    parser.add_argument("--mvi", required=True, help="Path to MVI ranked JSON")
    parser.add_argument("--gold-standard", required=True, help="Path to gold standard CSV")
    parser.add_argument("--output", required=True, help="Output JSON path")
    args = parser.parse_args()
    validate(args.mvi, args.gold_standard, args.output)


if __name__ == "__main__":
    main()
```

Write to `scripts/sensitivity_analysis.py`:

```python
#!/usr/bin/env python3
"""Sensitivity analysis: MVI rank stability under ±20% weight perturbations."""
import argparse
import itertools
import json
import os

import numpy as np
from scipy.stats import spearmanr


BASE_WEIGHTS = (0.5, 0.3, 0.2)
MULTIPLIERS = [0.8, 1.0, 1.2]
TOP_N = 100


def rerank_with_weights(genes: list, w1: float, w2: float, w3: float) -> list:
    """Recompute MVI with new weights and return gene IDs sorted by new MVI."""
    scored = []
    for g in genes:
        new_mvi = w1 * g["delta_growth"] + w2 * g["flux_centrality"] + w3 * g["pathway_essentiality"]
        scored.append((g["gene_id"], new_mvi))
    scored.sort(key=lambda x: x[1], reverse=True)
    return [gid for gid, _ in scored]


def run_sensitivity(ecoli_mvi_path: str, output_path: str) -> None:
    if not os.path.exists(ecoli_mvi_path):
        raise FileNotFoundError(f"MVI file not found: {ecoli_mvi_path}")

    with open(ecoli_mvi_path) as f:
        mvi_data = json.load(f)

    genes = sorted(mvi_data["genes"], key=lambda g: g["mvi"], reverse=True)
    baseline_ranking = [g["gene_id"] for g in genes]
    baseline_top100 = baseline_ranking[:TOP_N]
    # Rank positions for baseline top-100 (1-indexed)
    baseline_ranks = {gid: i + 1 for i, gid in enumerate(baseline_top100)}

    perturbations = []
    rho_values = []

    for mult in itertools.product(MULTIPLIERS, repeat=3):
        # Skip pure baseline (1.0, 1.0, 1.0)
        if mult == (1.0, 1.0, 1.0):
            continue

        # Scale and renormalize weights
        raw = [BASE_WEIGHTS[i] * mult[i] for i in range(3)]
        total = sum(raw)
        w1, w2, w3 = [r / total for r in raw]

        perturbed_ranking = rerank_with_weights(genes, w1, w2, w3)
        perturbed_top100 = perturbed_ranking[:TOP_N]

        # Spearman rho on rank positions of baseline top-100 genes
        base_pos = list(range(1, TOP_N + 1))
        perturbed_pos = []
        for gid in baseline_top100:
            try:
                pos = perturbed_top100.index(gid) + 1
            except ValueError:
                pos = TOP_N + 1  # gene fell out of top-100
            perturbed_pos.append(pos)

        rho, _ = spearmanr(base_pos, perturbed_pos)
        rho_values.append(float(rho))
        perturbations.append({
            "multipliers": list(mult),
            "weights": [round(w1, 4), round(w2, 4), round(w3, 4)],
            "spearman_rho": round(float(rho), 4),
        })

    result = {
        "base_weights": list(BASE_WEIGHTS),
        "multiplier_grid": MULTIPLIERS,
        "top_n": TOP_N,
        "n_perturbations": len(perturbations),
        "mean_spearman_rho": round(float(np.mean(rho_values)), 4),
        "min_spearman_rho": round(float(np.min(rho_values)), 4),
        "max_spearman_rho": round(float(np.max(rho_values)), 4),
        "std_spearman_rho": round(float(np.std(rho_values)), 4),
        "robust": float(np.min(rho_values)) > 0.85,
        "perturbations": perturbations,
    }

    out_dir = os.path.dirname(output_path)
    if out_dir:
        os.makedirs(out_dir, exist_ok=True)
    with open(output_path, "w") as f:
        json.dump(result, f, indent=2)

    print(f"Sensitivity analysis: {len(perturbations)} perturbations")
    print(f"Spearman rho -- mean: {result['mean_spearman_rho']:.4f}, min: {result['min_spearman_rho']:.4f}")
    print(f"Robust (min rho > 0.85): {result['robust']}")
    print(f"Saved to: {output_path}")


def main():
    parser = argparse.ArgumentParser(description="MVI sensitivity analysis")
    parser.add_argument("--ecoli-mvi", required=True, help="Path to E. coli MVI ranked JSON")
    parser.add_argument("--ecoli-knockouts", help="Unused (reserved for future)")
    parser.add_argument("--ecoli-baseline", help="Unused (reserved for future)")
    parser.add_argument("--ecoli-model", help="Unused (reserved for future)")
    parser.add_argument("--output", required=True, help="Output JSON path")
    args = parser.parse_args()
    run_sensitivity(args.ecoli_mvi, args.output)


if __name__ == "__main__":
    main()
```

Write to `scripts/verify_hashes.py`:

```python
#!/usr/bin/env python3
"""Verify SHA-256 hashes of all pipeline artifacts. PASS/FAIL per artifact."""
import argparse
import hashlib
import json
import os
import sys

ARTIFACTS = [
    "artifacts/ecoli_baseline.json",
    "artifacts/mtb_baseline.json",
    "artifacts/ecoli_knockouts.json",
    "artifacts/mtb_knockouts.json",
    "artifacts/ecoli_mvi_ranked.json",
    "artifacts/mtb_mvi_ranked.json",
    "artifacts/ecoli_validation.json",
    "artifacts/mtb_validation.json",
    "artifacts/sensitivity_analysis.json",
]

REFERENCE_PATH = "artifacts/reference_hashes.json"


def sha256_file(path: str) -> str:
    h = hashlib.sha256()
    with open(path, "rb") as f:
        for chunk in iter(lambda: f.read(65536), b""):
            h.update(chunk)
    return h.hexdigest()


def generate_hashes(base_dir: str) -> dict:
    hashes = {}
    for artifact in ARTIFACTS:
        full_path = os.path.join(base_dir, artifact)
        if not os.path.exists(full_path):
            print(f"MISSING: {artifact}")
            continue
        hashes[artifact] = sha256_file(full_path)
        print(f"Hashed: {artifact}")
    return hashes


def verify_hashes(base_dir: str) -> bool:
    ref_path = os.path.join(base_dir, REFERENCE_PATH)
    if not os.path.exists(ref_path):
        print(f"ERROR: Reference hashes not found at {ref_path}", file=sys.stderr)
        return False

    with open(ref_path) as f:
        reference = json.load(f)

    all_pass = True
    for artifact in ARTIFACTS:
        full_path = os.path.join(base_dir, artifact)
        if artifact not in reference:
            print(f"SKIP (not in reference): {artifact}")
            continue
        if not os.path.exists(full_path):
            print(f"FAIL (missing): {artifact}")
            all_pass = False
            continue
        actual = sha256_file(full_path)
        expected = reference[artifact]
        if actual == expected:
            print(f"PASS: {artifact}")
        else:
            print(f"FAIL: {artifact}")
            print(f"  expected: {expected}")
            print(f"  actual:   {actual}")
            all_pass = False

    return all_pass


def main():
    parser = argparse.ArgumentParser(description="Verify artifact SHA-256 hashes")
    parser.add_argument("--generate", action="store_true",
                        help="Generate reference hashes from current artifacts")
    parser.add_argument("--base-dir", default=".",
                        help="Base directory (default: current directory)")
    args = parser.parse_args()

    if args.generate:
        hashes = generate_hashes(args.base_dir)
        ref_path = os.path.join(args.base_dir, REFERENCE_PATH)
        with open(ref_path, "w") as f:
            json.dump(hashes, f, indent=2)
        print(f"\nReference hashes saved to {ref_path}")
        print(f"Hashed {len(hashes)} artifacts")
    else:
        passed = verify_hashes(args.base_dir)
        if passed:
            print("\nAll artifacts PASS hash verification.")
            sys.exit(0)
        else:
            print("\nHash verification FAILED.", file=sys.stderr)
            sys.exit(1)


if __name__ == "__main__":
    main()
```

---

## Step 3: Download metabolic models

```bash
python scripts/download_models.py --organisms ecoli mtb --outdir data/models/
```

**Expected output:**
```
[download] http://bigg.ucsd.edu/static/models/iJO1366.json -> data/models/iJO1366.json
[ok] iJO1366: 2583 reactions, 1367 genes
[download] http://bigg.ucsd.edu/static/models/iEK1008.json -> data/models/iEK1008.json
[ok] iEK1008: 1226 reactions, 1008 genes
[done] All models ready.
```

Files produced: `data/models/iJO1366.json` (2,583 reactions, 1,367 genes), `data/models/iEK1008.json` (1,226 reactions, 1,008 genes)

---

## Step 4: Run FBA baseline

```bash
python scripts/run_fba_baseline.py --model data/models/iJO1366.json --output artifacts/ecoli_baseline.json
python scripts/run_fba_baseline.py --model data/models/iEK1008.json --output artifacts/mtb_baseline.json
```

**Expected output:**
```
Model: iJO1366
Growth rate: 0.9824
Genes: 1367, Reactions: 2583

Model: iEK1008
Growth rate: 0.0582
Genes: 1008, Reactions: 1226
```

Expected: *E. coli* growth = 0.9824 h⁻¹, *M. tuberculosis* growth = 0.0582 h⁻¹

---

## Step 5: Gene knockout screen

```bash
python scripts/run_knockout_screen.py \
    --model data/models/iJO1366.json \
    --baseline artifacts/ecoli_baseline.json \
    --output artifacts/ecoli_knockouts.json

python scripts/run_knockout_screen.py \
    --model data/models/iEK1008.json \
    --baseline artifacts/mtb_baseline.json \
    --output artifacts/mtb_knockouts.json
```

**Expected output:**
```
Running single-gene knockout screen for iJO1366 (1367 genes)...
Screened: 1367 genes, Essential: 208

Running single-gene knockout screen for iEK1008 (1008 genes)...
Screened: 1008 genes, Essential: 270
```

Expected: *E. coli* 208 essential / 1,367 screened; MTB 270 / 1,008 screened

---

## Step 6: MVI scoring

```bash
python scripts/compute_mvi.py \
    --knockouts artifacts/ecoli_knockouts.json \
    --model data/models/iJO1366.json \
    --baseline artifacts/ecoli_baseline.json \
    --weights 0.5 0.3 0.2 \
    --output artifacts/ecoli_mvi_ranked.json

python scripts/compute_mvi.py \
    --knockouts artifacts/mtb_knockouts.json \
    --model data/models/iEK1008.json \
    --baseline artifacts/mtb_baseline.json \
    --weights 0.5 0.3 0.2 \
    --output artifacts/mtb_mvi_ranked.json
```

**Expected output:**
```
MVI computed for 1367 genes in iJO1366
Top gene: b0720 (MVI=0.7020)

MVI computed for 1008 genes in iEK1008
Top gene: Rv1310 (MVI=0.7174)
```

Expected: *E. coli* top gene b0720 (gltA), MVI = 0.7020; MTB top gene Rv1310 (atpD), MVI = 0.7174

---

## Step 7: Validation against gold-standard datasets

```bash
python scripts/validate_mvi.py \
    --mvi artifacts/ecoli_mvi_ranked.json \
    --gold-standard data/keio_essential_genes.csv \
    --output artifacts/ecoli_validation.json

python scripts/validate_mvi.py \
    --mvi artifacts/mtb_mvi_ranked.json \
    --gold-standard data/tb_drug_targets.csv \
    --output artifacts/mtb_validation.json
```

**Expected output:**
```
Model: iJO1366
AUC-ROC: 0.5569
Precision@10/20/50: 0.30 / ...
Novel top-20 candidates: 17

Model: iEK1008
AUC-ROC: 0.7806
Precision@10/20/50: 0.10 / ...
Novel top-20 candidates: 19
```

Expected: *E. coli* AUC-ROC = 0.5569, Precision@10 = 0.30, novel top-20 = 17 genes; MTB AUC-ROC = 0.7806, Precision@10 = 0.10, novel top-20 = 19 genes

---

## Step 8: Sensitivity analysis

```bash
python scripts/sensitivity_analysis.py \
    --ecoli-mvi artifacts/ecoli_mvi_ranked.json \
    --output artifacts/sensitivity_analysis.json
```

**Expected output:**
```
Sensitivity analysis: 26 perturbations
Spearman rho -- mean: 0.9997, min: 0.9997
Robust (min rho > 0.85): True
```

Expected: mean Spearman ρ = 0.9997, min = 0.9997, robust = True

---

## Step 9: Verify reproducibility

On first clean run, generate reference hashes:

```bash
python scripts/verify_hashes.py --generate --base-dir .
```

This saves `artifacts/reference_hashes.json`. On subsequent runs, verify against them:

```bash
python scripts/verify_hashes.py --base-dir .
```

**Expected output:**
```
PASS: artifacts/ecoli_baseline.json
PASS: artifacts/mtb_baseline.json
PASS: artifacts/ecoli_knockouts.json
PASS: artifacts/mtb_knockouts.json
PASS: artifacts/ecoli_mvi_ranked.json
PASS: artifacts/mtb_mvi_ranked.json
PASS: artifacts/ecoli_validation.json
PASS: artifacts/mtb_validation.json
PASS: artifacts/sensitivity_analysis.json

All artifacts PASS hash verification.
```

Expected: All 9 artifacts PASS hash verification

---

## Summary of expected outputs

| Artifact | Key metric |
|---|---|
| `artifacts/ecoli_baseline.json` | growth_rate = 0.9824 h⁻¹, 1,367 genes, 2,583 reactions |
| `artifacts/mtb_baseline.json` | growth_rate = 0.0582 h⁻¹, 1,008 genes, 1,226 reactions |
| `artifacts/ecoli_knockouts.json` | 208 essential / 1,367 screened |
| `artifacts/mtb_knockouts.json` | 270 essential / 1,008 screened |
| `artifacts/ecoli_mvi_ranked.json` | top gene b0720 (gltA), MVI = 0.7020 |
| `artifacts/mtb_mvi_ranked.json` | top gene Rv1310 (atpD), MVI = 0.7174 |
| `artifacts/ecoli_validation.json` | AUC-ROC = 0.5569, Precision@10 = 0.30, novel top-20 = 17 |
| `artifacts/mtb_validation.json` | AUC-ROC = 0.7806, Precision@10 = 0.10, novel top-20 = 19 |
| `artifacts/sensitivity_analysis.json` | mean Spearman ρ = 0.9997, min = 0.9997, robust = True |