Engineering

Iteration 1 — Build Phase

Mutations were generated through orthogonal DNA replication using EcORep's linear O-replicon. To analyze whether our luminescence reporter could quantitatively reflect protein performance under this biological complexity, we constructed an in vivo ODE model [View Code in Google Colab].

This ODE described:

• Logistic cell growth \( \dfrac{dE}{dt} = rE\!\left(1-\dfrac{E}{K}\right) \)
• Expression and dilution of SpyCatcher and SpyTag proteins²
• Their binding into a luminescent complex \(AB\) with rate \(k_{\text{bind}}\)
• Luminescence activation \( \dfrac{dL}{dt} = k_{\text{act}}(AB) - k_{\text{decay}}\,L \)

By simulating how growth, expression, and binding co-evolve, we could test whether our cellular system would allow distinct k_on values to be distinguished experimentally.

Iteration 1 — Test Phase

We ran parameter scans varying r (bacterial growth rate constant), α_C (synthesis rate of SpyCatcher), α_T (synthesis rate of SpyTag), k_bind (binding rate constant), and K (carrying capacity).

The results revealed that even when k_bind differed by an order of magnitude, the luminescence difference remained minor due to dilution and growth effects. This meant our in vivo system could not resolve kinetic differences between mutants.

Figure 1. In vivo ODE simulation of luciferase luminescence dynamics for two kinetic regimes (k = 60.00 vs 20.00 µM⁻¹·h⁻¹).

Iteration 1 — Learn Phase

Although the in vivo ODE confirmed stable mutation propagation, it also showed that the system was unidentifiable for kinetic parameters — luminescence reflected cell physiology more than binding kinetics.

Furthermore, it was impossible to evaluate the rate constant under different pH conditions inside cells.

Thus, the model directly informed our decision to decouple mutation from kinetic validation, motivating a move toward in vitro testing.

Iteration 2 — Design Phase

To bridge biological relevance with measurable precision, we designed a hybrid cycle: mutations still occurred in vivo, but validation was shifted in vitro.

Iteration 2 — Build Phase

During this stage, we continued using the EcORep orthogonal replication system to generate mutations inside E. coli.

However, sequence analysis revealed that the in vivo mutation process was too extensive — the entire O-replicon, including the luciferase reporter³ and signal peptide regions⁴, accumulated random mutations rather than remaining confined to the intended target gene.

We therefore implemented in silico structure prediction (RF Diffusion, AlphaFold2) to investigate whether such global mutations could compromise protein folding and complex formation.

Iteration 2 — Test Phase

Simulation results showed that these uncontrolled mutations frequently disrupted key structural motifs, leading to a collapse of the luciferase and signal peptide architecture.

This indicated that the in vivo mutagenic environment itself was unstable and that excessive replication errors could compromise the integrity of our reporter system.

Iteration 2 — Learn Phase

From this, we learned that performing mutation in vivo risks unwanted propagation across the entire O-replicon, damaging both expression and functional modules.

To ensure precise control, we decided to move mutation generation in vitro using error-prone PCR,⁵⁶⁷⁸⁹ where the mutational rate and target region could be strictly limited. This insight set the foundation for our fully in vitro epPCR-based workflow in the next cycle.

Iteration 3 — Design Phase

After recognizing the limits of cellular context, we transitioned to a fully in vitro workflow. Both mutation generation and validation were moved outside the cell for full control. Mutants were now created via error-prone PCR, ensuring adjustable mutation rates under defined Mn²⁺ and dNTP conditions.

Iteration 3 — Build Phase

To design the validation workflow, we constructed an in vitro kinetics ODE [View Code in Google Colab] describing the reaction:

A (SpyTag) + B (SpyCatcher) → AB

as an irreversible second-order process:

and mapped AB(t) to plate-reader luminescence through a saturating transfer function and activation/decay kinetics:

By scanning k_bind (≈ 6 vs 72 µM⁻¹ h⁻¹), we could determine the optimal starting concentrations and time points to maximize signal discrimination.

Iteration 3 — Test Phase

