The complete collateral sensitivity graph for WHO priority pathogens — where nodes are antibiotics and directed edges encode resistance to drug A causing increased susceptibility to drug B — contains strongly connected components that define closed evolutionary traps: resistance mutation cycles from which no clinically viable escape route exists. An algorithmic method applied to GLASS genomic data will identify these traps, and sequential therapy protocols cycling through the trap sequence will maintain pathogen susceptibility indefinitely because each resistance mutation obligately returns the organism to susceptibility to a prior drug in the cycle.

For at least one WHO priority pathogen (e.g., K. pneumoniae, A. baumannii, P. aeruginosa, E. coli), a directed collateral sensitivity graph (CSG) constructed from published and GLASS-derived genomic resistance data contains at least one strongly connected component (SCC) of size ≥ 3 antibiotics such that: (1) every directed edge A→B represents a statistically significant (p < 0.05, OR > 2.0) increase in susceptibility to drug B conditioned on resistance to drug A; (2) the SCC constitutes a "closed evolutionary trap" meaning no resistance mutation accessible from any node in the SCC leads to a state outside the SCC with clinical viability (MIC below EUCAST breakpoint); and (3) a sequential cycling protocol following the SCC order maintains pathogen MIC below clinical resistance breakpoints for ≥ 30 serial passages in vitro, compared to ≥ 50% of monotherapy controls reaching resistance by passage 15.

