Evolutionary trade-offs observed in antibiotic resistance can be mathematically modeled using complex interpolation of matrices to predict compensatory mutation pathways.

In a defined, contained Escherichia coli K-12 experimental-evolution system, compensatory evolution after acquisition of antibiotic resistance can be represented as transitions among genotype–phenotype states, and a complex-spectral interpolation model of environment-specific transition matrices will predict the next compensatory gene or pathway in previously unseen lineages more accurately than empirical-frequency, fitness-only, real-valued linear-interpolation, generalized linear, and sequence-based machine-learning baselines.

The primary falsifiable comparison is prospective pathway-level prediction. Before each new evolutionary interval is observed, the model must assign probabilities to the compensatory pathway that subsequently emerges. The prespecified effect is:

• At least 0.60 top-5 recall for the next validated compensatory pathway.
• At least a 20% reduction in multiclass log loss relative to the best prespecified baseline.
• At least a 10% reduction in log loss relative to a parameter-matched real-valued interpolation ablation.
• Calibration error no greater than 0.10.
• Reproduction of the effect in at least two antibiotic classes and at least three independent primary-resistance backgrounds per class.

A “compensatory mutation” is defined operationally as a secondary mutation or mutation module that, relative to its immediate resistant ancestor:

1. Increases antibiotic-free growth rate or competitive fitness by at least 10%;
2. Retains at least 80% of the ancestor’s resistance phenotype, measured on the same standardized susceptibility scale; and
3. Is supported causally by reconstruction in a nonpathogenic background, reversion, or an equivalently strong isogenic comparison.

The matrix claim is not satisfied merely by fitting retrospective observations. Predictions must be generated from a frozen model before the corresponding held-out lineage interval or prospective experiment is decoded.

The minimum viable experiment combines retrospective model development with a prospectively frozen experimental-evolution test.

A. Biological system

• Organism: a nonpathogenic E. coli K-12 reference background.
• Antibiotic coverage: three mechanistically distinct antibiotic classes selected after institutional biosafety review and baseline susceptibility testing.
• Primary-resistance backgrounds: four independently derived or archived resistant founders per antibiotic class, for 12 founders total.
• Exposure environments:
  1. Constant low selective pressure;
  2. Constant moderate selective pressure;
  3. Predefined fluctuating pressure.
• Replication: five independent lineages per founder–environment combination.
• Prospective evolving lineages: 12 founders × 3 environments × 5 replicates = 180 lineages.
• Controls:
  • 15 no-antibiotic control lineages;
  • 15 resistant-founder controls maintained without antibiotic;
  • Frozen immediate-ancestor samples for every sequenced interval;
  • Technical sequencing duplicates for 10% of samples;
  • Blinded phenotype repeats for 10% of isolates.

Exact culture and handling parameters should be established by the approved institutional SOP rather than optimized for maximal resistance. The design requires consistent transfer intervals and approximately 30–50 temporally ordered observations per lineage where feasible.

B. Sampling

• Measure bulk growth at every transfer or observation interval.
• Archive populations at least every fifth transfer interval.
• Sequence population samples at approximately 8–10 time points per lineage.
• Isolate and sequence representative clones at the resistant-founder stage, at candidate compensatory sweeps, and at the endpoint.
• Measure susceptibility and drug-free growth for all founders, endpoints, and candidate transition isolates.
• Use long-read sequencing on all founders and at least one endpoint per founder–environment combination to resolve structural variants.

Expected minimum sequencing volume:

• 180 lineages × 9 population time points = 1,620 population genomes;
• Approximately 300 clonal genomes;
• Approximately 50 long-read genomes;
• Total: approximately 1,970 sequencing libraries.

C. Retrospective training data

Before prospective decoding, construct a training corpus from:

• Qualifying serial-evolution datasets in NCBI Sequence Read Archive and European Nucleotide Archive;
• Publicly deposited E. coli resistance-evolution experiments with ordered samples and recorded antibiotic environments;
• Internal pilot lineages, if available, kept strictly separate from the prospective test;
• Resistance annotation from CARD and MEGARes;
• Gene and pathway annotation from EcoCyc, UniProt, Gene Ontology, and BV-BRC.

Only datasets containing an identifiable ancestor, at least three ordered time points, phenotype or exposure metadata, and sufficient sequencing depth are eligible.

D. Prediction freeze

After the prospective experiment is started but before later samples are sequenced:

