Utilizing subgraph isomorphism algorithms on protein-protein interaction networks derived from antibiotic-resistant bacterial strains will reveal conserved structural motifs that correlate with specific evolutionary fitness trade-offs.

Subgraph isomorphism algorithms applied to protein-protein interaction (PPI) networks constructed from antibiotic-resistant bacterial strains will identify statistically enriched conserved structural motifs (subgraph patterns recurring at frequency ≥2× background in ≥3 independent resistant strains) that show significant correlation (Spearman |ρ| ≥ 0.40, p < 0.05) with quantifiable evolutionary fitness trade-offs — specifically: (1) growth rate penalty (doubling time increase ≥10% vs. susceptible isogenic controls), (2) competitive fitness coefficient in mixed-culture assays, and (3) cross-resistance or collateral sensitivity profiles across ≥2 antibiotic classes. The hypothesis is falsifiable: if no motif class shows frequency enrichment beyond random expectation (permutation-corrected p > 0.05) or if enriched motifs fail to correlate with any measured fitness phenotype, the hypothesis is rejected.

MINIMUM VIABLE TEST (MVT) — 3-phase design targeting 90-day completion:

PHASE A — Network Construction (Days 1–30): Collect whole-genome sequences + experimentally validated PPI data for 20 strains: 10 antibiotic-resistant (≥2 resistance classes each), 10 isogenic or near-isogenic susceptible controls. Species: E. coli (n=8), K. pneumoniae (n=6), P. aeruginosa (n=6). Source networks from STRING v12 (confidence ≥700) filtered to organism-specific experimental evidence. Augment with co-immunoprecipitation data from BioGRID (release 4.4+). Construct strain-specific PPI networks as undirected weighted graphs. Minimum network size threshold: 500 nodes, 1,500 edges per strain.

PHASE B — Subgraph Isomorphism & Motif Enumeration (Days 15–60): Apply VF2++ algorithm (igraph 0.10+ implementation) for exact subgraph isomorphism at k=3,4,5 nodes. For k=6,7, apply FANMOD approximate enumeration (10^6 random subgraph samples). Enumerate all non-isomorphic connected subgraph patterns. Compute motif significance profile (MSP) for each strain: Z-score of observed vs. 1,000 degree-sequence-preserving random networks (Erdős–Rényi null with matched degree distribution). Identify motifs with Z > 2.0 in ≥60% of resistant strains and Z < 1.0 in ≥60% of susceptible controls.

PHASE C — Fitness Phenotyping & Correlation (Days 1–90, parallel): Measure three fitness proxies for all 20 strains: (1) growth rate in LB at 37°C (OD600 kinetics, 24h, n=3 biological replicates); (2) competitive fitness coefficient vs. fluorescently labeled susceptible reference strain (1:1 co-culture, 24h, flow cytometry ratio); (3) cross-resistance/collateral sensitivity profile (MIC panel: 8 antibiotics across 4 classes, EUCAST breakpoints). Compute fitness trade-off vector per strain. Correlate motif frequency matrix (strains × motif classes) with fitness vector using Spearman ρ, FDR correction (Benjamini-Hochberg).

# ============================================================
# BACTERIAL PPI MOTIF-FITNESS PIPELINE
# Subgraph Isomorphism → Motif Enrichment → Fitness Correlation
# ============================================================

# --- CONFIGURATION ---
CONFIG = {
    "species": ["Escherichia_coli", "Klebsiella_pneumoniae", "Pseudomonas_aeruginosa"],
    "n_strains_resistant": 10,
    "n_strains_susceptible": 10,
    "string_confidence_threshold": 700,
    "motif_sizes": [3, 4, 5, 6, 7],  # k values
    "n_null_permutations": 1000,
    "fanmod_samples": 1_000_000,  # for k >= 6
    "enrichment_fdr_threshold": 0.05,
    "correlation_threshold": 0.40,
    "random_seed": 42
}

# --- PHASE 1: NETWORK CONSTRUCTION ---
def build_strain_ppi_networks(strain_list, config):
    """
    For each strain: load STRING + BioGRID interactions,
    filter by confidence, map to strain-specific proteome,
    return dict of igraph Graph objects.
    """
    networks = {}
    for strain in strain_list:
        # Load STRING experimental interactions
        string_edges = load_string_interactions(
            organism=strain.species,
            confidence_min=config["string_confidence_threshold"],
            evidence_channels=["experimental", "database"]
        )
        # Map to strain-specific proteins via DIAMOND BLASTp
        strain_proteins = map_proteome_to_string(
            genome_fasta=strain.genome_path,
            reference_proteome=strain.species,
            identity_threshold=0.40,
            coverage_threshold=0.80
        )
        # Filter edges to strain-specific proteins
        strain_edges = [(u, v, w) for u, v, w in string_edges
                        if u in strain_proteins and v in strain_proteins]
        # Augment with BioGRID
        biogrid_edges = load_biogrid_interactions(organism=strain.species)
        all_edges = merge_deduplicate(strain_edges, biogrid_edges)
        # Build igraph Graph
        G = igraph.Graph.TupleList(all_edges, weights=True, directed=False)
        # Quality check
        assert G.vcount() >= 500, f"Network too small: {G.vcount()} nodes"
        assert G.ecount() >= 1500, f"Network too sparse: {G.ecount()} edges"
        networks[strain.id] = G
    return networks

