Evolutionary trade-offs in antibiotic resistance can be modeled and predicted using computational equilibrium concepts from game theory.

Computational game-theoretic equilibrium models (specifically Nash Equilibrium and Evolutionary Stable Strategy frameworks) can predict the frequency and magnitude of fitness trade-offs in antibiotic-resistant bacterial populations with statistically significant accuracy (R² ≥ 0.70, p < 0.05) compared to empirically measured resistance-fitness cost data across ≥3 antibiotic classes and ≥5 bacterial species, outperforming null models (random drift, simple logistic growth) by a margin of ≥15% in predictive accuracy.

PHASE 1 — Model Construction (Weeks 1–4): Build payoff matrices for 3 bacterial species (E. coli, S. aureus, P. aeruginosa) × 3 antibiotic classes (beta-lactams, fluoroquinolones, aminoglycosides) using published MIC and growth rate data. Implement Nash Equilibrium solver and Replicator Dynamics simulation. Parameterize fitness costs from literature (n ≥ 50 published competition assays per species).

PHASE 2 — In Silico Validation (Weeks 5–8): Run replicator dynamics simulations (10,000 iterations per condition) to predict equilibrium frequencies. Compare predictions against 3 published longitudinal evolution datasets (LTEE, PATRIC database, EUCAST surveillance data). Compute R², MAPE, and Brier scores.

PHASE 3 — Prospective Experimental Validation (Weeks 9–20): Conduct 96-well plate competition assays: 3 species × 3 antibiotics × 4 resistance mechanisms × 8 replicates = 288 experimental units. Measure growth rates at 0, 24, 48, 72 hours. Sequence resistant isolates (Illumina short-read, n=96) to confirm resistance genotypes. Compare observed equilibrium frequencies at generation 50 and 100 against model predictions.

PHASE 4 — Predictive Generalization Test (Weeks 21–24): Apply model trained on E. coli/beta-lactam data to predict outcomes in K. pneumoniae/carbapenem system (held-out validation). Assess transfer accuracy.

# Evolutionary Game Theory - Antibiotic Resistance Prediction Framework
# Architecture Overview

import numpy as np
import nashpy as nash
from scipy.integrate import odeint
import pymc3 as pm
from sklearn.metrics import r2_score, mean_absolute_percentage_error

# ============================================================
# MODULE 1: PAYOFF MATRIX CONSTRUCTION
# ============================================================
class ResistancePayoffMatrix:
    """
    Constructs payoff matrices from empirical fitness cost data.
    Strategies: [Susceptible, Resistant, Compensatory_Mutant, Multi_Resistant]
    """
    def __init__(self, species, antibiotic_class):
        self.species = species
        self.antibiotic_class = antibiotic_class
        self.strategies = ['S', 'R', 'CM', 'MR']
        
    def build_from_literature(self, fitness_cost_db):
        # fitness_cost_db: DataFrame with columns [species, mechanism, growth_rate, MIC]
        # Returns: NxN payoff matrix where entry [i,j] = fitness of strategy i vs strategy j
        payoff_matrix = np.zeros((4, 4))
        for i, strat_i in enumerate(self.strategies):
            for j, strat_j in enumerate(self.strategies):
                payoff_matrix[i,j] = self._compute_payoff(
                    strat_i, strat_j, fitness_cost_db
                )
        return payoff_matrix
    
    def _compute_payoff(self, s_i, s_j, db):
        # Payoff = growth_rate_i * (1 - fitness_cost_i) * antibiotic_effect(s_i, s_j)
        # antibiotic_effect depends on MIC ratio and environmental antibiotic concentration
        base_growth = db.loc[(db.species==self.species) & 
                              (db.strategy==s_i), 'growth_rate'].mean()
        fitness_cost = db.loc[(db.species==self.species) & 
                               (db.strategy==s_i), 'fitness_cost'].mean()
        mic_ratio = self._get_mic_ratio(s_i, s_j, db)
        return base_growth * (1 - fitness_cost) * mic_ratio