Simulations predicted that using 0.1 µM of each reactant and measuring at 5, 15, 30, 60 min would clearly separate mutants by k_on.

This ODE-guided design minimized trial-and-error in wet-lab testing.

Figure 2. In vitro simulation of split-luciferase luminescence under different binding rate constants (k₁ = 60 µM⁻¹·h⁻¹, k₂ = 20 µM⁻¹·h⁻¹).

Iteration 3 — Learn Phase

This cycle validated that a fully in vitro strategy yields quantitative and reproducible kinetic data under tunable pH.

The modeling framework now serves as a pre-experiment planning tool: each new mutant library first undergoes ODE simulation to set experimental parameters.

Iteration 1 — Build Phase

We model base changes with the Jukes–Cantor (JC69) 4-state substitution model (equal base frequencies/rates). Under JC69, the probability a nucleotide remains the same after c cycles at per-cycle rate μ is:

and P_diff(c) = 1 - P_same(c). Matching a measured per-reaction error p to P_diff(c) and solving for μ yields:

We use Eq. 2 to calibrate μ for any chosen p and cycle count c.³

# --- Code A: JC69 essentials (calibration + transition probabilities) ---
import numpy as np

def per_reaction_to_per_cycle_error(p_fraction: float, cycles: int) -> float:
    """
    Convert a per-reaction error fraction p (e.g., 0.02 for 2%) into
    the per-cycle substitution rate μ using JC69. 
    """
    if cycles <= 0:
        return 0.0
    return 0.75 * (1.0 - (1.0 - (4.0 * p_fraction) / 3.0) ** (1.0 / cycles))

def jc69_psame_pdiff(mu: float, cycles: int):
    """
    JC69 transition probabilities after 'cycles' steps: 
      P_same(c) = 1/4 + 3/4 * (1 - 4μ/3)^c
      P_diff(c) = 1 - P_same(c)
    """
    p_same = 0.25 + 0.75 * (1.0 - 4.0 * mu / 3.0) ** cycles
    return p_same, 1.0 - p_same

def jc69_transition_matrix(p_same: float):
    """
    4×4 single-base transition matrix with equal off-diagonals (A,C,G,T). 
    """
    P = np.full((4, 4), (1.0 - p_same) / 3.0)
    np.fill_diagonal(P, p_same)
    bases = ['A', 'C', 'G', 'T']
    return P, bases

Iteration 1 — Test Phase

We validated Eq. 2 numerically: plugging in p = {0.5, 1, 2, 3.5}% and c = 10, 20, 40 gave self-consistent P_diff(c) when re-composed via Eq. 1.

Iteration 1 — Learn Phase

We confirmed the need to separate per-reaction p from per-cycle μ; without this calibration, downstream stop-codon risk is badly mis-estimated at high cycle counts. We proceed to codon-level predictions.

Iteration 2 — Design Phase

Goal: compute the probability an ORF accumulates new stop codons after epPCR. Literature explains that mutagenic conditions increase substitution frequency broadly (not only transitions), consistent with a symmetric model at planning time.¹

We define two usability metrics used throughout:

• Usability-A (no-new-stops): fraction of molecules that introduce no additional stop codon relative to the original ORF (the strict screening criterion).
• Usability-B (final-no-stop): fraction of molecules that end stop-free (informative but less strict).

Iteration 2 — Build Phase

For each codon b₁b₂b₃, we compute the probability it becomes any of {TAA, TAG, TGA} after c cycles using the single-base transition probabilities from JC69 (independence across positions):

We then multiply the no-stop probabilities across all relevant codons (log-sum for stability) to get Usability-A, and similarly compute Usability-B.³

Figure 1. Usability-A (%) surface versus per-reaction error and cycles; red points = 10%-usability ridge (design threshold).

Figure 2. Expected mutated sites; surface colored by concentration score 1/(1+CV); red points = 10% ridge.

# --- Code B: codon-level stop risk + usability definitions (uses Code A) ---
import numpy as np

STOP_CODONS = {"TAA", "TAG", "TGA"}