# --- PHASE 2: NULL MODEL GENERATION ---
def generate_null_ensemble(G, n_permutations=1000, n_swaps_multiplier=100):
    """
    Maslov-Sneppen edge swap preserving degree sequence.
    Returns list of n_permutations random graphs.
    """
    null_graphs = []
    for i in range(n_permutations):
        G_null = G.copy()
        n_swaps = G.ecount() * n_swaps_multiplier
        G_null = edge_swap_rewire(G_null, n_swaps=n_swaps)
        # Verify degree sequence preserved
        assert sorted(G_null.degree()) == sorted(G.degree())
        null_graphs.append(G_null)
    return null_graphs

# --- PHASE 3: SUBGRAPH ENUMERATION ---
def enumerate_motifs_exact(G, k_values=[3, 4, 5]):
    """
    VF2++ exact subgraph isomorphism for k=3,4,5.
    Returns dict: {canonical_label: count}
    """
    motif_counts = {}
    for k in k_values:
        # Generate all non-isomorphic connected graphs of size k
        reference_motifs = get_canonical_motifs(k)  # nauty-generated
        for motif_template in reference_motifs:
            canonical_label = f"k{k}_{motif_template.canonical_hash}"
            count = count_subgraph_isomorphisms_vf2pp(G, motif_template)
            motif_counts[canonical_label] = count
    return motif_counts

def enumerate_motifs_approximate(G, k_values=[6, 7], n_samples=1_000_000):
    """
    FANMOD-style random subgraph sampling for k=6,7.
    Returns dict: {canonical_label: frequency}
    """
    motif_frequencies = {}
    for k in k_values:
        # Random subgraph sampling
        sampled_subgraphs = random_subgraph_sample(G, k=k, n_samples=n_samples)
        # Canonicalize each sample
        for sg in sampled_subgraphs:
            label = nauty_canonical_form(sg)
            motif_frequencies[label] = motif_frequencies.get(label, 0) + 1
        # Normalize to frequency
        total = sum(motif_frequencies.values())
        motif_frequencies = {k: v/total for k, v in motif_frequencies.items()}
    return motif_frequencies

# --- PHASE 4: Z-SCORE COMPUTATION ---
def compute_motif_zscore_profile(G, null_ensemble, k_values):
    """
    For each motif class: Z = (observed - mean_null) / std_null
    """
    observed = enumerate_motifs_exact(G, k_values=[k for k in k_values if k <= 5])
    observed.update(enumerate_motifs_approximate(G, k_values=[k for k in k_values if k > 5]))
    
    null_distributions = {}
    for null_G in null_ensemble:
        null_counts = enumerate_motifs_exact(null_G, k_values=[k for k in k_values if k <= 5])
        null_counts.update(enumerate_motifs_approximate(null_G, k_values=[k for k in k_values if k > 5]))
        for motif, count in null_counts.items():
            null_distributions.setdefault(motif, []).append(count)
    
    z_scores = {}
    for motif, obs_count in observed.items():
        null_vals = null_distributions.get(motif, [0] * len(null_ensemble))
        mu, sigma = np.mean(null_vals), np.std(null_vals)
        z_scores[motif] = (obs_count - mu) / sigma if sigma > 0 else 0.0
    return z_scores

# --- PHASE 5: ENRICHMENT ANALYSIS ---
def identify_resistance_enriched_motifs(resistant_zscores, susceptible_zscores, fdr_threshold=0.05):
    """
    Mann-Whitney U test + BH FDR correction.
    Returns list of resistance-enriched motif (REM) labels.
    """
    all_motifs = set(list(resistant_zscores[0].keys()))
    p_values = {}
    for motif in all_motifs:
        r_vals = [zs.get(motif, 0) for zs in resistant_zscores]
        s_vals = [zs.get(motif, 0) for zs in susceptible_zscores]
        _, p = scipy.stats.mannwhitneyu(r_vals, s_vals, alternative='greater')
        p_values[motif] = p
    
    # BH FDR correction
    motifs, pvals = zip(*p_values.items())
    _, fdr_corrected, _, _ = statsmodels.stats.multitest.multipletests(pvals, method='fdr_bh')
    
    # Filter: FDR < threshold AND Z_resistant > 2.0 AND Z_susceptible < 1.0
    REMs = [m for m, fdr, z_r, z_s in zip(
        motifs, fdr_corrected,
        [np.mean([zs.get(m, 0) for zs in resistant_zscores]) for m in motifs],
        [np.mean([zs.get(m, 0) for zs in susceptible_zscores]) for m in motifs]
    ) if fdr < fdr_threshold and z_r > 2.0 and z_s < 1.0]
    return REMs