PHASE 1 — Computational Graph Construction (Weeks 1–8): Mine GLASS WGS data + 47 published collateral sensitivity studies to construct species-specific CSGs. Apply SCC detection (Tarjan's algorithm). Validate edge significance with logistic regression controlling for genetic background and geographic confounders.

PHASE 2 — Evolutionary Accessibility Analysis (Weeks 6–12): For each candidate SCC, enumerate all single-nucleotide variants within 2 mutational steps of each node using DMS (deep mutational scanning) data and structural modeling. Identify whether any variant simultaneously confers resistance to all SCC drugs.

PHASE 3 — In Vitro Evolution Validation (Weeks 10–24): Conduct serial passage experiments in 96-well format using the top-ranked SCC cycling protocol vs. monotherapy and random cycling controls. Use 12 replicate lineages per condition. Perform WGS at passages 5, 10, 20, 30.

PHASE 4 — Mechanistic Confirmation (Weeks 20–30): Introduce specific resistance mutations via CRISPR-Cas9 into isogenic backgrounds to confirm directed edge causality (not correlation). Measure MIC changes for all SCC drugs after each engineered mutation.

PHASE 5 — Clinical Data Validation (Weeks 24–36): Retrospective analysis of GLASS patient isolate pairs (admission vs. discharge) to test whether patients treated with drugs approximating SCC cycling show lower resistance emergence rates than those on monotherapy.

# COLLATERAL SENSITIVITY GRAPH (CSG) EVOLUTIONARY TRAP FINDER
# Architecture: 5-module pipeline

# ============================================================
# MODULE 1: DATA INGESTION & STANDARDIZATION
# ============================================================
class AMRDataIngestion:
    def load_glass_wgs(self, path: str) -> pd.DataFrame:
        # Load GLASS assemblies + phenotype table
        # Columns: isolate_id, species, country, year, drug_1..drug_N (S/I/R)
        # Run AMRFinder+ on all assemblies -> genotype matrix
        # Standardize to EUCAST 2024 breakpoints
        return standardized_phenotype_matrix  # shape: [n_isolates x n_drugs]

    def load_published_cs_data(self, meta_analysis_db: str) -> pd.DataFrame:
        # Parse 47 published collateral sensitivity studies
        # Standardize drug names to ATC codes
        # Output: [study_id, drug_A, drug_B, OR, CI_lower, CI_upper, species]
        return cs_compendium  # shape: [n_pairs x 7]

# ============================================================
# MODULE 2: EDGE ESTIMATION (LOGISTIC REGRESSION)
# ============================================================
class EdgeEstimator:
    def compute_cs_edges(self, phenotype_matrix, genotype_matrix):
        edges = {}
        drugs = phenotype_matrix.columns
        # Compute population structure PCs
        pcs = PCA(n_components=10).fit_transform(genotype_matrix)
        
        for drug_A in drugs:
            for drug_B in drugs:
                if drug_A == drug_B:
                    continue
                # Subset to isolates resistant to drug_A
                mask_resistant_A = phenotype_matrix[drug_A] == 'R'
                # Logistic regression: P(susceptible_B) ~ resistant_A + PCs + country + year
                X = np.column_stack([
                    (phenotype_matrix[drug_A] == 'R').astype(int),
                    pcs, country_dummies, year_covariate
                ])
                y = (phenotype_matrix[drug_B] == 'S').astype(int)
                model = LogisticRegression().fit(X, y)
                OR = np.exp(model.coef_[0][0])
                p_value = compute_wald_pvalue(model)
                edges[(drug_A, drug_B)] = {'OR': OR, 'p': p_value}
        
        # Apply Benjamini-Hochberg FDR correction
        edges_df = pd.DataFrame(edges).T
        edges_df['p_adj'] = multipletests(edges_df['p'], method='fdr_bh')[1]
        # Filter: OR > 2.0, p_adj < 0.05
        significant_edges = edges_df[(edges_df['OR'] > 2.0) & (edges_df['p_adj'] < 0.05)]
        return significant_edges

# ============================================================
# MODULE 3: GRAPH CONSTRUCTION & SCC DETECTION
# ============================================================
class CSGraphAnalyzer:
    def build_graph(self, significant_edges: pd.DataFrame) -> nx.DiGraph:
        G = nx.DiGraph()
        for (drug_A, drug_B), row in significant_edges.iterrows():
            G.add_edge(drug_A, drug_B, weight=row['OR'], p_adj=row['p_adj'])
        return G

    def find_evolutionary_traps(self, G: nx.DiGraph) -> list:
        # Tarjan's algorithm via NetworkX
        sccs = list(nx.strongly_connected_components(G))
        # Filter: size >= 3, all edges significant
        traps = []
        for scc in sccs:
            if len(scc) >= 3:
                subgraph = G.subgraph(scc)
                # Check all edges within SCC are significant
                min_OR = min([d['weight'] for _, _, d in subgraph.edges(data=True)])
                if min_OR > 2.0:
                    traps.append({
                        'nodes': list(scc),
                        'size': len(scc),
                        'mean_OR': np.mean([d['weight'] for _, _, d in subgraph.edges(data=True)]),
                        'subgraph': subgraph
                    })
        # Rank by mean OR * size
        traps.sort(key=lambda x: x['mean_OR'] * x['size'], reverse=True)
        return traps

    def null_model_test(self, G: nx.DiGraph, n_permutations: int = 10000) -> float:
        observed_max_scc = max(len(s) for s in nx.strongly_connected_components(G))
        null_distribution = []
        for _ in range(n_permutations):
            G_rand = nx.directed_configuration_model(
                [d for _, d in G.in_degree()],
                [d for _, d in G.out_degree()]
            )
            null_max = max(len(s) for s in nx.strongly_connected_components(G_rand))
            null_distribution.append(null_max)
        empirical_p = np.mean(np.array(null_distribution) >= observed_max_scc)
        return empirical_p

# ============================================================
# MODULE 4: ESCAPE ROUTE ANALYSIS
# ============================================================
class EscapeRouteAnalyzer:
    def enumerate_escape_mutations(self, scc_nodes: list, dms_data: dict) -> pd.DataFrame:
        # For each drug in SCC, get all resistance mutations from DMS data
        # Check if any single mutation confers resistance to ALL SCC drugs
        all_mutations = {}
        for drug in scc_nodes:
            all_mutations[drug] = dms_data[drug]['resistance_mutations']  # {mutation: MIC_fold_change}
        
        # Find mutations appearing in resistance sets for multiple drugs
        escape_candidates = []
        for mutation in set.union(*[set(m.keys()) for m in all_mutations.values()]):
            drugs_resisted = [drug for drug in scc_nodes 
                            if mutation in all_mutations[drug] 
                            and all_mutations[drug][mutation] > 2.0]  # >2x MIC increase
            if len(drugs_resisted) == len(scc_nodes):
                escape_candidates.append({'mutation': mutation, 'full_escape': True})
            elif len(drugs_resisted) >= 2:
                escape_candidates.append({'mutation': mutation, 'full_escape': False,
                                         'drugs_resisted': drugs_resisted})
        return pd.DataFrame(escape_candidates)

# ============================================================
# MODULE 5: IN VITRO CYCLING PROTOCOL OPTIMIZER
# ============================================================
class CyclingProtocolOptimizer:
    def design_cycling_protocol(self, trap: dict) -> list:
        # Find Hamiltonian cycle through SCC maximizing minimum edge OR
        # Use TSP solver (OR-Tools) with edge weights = 1/OR
        from ortools.constraint_solver import routing_enums_pb2, pywrapcp
        nodes = trap['nodes']
        G = trap['subgraph']
        # Maximize minimum OR along cycle path
        # Returns ordered drug sequence for cycling protocol
        optimal_cycle = self._solve_max_min_hamiltonian(G, nodes)
        return optimal_cycle

    def simulate_evolution(self, cycling_protocol: list, 
                          n_passages: int = 30,
                          n_replicates: int = 12) -> pd.DataFrame:
        # Stochastic simulation of resistance evolution under cycling
        # Uses Wright-Fisher model with mutation rates from literature
        # Returns MIC trajectories for all replicates
        results = []
        for rep in range(n_replicates):
            mic_trajectory = self._wright_fisher_simulation(
                cycling_protocol, n_passages,
                mutation_rate=1e-9, population_size=1e8
            )
            results.append(mic_trajectory)
        return pd.DataFrame(results)

# ============================================================
# MAIN PIPELINE
# ============================================================
def run_evp_pipeline(glass_path, published_cs_path, output_dir):
    # Step 1: Load data
    ingestion = AMRDataIngestion()
    phenotype_matrix = ingestion.load_glass_wgs(glass_path)
    cs_compendium = ingestion.load_published_cs_data(published_cs_path)
    
    # Step 2: Estimate edges
    estimator = EdgeEstimator()
    significant_edges = estimator.compute_cs_edges(phenotype_matrix, genotype_matrix)
    
    # Step 3: Build graph and find traps
    analyzer = CSGraphAnalyzer()
    G = analyzer.build_graph(significant_edges)
    traps = analyzer.find_evolutionary_traps(G)
    null_p = analyzer.null_model_test(G)
    
    # Step 4: Escape analysis
    escape_analyzer = EscapeRouteAnalyzer()
    for trap in traps[:3]:  # Top 3 traps
        escapes = escape_analyzer.enumerate_escape_mutations(trap['nodes'], dms_data)
        trap['escape_routes'] = escapes
    
    # Step 5: Design cycling protocols
    optimizer = CyclingProtocolOptimizer()
    for trap in traps[:3]:
        trap['cycling_protocol'] = optimizer.design_cycling_protocol(trap)
        trap['simulation'] = optimizer.simulate_evolution(trap['cycling_protocol'])
    
    # Output: ranked trap list, cycling protocols, escape analysis, null model p-value
    save_results(traps, null_p, output_dir)
    return traps, null_p

# Expected outputs:
# - CSG visualization (Gephi-compatible .gexf)
# - SCC report with OR, p-values, clinical drug availability
# - Escape route analysis table
# - Cycling protocol sequences (ordered drug lists)
# - Simulation MIC trajectories
# - Statistical validation (null model p-values)