The co-occurrence network of resistance genes across WHO GLASS genomic isolates exhibits scale-free topology with identifiable hub genes — those whose presence predicts the broadest range of co-occurring resistance phenotypes. Disrupting the protein products of the top-5 hub genes with small-molecule inhibitors will collapse multi-drug resistance phenotypes across phylogenetically unrelated pathogens simultaneously, because these hubs represent regulatory or assembly factors for resistance gene cassettes rather than individual resistance mechanisms; this network-medicine approach will be validated by showing that hub-gene knockdowns restore susceptibility to ≥3 structurally unrelated antibiotic classes in ESKAPEE organisms.

The co-occurrence network of resistance genes (RGs) constructed from WHO GLASS genomic isolates (≥10,000 isolates, ≥5 ESKAPEE species) exhibits statistically significant scale-free topology (power-law degree distribution, exponent γ between 2 and 3, R² ≥ 0.85 by maximum-likelihood fit). The top-5 hub genes by betweenness centrality and degree are identifiable and encode proteins that function as regulatory or assembly factors for resistance gene cassettes rather than direct resistance effectors. Small-molecule inhibitors targeting these 5 hub gene products will restore susceptibility (MIC reduction ≥4-fold) to ≥3 structurally unrelated antibiotic classes simultaneously in ≥4 of 7 ESKAPEE organisms, as measured by broth microdilution and confirmed by RNA-seq showing downregulation of ≥50% of co-occurring resistance genes within the hub gene's network neighborhood.

PHASE 1 — NETWORK CONSTRUCTION AND TOPOLOGICAL ANALYSIS (Computational, Weeks 1–8): Construct binary presence/absence matrix of resistance genes from WHO GLASS WGS data using AMRFinder+ and CARD. Build co-occurrence network with phi-coefficient edge weights (threshold φ ≥ 0.3, p < 0.001 Bonferroni-corrected). Test scale-free topology by maximum-likelihood power-law fitting (powerlaw Python package). Identify top-5 hub genes by composite score: (betweenness centrality × degree × eigenvector centrality). Validate hub stability across bootstrap resamples (n=1,000) and two GLASS release years.

PHASE 2 — HUB GENE FUNCTIONAL CHARACTERIZATION (Wet lab, Weeks 6–16): Annotate hub gene products using HHpred, InterPro, AlphaFold2 structural prediction, and literature. Clone hub genes into expression vectors. Perform bacterial two-hybrid assays to confirm protein-protein interactions with resistance gene products. Use CRISPRi knockdown in 7 ESKAPEE representative strains to confirm essentiality and resistance gene co-regulation.

PHASE 3 — INHIBITOR IDENTIFICATION AND VALIDATION (Wet lab + Computational, Weeks 12–28): Virtual screen FDA-approved compound library + ChemBridge diversity set (≥500,000 compounds) against AlphaFold2 hub protein structures using AutoDock-GPU. Select top-50 compounds per hub. Validate binding by SPR (KD < 1 μM threshold). Test MIC restoration in ESKAPEE panel (broth microdilution, CLSI guidelines). Confirm mechanism by RNA-seq (n=3 biological replicates per organism per inhibitor).

PHASE 4 — IN VIVO PROOF-OF-CONCEPT (Weeks 24–40): Murine infection models (neutropenic mouse, K. pneumoniae and A. baumannii) with hub inhibitor + antibiotic combination. CFU reduction endpoint. Toxicity panel (ALT, AST, creatinine, CBC).

# PHASE 1: Network Construction Pipeline
# =======================================

import numpy as np
import pandas as pd
import networkx as nx
from scipy.stats import pearsonr
import powerlaw