def clean_sequence(seq: str) -> str:
    return ''.join(b for b in seq.upper() if b in "ACGT")

def get_codons(seq: str, frame: int = 0):
    if frame:
        seq = seq[frame:]
    k = len(seq) % 3
    if k:
        seq += "N" * (3 - k)
    return [c for c in (seq[i:i+3] for i in range(0, len(seq), 3)) if "N" not in c]

def codon_to_stop_probability(codon: str, P, bases) -> float:
    """
    P(codon → any STOP) after c cycles under JC69:
        sum_{stop ∈ {TAA,TAG,TGA}} ∏_{i=1..3} P[b_i → stop_i]. 
    """
    idx = {b: i for i, b in enumerate(bases)}
    total = 0.0
    for stop in STOP_CODONS:
        p = 1.0
        for a, b in zip(codon, stop):
            p *= P[idx[a], idx[b]]
        total += p
    return total

def usability_fraction(sequence: str, P, bases, frame: int = 0, mode: str = "no_new_stops") -> float:
    """
    Usability-A ('no_new_stops'): fraction with NO added stops (strict).
    Usability-B ('final_no_stop'): fraction that ends stop-free (relaxed).
    """
    codons = get_codons(clean_sequence(sequence), frame)
    if not codons:
        return 1.0

    idx_stop0 = {i for i, c in enumerate(codons) if c in STOP_CODONS}

    if mode == "no_new_stops":
        relevant = (c for i, c in enumerate(codons) if i not in idx_stop0)
    else:  
        relevant = (c for c in codons)

    log_ok = 0.0
    for codon in relevant:
        p_stop = codon_to_stop_probability(codon, P, bases)
        p_ok = max(1e-300, 1.0 - p_stop)  # numerical safety
        log_ok += np.log(p_ok)
    return float(np.exp(log_ok))

Iteration 2 — Test Phase

Using your target ORF and frames 0/1/2:

• At 1% error, 20 cycles, Usability-A > 90%; Usability-B slightly higher (as expected).
• At 3% error, 50 cycles, Usability-A falls below 10%, indicating most variants would be unscreenable under strict criteria.
• Frame choice can subtly change per-codon risk near boundaries, visible in the heatmap of stop-probabilities.

Iteration 2 — Learn Phase

The strict metric (Usability-A) is what matters for protein-level screening; Usability-B is retained as an educational/diagnostic view. With this decision, we move to selecting cycles that keep Usability-A acceptable.

Iteration 3 — Design Phase

Goal: define a simple rule for choosing cycle counts: the ridge where Usability-A = 10% across the (p,c) grid. This aligns with practical guidance that highly mutagenic mixes (Mn²⁺, unbalanced dNTPs) should limit cycles to avoid unusable libraries.¹

Iteration 3 — Build Phase

We implemented a ridge finder that, for each p, selects the cycle c whose Usability-A is closest to 10% and exports a CSV summary.

# --- Code C: grid evaluation + 10% ridge + minimal plotting/export (uses A & B) ---
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def evaluate_grid(sequence: str,
                  p_min=0.5, p_max=4.0, p_step=0.05,   
                  c_min=1, c_max=100, c_step=1,        
                  frame=0, mode="no_new_stops"):
    rows = []
    p_vals = np.arange(p_min, p_max + 1e-12, p_step) / 100.0
    c_vals = np.arange(c_min, c_max + 1e-12, c_step)

    for p in p_vals:
        for c in c_vals:
            mu = per_reaction_to_per_cycle_error(p, int(c))
            p_same, p_diff = jc69_psame_pdiff(mu, int(c))
            P, bases = jc69_transition_matrix(p_same)

            u = usability_fraction(sequence, P, bases, frame=frame, mode=mode)
            rows.append({
                "per_reaction_error_pct": p * 100.0,
                "cycles": int(c),
                "per_cycle_mu": mu,
                "p_same": p_same,
                "p_diff": p_diff,
                "usability_pct": 100.0 * u,
                "expected_mutations": len(clean_sequence(sequence)) * p_diff
            })
    return pd.DataFrame(rows)