# ============================================================
# MODULE 2: NASH EQUILIBRIUM SOLVER
# ============================================================
class NashEquilibriumSolver:
    def __init__(self, payoff_matrix_A, payoff_matrix_B=None):
        # For symmetric games, B = A.T
        self.A = payoff_matrix_A
        self.B = payoff_matrix_B if payoff_matrix_B is not None else payoff_matrix_A.T
        self.game = nash.Game(self.A, self.B)
    
    def find_all_equilibria(self):
        equilibria = list(self.game.support_enumeration())
        return equilibria  # List of (sigma_1, sigma_2) mixed strategy pairs
    
    def find_ess(self, equilibria):
        # Evolutionary Stable Strategy: NE that satisfies second-order stability
        ess_candidates = []
        for eq in equilibria:
            sigma = eq[0]
            if self._is_ess(sigma):
                ess_candidates.append(sigma)
        return ess_candidates
    
    def _is_ess(self, sigma):
        # ESS condition: for all mutant strategies m != sigma,
        # u(sigma, sigma) > u(m, sigma) OR
        # [u(sigma, sigma) == u(m, sigma) AND u(sigma, m) > u(m, m)]
        n_strategies = len(sigma)
        for mutant_idx in range(n_strategies):
            if sigma[mutant_idx] < 1e-6:  # Not in support
                mutant = np.zeros(n_strategies)
                mutant[mutant_idx] = 1.0
                u_sigma_sigma = sigma @ self.A @ sigma
                u_mutant_sigma = mutant @ self.A @ sigma
                if u_mutant_sigma >= u_sigma_sigma:
                    return False
        return True

# ============================================================
# MODULE 3: REPLICATOR DYNAMICS SIMULATION
# ============================================================
class ReplicatorDynamics:
    def __init__(self, payoff_matrix, n_strategies):
        self.A = payoff_matrix
        self.n = n_strategies
    
    def dynamics(self, x, t):
        # dx_i/dt = x_i * [(Ax)_i - x^T A x]
        Ax = self.A @ x
        avg_fitness = x @ Ax
        dxdt = x * (Ax - avg_fitness)
        return dxdt
    
    def simulate(self, x0, t_span, n_points=1000):
        t = np.linspace(0, t_span, n_points)
        trajectory = odeint(self.dynamics, x0, t)
        equilibrium = trajectory[-1]  # Final state
        return trajectory, equilibrium
    
    def run_ensemble(self, n_initial_conditions=100, t_span=500):
        equilibria = []
        for _ in range(n_initial_conditions):
            # Random initial frequencies (Dirichlet distribution)
            x0 = np.random.dirichlet(np.ones(self.n))
            _, eq = self.simulate(x0, t_span)
            equilibria.append(eq)
        return np.array(equilibria)

# ============================================================
# MODULE 4: BAYESIAN PARAMETER ESTIMATION
# ============================================================
class BayesianPayoffEstimator:
    def __init__(self, observed_frequencies, prior_fitness_costs):
        self.observed = observed_frequencies
        self.priors = prior_fitness_costs
    
    def fit(self, n_samples=2000, n_chains=4):
        with pm.Model() as model:
            # Priors on fitness costs from literature
            fitness_cost = pm.Normal('fitness_cost', 
                                      mu=self.priors['mean'], 
                                      sigma=self.priors['std'],
                                      shape=len(self.priors['mean']))
            # Likelihood: observed frequencies ~ Dirichlet(predicted_equilibrium)
            predicted_eq = pm.Deterministic('predicted_eq', 
                                             self._compute_equilibrium(fitness_cost))
            obs = pm.Dirichlet('obs', a=predicted_eq * 100, 
                                observed=self.observed)
            trace = pm.sample(n_samples, chains=n_chains, 
                               return_inferencedata=True)
        return trace
    
    def _compute_equilibrium(self, fitness_cost):
        # Differentiable equilibrium computation for gradient-based sampling
        # Simplified: returns softmax of fitness-adjusted payoffs
        adjusted_payoffs = pm.math.softmax(-fitness_cost)
        return adjusted_payoffs