def build_resistance_network(glass_wgs_dir, amrfinder_db):
    """
    Build co-occurrence network from WHO GLASS WGS data.
    Returns: NetworkX graph G, hub_scores DataFrame
    """
    
    # Step 1: AMR gene detection
    # Run AMRFinder+ on all assemblies (parallelized, 200 CPU cores)
    # amrfinder -n {assembly.fasta} -d {amrfinder_db} --plus -o {output.tsv}
    
    presence_matrix = load_amrfinder_results(glass_wgs_dir)
    # Shape: (n_isolates, n_genes), binary 0/1
    # Expected: ~15,000 isolates × ~2,000 resistance genes
    
    # Step 2: Phylogenetic dereplication
    ani_matrix = compute_fastani_pairwise(glass_wgs_dir)  # GPU-accelerated
    derep_indices = greedy_dereplication(ani_matrix, threshold=0.995)
    presence_matrix = presence_matrix[derep_indices]
    # Expected after dereplication: ~8,000–10,000 isolates
    
    # Step 3: Filter invariant genes
    gene_freq = presence_matrix.mean(axis=0)
    variable_genes = (gene_freq >= 0.01) & (gene_freq <= 0.99)
    presence_matrix = presence_matrix[:, variable_genes]
    
    # Step 4: Compute phi-coefficients (vectorized)
    n = presence_matrix.shape[0]
    phi_matrix = np.zeros((presence_matrix.shape[1], presence_matrix.shape[1]))
    for i in range(presence_matrix.shape[1]):
        for j in range(i+1, presence_matrix.shape[1]):
            phi, p = phi_coefficient(presence_matrix[:, i], presence_matrix[:, j])
            p_bonferroni = p * (presence_matrix.shape[1] ** 2 / 2)
            if phi >= 0.3 and p_bonferroni < 0.05:
                phi_matrix[i, j] = phi_matrix[j, i] = phi
    
    # Step 5: Build NetworkX graph
    G = nx.from_numpy_array(phi_matrix)
    G = nx.relabel_nodes(G, {i: gene_names[i] for i in range(len(gene_names))})
    
    return G

def test_scale_free_topology(G):
    """
    Test power-law degree distribution.
    Returns: fit object, R², gamma, p_value
    """
    degrees = [d for n, d in G.degree()]
    fit = powerlaw.Fit(degrees, discrete=True)
    
    # Likelihood ratio test: power-law vs. exponential
    R, p = fit.distribution_compare('power_law', 'exponential', normalized_ratio=True)
    
    # KS goodness-of-fit
    ks_stat, ks_p = fit.power_law.KS()
    
    gamma = fit.power_law.alpha
    xmin = fit.power_law.xmin
    
    # Bootstrap confidence intervals (n=1000)
    gamma_bootstrap = bootstrap_powerlaw_fit(degrees, n_iter=1000)
    gamma_ci = np.percentile(gamma_bootstrap, [2.5, 97.5])
    
    return {
        'gamma': gamma,
        'gamma_ci': gamma_ci,
        'xmin': xmin,
        'LR_vs_exponential': R,
        'LR_p_value': p,
        'KS_stat': ks_stat,
        'KS_p': ks_p,
        'scale_free_confirmed': (R > 0 and p < 0.05 and 2.0 <= gamma <= 3.0)
    }

def compute_hub_scores(G):
    """
    Compute composite hub scores for all nodes.
    Returns: DataFrame sorted by hub_score descending
    """
    degree_cent = nx.degree_centrality(G)
    between_cent = nx.betweenness_centrality(G, k=500, normalized=True)  # k=500 pivots
    eigen_cent = nx.eigenvector_centrality(G, max_iter=1000)
    
    # Normalize each metric to [0,1]
    def normalize(d):
        vals = np.array(list(d.values()))
        vals_norm = (vals - vals.min()) / (vals.max() - vals.min() + 1e-10)
        return dict(zip(d.keys(), vals_norm))
    
    deg_n = normalize(degree_cent)
    bet_n = normalize(between_cent)
    eig_n = normalize(eigen_cent)
    
    hub_scores = {}
    for gene in G.nodes():
        hub_scores[gene] = (deg_n[gene] * bet_n[gene] * eig_n[gene]) ** (1/3)
    
    df = pd.DataFrame.from_dict(hub_scores, orient='index', columns=['hub_score'])
    df['degree'] = [degree_cent[g] * (G.number_of_nodes()-1) for g in df.index]
    df['betweenness'] = [between_cent[g] for g in df.index]
    df['eigenvector'] = [eigen_cent[g] for g in df.index]
    df = df.sort_values('hub_score', ascending=False)
    
    return df