def find_usability_ridge(df: pd.DataFrame, target=10.0):
    """For each p, find the cycle where usability_pct is closest to 'target' (e.g., 10%)."""
    out = []
    for p in sorted(df["per_reaction_error_pct"].unique()):
        sub = df[df["per_reaction_error_pct"] == p].copy()
        k = (sub["usability_pct"] - target).abs().idxmin()
        r = sub.loc[k]
        out.append({"per_reaction_error_pct": p,
                    "cycles_at_target": int(r["cycles"]),
                    "usability_at_target": float(r["usability_pct"]),
                    "expected_mutations_at_target": float(r["expected_mutations"])})
    return pd.DataFrame(out)

def plot_usability_surface(df: pd.DataFrame, ridge: pd.DataFrame):
    p_unique = sorted(df["per_reaction_error_pct"].unique())
    c_unique = sorted(df["cycles"].unique())
    P, C = np.meshgrid(p_unique, c_unique)

    Z = np.zeros_like(P, dtype=float)
    for i, c in enumerate(c_unique):
        for j, p in enumerate(p_unique):
            v = df[(df["per_reaction_error_pct"] == p) & (df["cycles"] == c)]["usability_pct"]
            if not v.empty:
                Z[i, j] = v.iloc[0]

    fig = plt.figure(figsize=(9, 7))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(P, C, Z, cmap="viridis", alpha=0.8)
    if not ridge.empty:
        ax.scatter(ridge["per_reaction_error_pct"], ridge["cycles_at_target"],
                   ridge["usability_at_target"], color="red", s=30, label="10% ridge")
        ax.legend()
    ax.set_xlabel("Per-reaction error (%)")
    ax.set_ylabel("PCR cycles")
    ax.set_zlabel("Usability-A (%)")
    ax.set_title("Figure 2: Usability-A surface with 10% ridge")
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    TARGET_SEQUENCE = "ACGT" * 150  # replace with your ORF
    df = evaluate_grid(TARGET_SEQUENCE, p_min=0.5, p_max=4.0, p_step=0.1,
                       c_min=1, c_max=100, c_step=2, frame=0, mode="no_new_stops")
    ridge = find_usability_ridge(df, target=10.0)
    df.to_csv("pcr_error_grid.csv", index=False)
    ridge.to_csv("usability_ridge_10pct.csv", index=False)
    plot_usability_surface(df, ridge)

Iteration 3 — Test Phase

Representative results from our grid (sequence-specific but robust):

• 0.5% error → ~80–90 cycles at 10% Usability-A.
• 2.0% error → ~25 cycles.
• 3.5% error → <15 cycles.

This matches practical experience: strongly error-prone conditions demand fewer cycles to keep usable ORFs.¹

Iteration 3 — Learn Phase

Actionable rule for Build:

• For kit-like recipes near ~2% error, cap epPCR at ≈20–25 cycles.¹
• For milder recipes (~0.6–1.0%), cycles can extend ≥60 while maintaining a reasonable usability band.¹

These limits carry directly into our wet-lab rounds and sample-throughput calculations (starting copy number × amplification factor × Usability-A).

Iteration 1 — Build Phase

Each construct used a (GGGGS)₃ linker to maintain flexibility between domains.

Expression in E. coli and western blot analysis verified full-length fusion proteins, indicating structural compatibility.

Iteration 1 — Test Phase

AlphaFold2 modeling and prior structural analyses support that SpyCatcher at the N'-luciferase C terminus minimizes steric hindrance and preserves the catalytic region's architecture.³⁴

Furthermore, the exposed positioning of SpyCatcher's reactive lysine allows spontaneous covalent ligation upon SpyTag encounter, offering a structurally and kinetically favorable orientation for reconstitution.

N-terminal Split Luciferase fused with SpyCatcher

C-terminal Split Luciferase fused with SpyTag

Iteration 1 — Learn Phase

This cycle established that N'-Luciferase–SpyCatcher and C-'Luciferase–SpyTag fusions provide a stable, functionally rational design for split-luciferase reconstitution, preserving both luciferase folding and SpyCatcher reactivity.