1. Freeze the genotype-state vocabulary and pathway mapping.
2. Freeze matrix-estimation rules and hyperparameters.
3. Train the model only on retrospective data and already decoded pilot intervals.
4. Hash and timestamp the model container, configuration, pathway dictionary, and prediction files.
5. Generate ranked pathway probabilities for each blinded lineage interval.
6. Decode subsequent samples only after predictions are locked.

E. Primary analysis unit

The primary unit is one eligible ancestor-to-descendant interval. Intervals descending from the same founder remain in the same split. The outcome is the first compensatory pathway crossing a prespecified detection threshold. Intervals with no compensatory event are retained as a “none detected” class or handled by a competing-risk survival formulation.

F. Causal validation

Select:

• The ten highest-confidence predicted pathways that subsequently occur;
• Five high-confidence predictions that do not occur;
• Five unpredicted recurrent pathways.

Test these in matched, nonpathogenic resistant backgrounds using institutionally approved isogenic reconstruction, reversion, or equivalent causal comparison. Each comparison should include at least three independent biological replicates and blinded phenotype measurement.

SYSTEM NAME
    COMPASS-AMR:
    Complex Operator Modeling for Prediction of Antibiotic
    Secondary and Suppressor trajectories

HIGH-LEVEL COMPONENTS
    1. Immutable data registry
    2. Sequencing and phenotype ingestion
    3. Variant-trajectory engine
    4. Gene/pathway ontology mapper
    5. Transition-state constructor
    6. Complex spectral interpolation model
    7. Baseline-model suite
    8. Prospective prediction freezer
    9. Evaluation and statistical inference
   10. Wet-laboratory validation tracker
   11. Reproducibility and audit layer


A. PROJECT SETUP

function initialize_project(config):
    assert config.reference_genome_version is fixed
    assert config.CARD_version is fixed
    assert config.EcoCyc_version is fixed
    assert config.primary_metric == "prospective_multiclass_log_loss"
    assert config.split_unit == "resistant_founder"
    assert config.bootstrap_replicates >= 2000

    create_read_only_directory("/registry/raw")
    create_versioned_directory("/registry/processed")
    create_versioned_directory("/models")
    create_append_only_directory("/predictions_frozen")
    create_append_only_directory("/audit")

    write_preregistration_hash(config)
    build_containers([
        "read_qc",
        "variant_calling",
        "structural_variant_calling",
        "phenotype_processing",
        "matrix_model",
        "baseline_models",
        "evaluation"
    ])

    validate_biosafety_approval(config.strain_ids)
    validate_no_clinical_enhancement_work(config.protocol)


B. DATA MODEL

SampleRecord:
    sample_id
    lineage_id
    founder_id
    parent_sample_id
    antibiotic_class
    environment_id
    observation_index
    collection_timestamp
    sequencing_batch
    phenotype_batch
    blinded_status
    raw_read_uris
    checksums

PhenotypeRecord:
    sample_id
    drug_free_growth_rate
    lag_time
    carrying_capacity_proxy
    susceptibility_measure
    replicate_count
    measurement_uncertainty

VariantRecord:
    sample_id
    chromosome_or_contig
    normalized_position
    reference_allele
    alternate_allele
    variant_type
    allele_frequency_mean
    allele_frequency_interval
    read_depth
    quality_flags

TransitionRecord:
    lineage_id
    founder_id
    source_observation
    target_observation
    source_state
    target_state
    candidate_event_set
    environmental_vector
    compensation_label
    interval_censored
    observation_weight


C. INGESTION AND QUALITY CONTROL

function ingest_sample(sample_manifest):
    for sample in sample_manifest:
        verify_checksum(sample.raw_reads)
        verify_sample_id_unique(sample.sample_id)
        verify_parent_link(sample.parent_sample_id)
        register_sample(sample)

function run_sequence_pipeline(sample):
    qc = quality_control_reads(sample.raw_reads)

    if qc.contamination_score > CONFIG.max_contamination:
        quarantine(sample, reason="possible contamination")
        return

    aligned_reads = map_to_frozen_reference(
        qc.cleaned_reads,
        reference=CONFIG.reference_genome
    )

    small_variants = call_small_variants(aligned_reads)
    copy_number = estimate_copy_number(aligned_reads)
    structural_variants = call_structural_variants(aligned_reads)

    normalized = normalize_variants(
        small_variants,
        copy_number,
        structural_variants
    )

    save_with_provenance(normalized)