def bootstrap_hub_stability(presence_matrix, gene_names, n_iter=1000, sample_frac=0.8):
    """
    Bootstrap hub score stability assessment.
    Returns: Spearman rho for top-20 hub rank order
    """
    n_isolates = presence_matrix.shape[0]
    hub_rank_matrix = []
    
    for i in range(n_iter):
        sample_idx = np.random.choice(n_isolates, 
                                       size=int(n_isolates * sample_frac), 
                                       replace=True)
        G_boot = build_resistance_network_from_matrix(
            presence_matrix[sample_idx], gene_names)
        hub_df = compute_hub_scores(G_boot)
        hub_rank_matrix.append(hub_df.index[:20].tolist())
    
    # Compute pairwise Spearman correlations
    rho_values = compute_rank_correlations(hub_rank_matrix)
    return np.mean(rho_values), np.std(rho_values)

# PHASE 2: Virtual Screening Pipeline
# =====================================

def virtual_screen_hub_proteins(hub_gene_list, compound_library_smiles):
    """
    AutoDock-GPU virtual screening against hub protein structures.
    """
    results = {}
    
    for hub_gene in hub_gene_list:
        # Load AlphaFold2 structure (pre-downloaded from EBI)
        structure_pdb = f"alphafold/{hub_gene}_AF2.pdb"
        
        # Prepare receptor (AutoDock Tools)
        # autodock_tools_prepare_receptor -r {structure_pdb} -o {hub_gene}.pdbqt
        
        # Identify binding pocket (fpocket)
        pockets = run_fpocket(structure_pdb)
        best_pocket = pockets[0]  # Highest druggability score
        
        # Run AutoDock-GPU (GPU-accelerated, ~500k compounds)
        # autodock_gpu --ffile {hub_gene}.maps.fld 
        #              --lfile {compound_library}.pdbqt
        #              --nrun 100 --heurmax 12000000
        
        docking_results = parse_autodock_results(f"{hub_gene}_docking_output/")
        
        # ADMET filtering (SwissADME API)
        filtered = apply_admet_filters(docking_results, 
                                        logP_range=(1, 5),
                                        MW_max=500,
                                        remove_PAINS=True)
        
        # Select top-50 by docking score
        top50 = filtered.nsmallest(50, 'docking_score')
        results[hub_gene] = top50
    
    return results

# PHASE 3: MIC Restoration Analysis
# ===================================

def analyze_mic_restoration(mic_data):
    """
    Analyze MIC fold-changes after hub gene knockdown or inhibitor treatment.
    mic_data: DataFrame with columns [organism, hub_gene, antibiotic, 
                                       MIC_control, MIC_treatment]
    Returns: Summary statistics and success/failure determination
    """
    mic_data['fold_change'] = mic_data['MIC_control'] / mic_data['MIC_treatment']
    mic_data['restored'] = mic_data['fold_change'] >= 4.0
    
    # Count antibiotic classes restored per organism per hub
    class_map = {
        'meropenem': 'beta_lactam', 'ampicillin': 'beta_lactam',
        'gentamicin': 'aminoglycoside', 'amikacin': 'aminoglycoside',
        'ciprofloxacin': 'fluoroquinolone',
        'colistin': 'polymyxin',
        'vancomycin': 'glycopeptide',
        'tigecycline': 'tetracycline',
        'linezolid': 'oxazolidinone'
    }
    mic_data['antibiotic_class'] = mic_data