# --- PHASE 6: FITNESS CORRELATION ---
def correlate_motifs_with_fitness(motif_matrix, fitness_matrix, REMs, alpha=0.05):
    """
    Spearman correlation with Bonferroni correction.
    motif_matrix: (n_strains x n_motifs) DataFrame
    fitness_matrix: (n_strains x 3) DataFrame [growth_rate, competitive_fitness, resistance_breadth]
    """
    results = []
    n_tests = len(REMs) * fitness_matrix.shape[1]
    for motif in REMs:
        for fitness_metric in fitness_matrix.columns:
            rho, p = scipy.stats.spearmanr(
                motif_matrix[motif], fitness_matrix[fitness_metric]
            )
            p_corrected = min(p * n_tests, 1.0)  # Bonferroni
            results.append({
                "motif": motif,
                "fitness_metric": fitness_metric,
                "spearman_rho": rho,
                "p_raw": p,
                "p_bonferroni": p_corrected,
                "significant": p_corrected < alpha and abs(rho) >= 0.40
            })
    return pd.DataFrame(results)

# --- PHASE 7: PHYLOGENETIC CORRECTION (R subprocess) ---
def run_pgls_correction(motif_vector, fitness_vector, phylo_tree_newick):
    """
    Calls R script for PGLS via subprocess.
    Returns p-value and lambda estimate.
    """
    r_script = f"""
    library(ape); library(nlme)
    tree <- read.tree(text="{phylo_tree_newick}")
    data <- data.frame(motif={list(motif_vector)}, fitness={list(fitness_vector)})
    rownames(data) <- tree$tip.label
    pgls_model <- gls(fitness ~ motif, data=data,
                      correlation=corBrownian(phy=tree),
                      method="ML")
    summary(pgls_model)$tTable["motif", "p-value"]
    """
    result = subprocess.run(["Rscript", "-e", r_script], capture_output=True, text=True)
    p_pgls = float(result.stdout.strip())
    return p_pgls

# --- MAIN PIPELINE ---
def main():
    np.random.seed(CONFIG["random_seed"])
    
    # Load strain metadata
    resistant_strains = load_strain_metadata("resistant", n=10)
    susceptible_strains = load_strain_metadata("susceptible", n=10)
    all_strains = resistant_strains + susceptible_strains
    
    # Build networks
    print("Building PPI networks...")
    networks = build_strain_ppi_networks(all_strains, CONFIG)
    
    # Generate null ensembles (parallelized)
    print("Generating null models...")
    null_ensembles = {sid: generate_null_ensemble(G, CONFIG["n_null_permutations"])
                     for sid, G in networks.items()}
    
    # Compute Z-score profiles
    print("Enumerating motifs and computing Z-scores...")
    zscores = {sid: compute_motif_zscore_profile(G, null_ensembles[sid], CONFIG["motif_sizes"])
               for sid, G in networks.items()}
    
    # Identify REMs
    resistant_zscores = [zscores[s.id] for s in resistant_strains]
    susceptible_zscores = [zscores[s.id] for s in susceptible_strains]
    REMs = identify_resistance_enriched_motifs(resistant_zscores, susceptible_zscores)
    print(f"Identified {len(REMs)} resistance-enriched motifs")
    
    # Fitness phenotyping (external data loaded from lab measurements)
    fitness_df = load_fitness_measurements(all_strains)
    
    # Build motif frequency matrix
    motif_df = pd.DataFrame(zscores).T  # strains x motifs
    
    # Correlate
    correlation_results = correlate_motifs_with_fitness(motif_df, fitness_df, REMs)
    significant_pairs = correlation_results[correlation_results["significant"]]
    print(f"Significant motif-fitness correlations: {len(significant_pairs)}")
    
    # Phylogenetic correction for significant pairs
    phylo_tree = build_core_genome_phylogeny(all_strains)
    for _, row in significant_pairs.iterrows():
        p_pgls = run_pgls_correction(
            motif_df[row["motif"]], fitness_df[row["fitness_metric"]], phylo_tree
        )
        print(f"PGLS p-value for {row['motif']} ~ {row['fitness_metric']}: {p_pgls:.4f}")
    
    # Export results
    correlation_results.to_csv("results/motif_fitness_correlations.csv", index=False)
    export_network_visualizations(networks, REMs, "results/cytoscape/")
    print("Pipeline complete. Results in results/")

if __name__ == "__main__":
    main()

# ============================================================
# COMPUTE RESOURCE ESTIMATES:
# Network construction: ~2h CPU per strain × 20 = 40h CPU
# Null model generation: ~8h CPU per strain × 20 = 160h CPU
# Motif enumeration k=3-5: ~4h CPU per strain × 20 = 80h CPU
# Motif enumeration k=6-7: ~8h CPU per strain × 20 = 160h CPU
# Correlation + PGLS: ~2h CPU total = 2h CPU
# GPU: optional GNN extension only = 12h GPU
# Peak RAM: 128 GB (null ensemble storage for largest networks)
# ============================================================