Build Phase

Following the original design, we prepared to simultaneously introduce:

• a linear O-replicon as the mutational cargo;
• a replication plasmid encoding TP, O-DNAP, SSB, and DSB;
• an auxiliary λ-Gam plasmid to inhibit the host RecBCD nuclease;
• and an additional plasmid carrying mutagenic variants such as O-DNAP(N71D) or O-DNAP(Y127A).

This configuration required maintaining three to four plasmids with different copy numbers and antibiotic markers under IPTG induction.

Test Phase

Before attempting the full reconstruction, we consulted our instructors and molecular biology advisors to evaluate the practical feasibility of this setup.

They pointed out that:

"Although the original EcORep system is theoretically feasible, maintaining multiple plasmids with different replication origins and antibiotic resistances will impose severe metabolic burden and instability on E. coli. In addition, some modules require specific host strains and tightly controlled expression conditions to function properly."

Based on this expert feedback, we realized that rebuilding the complete EcORep architecture under our current lab conditions would not be achievable.

Attempting the design as-is would likely result in low transformation efficiency, high cellular stress, and rapid plasmid loss.

This discussion served as a theoretical validation of our design, allowing us to recognize its practical limitations before proceeding.

Learn Phase

Through this test and expert consultation, we learned that while the multi-layered EcORep system is conceptually powerful, its high complexity and plasmid burden make it difficult to reproduce and sustain in standard E. coli settings.

As a result, we decided to abandon the multi-plasmid reconstruction approach and instead redesign a simplified two-component version:

• a linear O-replicon and
• a single helper plasmid that integrates all essential replication modules.

This realization became the foundation for Cycles 4 and 5, where we began re-engineering EcORep into a more modular, stable, and accessible form.

Build Phase (Shared with Cycle 6)

A linear O-replicon was constructed, flanked by 18 bp PRD1 inverted terminal repeats (ITRs) (The paper mentioned that using 18 bp from each side works), and containing a resistance cassette and reporter gene (GFP in pretest) for selection.

The EcoRep plasmid carried both Operon 1 and Operon 2, allowing independent control by IPTG and rhamnose.

Co-transformation into E. coli established a dual-operon system capable of either faithful or mutagenic replication depending on induction conditions.

Cycle 5 — Replication Test

Cells were induced with IPTG to activate Operon 1.

Replication efficiency was assessed through kanamycin and ampicillin resistance as well as DNA retention after 48 hours of subculture.

The linear O-replicon was stably maintained with no significant loss, confirming that the three-gene module was sufficient for orthogonal DNA replication in E. coli.

Cycle 5 — Replication Learn

This test demonstrated that a minimal set of TP, O-DNAP, and SSB can sustain replication independently of the host genome.

The removal of DSB components reduced genetic load while preserving replication stability, providing a solid foundation for subsequent mutation control experiments.

Build Phase (Shared with Cycle 4)

A linear O-replicon was constructed, flanked by 18 bp PRD1 inverted terminal repeats (ITRs) (The paper mentioned that using 18 bp from each side works), and containing a resistance cassette and reporter gene (GFP in pretest) for selection.

The EcoRep plasmid carried both Operon 1 and Operon 2, allowing independent control by IPTG and rhamnose.

Co-transformation into E. coli established a dual-operon system capable of either faithful or mutagenic replication depending on induction conditions.

Cycle 6 — Mutation Test

To evaluate mutation control, samples were collected at 0 h, 24 h, 36 h, and 48 h after rhamnose induction to activate Operon 2.

The SpyCatcher region was sequenced to measure mutation accumulation over time.

Results suggested that mutation accumulation should show a positive correlation with induction time, but due to the low overall mutation rate or the relatively short induction period, the trend was not clearly distinguishable within the observed timeframe.

Cycle 6 — Mutation Learn

These findings confirmed that mutation strength can be finely tuned by adjusting rhamnose induction duration without compromising system stability.

The dual-operon architecture provides a reliable framework for controlled, gradual mutation — a crucial step toward continuous directed evolution in E. coli.