# ============================================================
# MODULE 5: VALIDATION AND METRICS
# ============================================================
class ModelValidator:
    def __init__(self, predictions, observations):
        self.pred = predictions
        self.obs = observations
    
    def compute_metrics(self):
        metrics = {
            'r2': r2_score(self.obs, self.pred),
            'mape': mean_absolute_percentage_error(self.obs, self.pred),
            'pearson_r': np.corrcoef(self.obs, self.pred)[0,1],
            'mean_abs_error': np.mean(np.abs(self.obs - self.pred))
        }
        return metrics
    
    def permutation_test(self, n_permutations=10000):
        observed_r2 = r2_score(self.obs, self.pred)
        null_r2s = []
        for _ in range(n_permutations):
            shuffled_pred = np.random.permutation(self.pred)
            null_r2s.append(r2_score(self.obs, shuffled_pred))
        p_value = np.mean(np.array(null_r2s) >= observed_r2)
        return p_value, observed_r2, np.array(null_r2s)
    
    def compare_to_null_model(self, null_predictions):
        null_r2 = r2_score(self.obs, null_predictions)
        model_r2 = r2_score(self.obs, self.pred)
        improvement = (model_r2 - null_r2) / max(abs(null_r2), 1e-6)
        return {'null_r2': null_r2, 'model_r2': model_r2, 
                'relative_improvement': improvement}

# ============================================================
# MODULE 6: MAIN PIPELINE
# ============================================================
def run_validation_pipeline(species_list, antibiotic_list, fitness_db, 
                              experimental_data):
    results = {}
    for species in species_list:  # ['E.coli', 'S.aureus', 'P.aeruginosa']
        for antibiotic in antibiotic_list:  # ['beta-lactam', 'fluoroquinolone', 'aminoglycoside']
            condition = f"{species}_{antibiotic}"
            
            # Step 1: Build payoff matrix
            rpm = ResistancePayoffMatrix(species, antibiotic)
            A = rpm.build_from_literature(fitness_db)
            
            # Step 2: Find Nash Equilibria
            nes = NashEquilibriumSolver(A)
            equilibria = nes.find_all_equilibria()
            ess = nes.find_ess(equilibria)
            
            # Step 3: Simulate replicator dynamics
            rd = ReplicatorDynamics(A, n_strategies=4)
            trajectories = rd.run_ensemble(n_initial_conditions=100)
            predicted_eq = np.median(trajectories, axis=0)
            
            # Step 4: Bayesian refinement
            obs_freq = experimental_data[condition]['frequencies']
            priors = fitness_db.get_priors(species, antibiotic)
            bpe = BayesianPayoffEstimator(obs_freq, priors)
            trace = bpe.fit()
            
            # Step 5: Validate
            validator = ModelValidator(predicted_eq, obs_freq)
            metrics = validator.compute_metrics()
            p_val, obs_r2, null_dist = validator.permutation_test()
            
            results[condition] = {
                'nash_equilibria': equilibria,
                'ess': ess,
                'predicted_frequencies': predicted_eq,
                'metrics': metrics,
                'p_value': p_val,
                'bayesian_trace': trace
            }
    
    return results

# ============================================================
# EXECUTION
# ============================================================
# species_list = ['E.coli', 'S.aureus', 'P.aeruginosa']
# antibiotic_list = ['beta-lactam', 'fluoroquinolone', 'aminoglycoside']
# results = run_validation_pipeline(species_list, antibiotic_list, 
#                                    fitness_db, experimental_data)
# Generate summary report with all metrics and visualizations