function process_phenotypes(raw_phenotype_table):
    randomization_check(raw_phenotype_table)
    estimate_plate_and_batch_effects()
    calculate_growth_parameters()
    calculate_susceptibility_with_uncertainty()

    if technical_CV > 0.10:
        flag_batch_for_repeat()

    return standardized_phenotypes


D. TRAJECTORY INFERENCE

function infer_lineage_trajectory(lineage_samples):
    sort lineage_samples by observation_index

    for variant in union_of_variants(lineage_samples):
        posterior_trajectory = estimate_frequency_posterior(
            variant_counts_over_time,
            read_depth_over_time,
            sequencing_error_model
        )

        establishment_interval = first_interval_where(
            posterior_probability(
                frequency > CONFIG.establishment_threshold
            ) >= CONFIG.posterior_threshold
        )

        if multiple_events_share_unresolved_interval:
            emit_set_valued_event(
                variants=events,
                interval_censored=True
            )
        else:
            emit_ordered_event(variant, establishment_interval)

    return ancestry_consistent_event_graph


E. PATHWAY AND STATE CONSTRUCTION

function map_variant_to_modules(variant):
    gene = annotate_against_reference(variant)
    resistance_labels = query_CARD_and_MEGARes(gene, variant)
    ecocyc_modules = query_EcoCyc(gene)
    go_terms = query_GO(gene)
    protein_features = query_UniProt(gene)

    return ModuleMapping(
        gene=gene,
        operons=ecocyc_modules.operons,
        pathways=ecocyc_modules.pathways,
        resistance_terms=resistance_labels,
        functional_terms=go_terms,
        provenance=all_database_versions
    )

function build_state(sample, ancestor, environment):
    resistance_modules = modules_with_established_variants(sample)
    resistance_bin = discretize_resistance(sample.susceptibility)
    cost_bin = discretize(
        ancestor.drug_free_growth - sample.drug_free_growth
    )

    return State(
        active_modules=resistance_modules,
        resistance_bin=resistance_bin,
        fitness_cost_bin=cost_bin,
        environment_vector=encode_environment(environment)
    )

function classify_compensation(source, target):
    fitness_change = (
        target.drug_free_growth - source.drug_free_growth
    ) / source.drug_free_growth

    resistance_retention = (
        target.resistance_measure / source.resistance_measure
    )

    if fitness_change >= 0.10 and resistance_retention >= 0.80:
        return "candidate_compensatory"
    elif resistance_retention < 0.80:
        return "reversion_or_resensitization"
    else:
        return "noncompensatory_or_uncertain"


F. TRANSITION MATRIX ESTIMATION

function construct_transition_counts(records, state_dictionary):
    for environment e:
        C_e = zero_matrix(number_of_states)

        for record in records where record.environment == e:
            i = state_dictionary.index(record.source_state)
            j = state_dictionary.index(record.target_state)

            weight = record.observation_weight

            if record.interval_censored:
                distribute weight across candidate target states
                according to uncertainty model
            else:
                C_e[i, j] += weight

        C_e = hierarchical_dirichlet_smooth(
            C_e,
            pathway_parent_counts,
            alpha=hyperparameter_selected_inside_training_only_CV
        )

        T_e = row_normalize(C_e)
        assert approximately_row_stochastic(T_e)

        save(T_e, uncertainty=cluster_bootstrap(records))

    return environment_transition_operators


G. TRADE-OFF OPERATOR

function construct_tradeoff_matrix(records):
    for transition i -> j:
        delta_fitness = standardized_change_in_drug_free_growth
        retained_resistance = standardized_resistance_retention
        exposure_growth = standardized_growth_under_exposure

        F[i, j] = weighted_combination(
            delta_fitness,
            retained_resistance,
            exposure_growth,
            weights_learned_on_training_data_only
        )

    return F

function augment_operator(T, F):
    # Keep transition probabilities and phenotype effects distinct
    # while allowing the interpolation model to use both.
    A = block_matrix([
        [regularize(T),        CONFIG.lambda_tradeoff * F],
        [zero_matrix_like(T),  identity_matrix(size(T))]
    ])

    return A


H. COMPLEX SPECTRAL INTERPOLATION

function stable_complex_representation(A):
    # Add only the regularization selected during inner CV.
    A_reg = regularize_for_condition_number(A)

    U, S = complex_schur_decomposition(A_reg)

    eigenmodes = extract_conjugate_paired_modes(U, S)

    if condition_number(U) > CONFIG.max_condition_number:
        return failure("unstable spectral basis")

    return SpectralRepresentation(U, S, eigenmodes)

function align_modes(rep_a, rep_b):
    cost_matrix = matrix()

    for mode p in rep_a.modes:
        for mode q in rep_b.modes:
            eigenvalue_cost = abs(p.eigenvalue - q.eigenvalue)
            vector_cost = 1 - abs(inner_product(p.vector, q.vector))
            conjugacy_penalty = incompatible_pair_penalty(p, q)

            cost_matrix[p, q] = (
                w1 * eigenvalue_cost +
                w2 * vector_cost +
                w3 * conjugacy_penalty
            )

    assignment = hungarian_algorithm(cost_matrix)
    return reorder_and_phase_align(rep_b, assignment)

function choose_log_branch(representations, training_transitions):
    candidates = enumerate_prespecified_branches(
        max_phase_wrap=CONFIG.max_phase_wrap
    )

    for branch in candidates:
        score = training_negative_log_likelihood(branch)
        instability = perturbation_instability(branch)
        imaginary_residual = reconstruction_imaginary_norm(branch)

        objective[branch] = (
            score +
            CONFIG.stability_penalty * instability +
            CONFIG.imaginary_penalty * imaginary_residual
        )

    return argmin(objective)

function interpolate_operator(anchor_operators, target_environment):
    reps = [
        stable_complex_representation(A)
        for A in anchor_operators
    ]

    reps = jointly_align_all_modes(reps)
    branch = choose_log_branch(reps, training_transitions)

    weights = environment_interpolation_weights(
        anchor_environment_vectors,
        target_environment,
        method=CONFIG.environment_kernel
    )

    interpolated_log_magnitudes = weighted_sum(
        weights,
        [log_mode_magnitudes(rep, branch) for rep in reps]
    )

    interpolated_phases = weighted_sum(
        weights,
        [unwrapped_mode_phases(rep, branch) for rep in reps]
    )

    interpolated_subspaces = unitary_geodesic_interpolation(
        [rep.invariant_subspaces for rep in reps],
        weights
    )

    A_complex = reconstruct_from_modes(
        interpolated_log_magnitudes,
        interpolated_phases,
        interpolated_subspaces
    )

    A_reconstructed = matrix_exponential(A_complex)

    imaginary_residual = norm(imaginary_part(A_reconstructed))
    A_real = real_part_after_conjugate_pair_reconstruction(A_reconstructed)

    T_raw = extract_transition_block(A_real)
    T_projected = nearest_row_stochastic_matrix(T_raw)

    projection_distance = normalized_frobenius_distance(
        T_raw,
        T_projected
    )

    diagnostics = {
        "imaginary_residual": imaginary_residual,
        "projection_distance": projection_distance,
        "condition_number": maximum_condition_number(reps),
        "branch": branch
    }

    if projection_distance > CONFIG.abort_projection_distance:
        mark_prediction_unreliable()

    return T_projected, diagnostics


I. PREDICTION

function predict_next_pathway(lineage_state, target_environment, model):
    T, diagnostics = model.interpolate_operator(target_environment)

    current_index = state_dictionary.index(lineage_state)
    state_probabilities = T[current_index, :]

    pathway_probabilities = aggregate_state_probabilities_to_pathways(
        state_probabilities,
        ontology=CONFIG.frozen_pathway_ontology
    )

    no_event_probability = competing_risk_no_event_model(
        lineage_state,
        target_environment
    )

    full_distribution = normalize(
        pathway_probabilities + no_event_probability
    )

    return Prediction(
        ranked_pathways=sort_descending(full_distribution),
        diagnostics=diagnostics
    )


J. BASELINE CONTROLS

function train_all_baselines(training_data):
    models = {}

    models["global_frequency"] = fit_global_frequency(training_data)
    models["class_frequency"] = fit_class_frequency(training_data)
    models["empirical_markov"] = fit_empirical_markov(training_data)
    models["fitness_only"] = fit_multinomial_fitness_model(training_data)
    models["real_linear"] = fit_real_convex_interpolation(training_data)
    models["real_log_euclidean"] = fit_real_log_interpolation(training_data)
    models["multinomial_glm"] = fit_regularized_multinomial(training_data)
    models["boosted_trees"] = fit_grouped_gradient_boosting(training_data)

    if effective_training_events >= CONFIG.minimum_neural_sample:
        models["pathway_neural"] = fit_pathway_feature_network(
            training_data
        )

    return models


K. SPLITTING AND LEAKAGE CONTROL

function create_grouped_splits(records):
    groups = unique(records.founder_id)

    outer_splits = grouped_cross_validation(
        groups=groups,
        stratify_by=[antibiotic_class, environment]
    )

    for split in outer_splits:
        assert no_founder_overlap(split.train, split.test)
        assert no_descendant_overlap(split.train, split.test)
        assert all_hyperparameters_fit_only_on(split.train)

    return outer_splits

function prospective_freeze(predictions, config, model):
    artifact = {
        "predictions": predictions,
        "model_hash": sha256(model),
        "config_hash": sha256(config),
        "ontology_hash": sha256(CONFIG.pathway_ontology),
        "timestamp": trusted_timestamp(),
        "allowed_evaluation_code_hash": sha256(evaluation_container)
    }

    write_once("/predictions_frozen/" + artifact.timestamp, artifact)
    lock_permissions()
    return artifact


L. EVALUATION

function evaluate_predictions(frozen_predictions, decoded_outcomes):
    verify_prediction_timestamp_precedes_outcome_decoding()

    metrics = {
        "log_loss": multiclass_log_loss(),
        "brier_score": multiclass_brier_score(),
        "top1_accuracy": top_k_recall(k=1),
        "top3_recall": top_k_recall(k=3),
        "top5_recall": top_k_recall(k=5),
        "mean_reciprocal_rank": MRR(),
        "expected_calibration_error": ECE(bins=10),
        "calibration_slope": calibration_slope()
    }

    confidence_intervals = founder_clustered_bootstrap(
        metrics,
        replicates=2000
    )

    paired_differences = compare_against_each_baseline(
        metric="log_loss",
        bootstrap_clusters="founder_id"
    )

    return metrics, confidence_intervals, paired_differences


M. ABLATION AND ADVERSARIAL TESTS

function run_ablation_suite():
    configurations = [
        "full_complex_model",
        "remove_complex_phase",
        "real_only_interpolation",
        "remove_tradeoff_matrix",
        "remove_environment_covariates",
        "remove_module_pooling",
        "remove_structural_variants",
        "remove_public_training_data"
    ]

    for configuration in configurations:
        retrain_inside_identical_grouped_splits()
        evaluate_on_same_frozen_test_intervals()
        store_paired_metric_difference()

function run_negative_controls():
    for iteration in 1..500:
        permute_outcomes_within_antibiotic_class()
        train_and_evaluate_complete_pipeline()

    simulate_real_spectrum_markov_systems()
    simulate_complex_conjugate_mode_systems()
    simulate_history_dominated_unique_pathways()

    assert model_does_not_claim_complex_advantage_on_real_only_systems()
    estimate_power_to_detect_known_complex_advantage()


N. CAUSAL VALIDATION TRACKER

function select_candidates_for_validation(predictions, outcomes):
    selected = stratified_select(
        10 correctly_predicted_high_confidence,
        5 predicted_but_not_observed,
        5 recurrent_unpredicted
    )

    return selected

function register_validation_result(candidate):
    require at_least_3_biological_replicates
    require blinded_labels
    require matched_resistant_ancestor
    require randomized_measurement_positions

    fitness_restoration = estimate_mixed_effect(
        outcome="drug_free_growth",
        fixed_effect="genotype",
        random_effect="batch"
    )

    resistance_retention = calculate_retention_ratio()

    candidate.causally_validated = (
        fitness_restoration.relative_effect >= 0.10 and
        resistance_retention >= 0.80 and
        confidence_interval_supports_effect
    )


O. OUTPUTS

Final output directory contains:
    1. Frozen prediction distributions
    2. Prospective metric table and confidence intervals
    3. Baseline and ablation comparison table
    4. Calibration plots
    5. Founder-grouped and leave-class-out analyses
    6. Spectral-mode stability report
    7. Matrix projection and branch diagnostics
    8. Causal phenotype-validation report
    9. Complete inclusion/exclusion manifest
   10. Containers, environment files, and random seeds
   11. Independent reproduction report
   12. Go/no-go decision against preregistered criteria