Model

Model Overview

From 0.03% to 87%: The Truth Behind the Efficiency Gap

Imagine a gene integration system where success rates can vary by nearly 3000-fold under different conditions. This isn't an experimental error—it’s a black box waiting to be opened.

In our GenOMe system, one-to-one means that one pair of attB and attP sites recombines once through Bxb1 integrase to form one integration event.

Traditional genetic engineering relies on trial and error, but we asked:

• Which step truly limits efficiency?
• Why is two-to-two deletion so much harder than one-to-one insertion?
• Can we mathematically predict optimal conditions instead of blind testing?

Our answer: Build an interpretable, predictive five-layer integration model.

The Efficiency Gap at a Glance

A single variable change produces a three-order-of-magnitude leap — this is the core mystery our model explains.

Modeling Strategy: Dissecting the Black Box Layer by Layer

We decomposed the Bxb1 gene integration system into five modelable biological layers, peeling away complexity like an onion:

Model1(M1): Transcriptional Dynamics and Induction Window Optimization

Quantifies how arabinose initiates Bxb1 gene transcription and the dynamics of mRNA production and degradation.

M2 - Translation Layer: The Protein Accumulation Race

Calculates how long it takes for Bxb1 protein to reach effective concentration and when it begins to lose activity.

M3 - Chromosomal Layer: The Hidden Major Bottleneck

Models how ssAP forms attB sites on the chromosome — the system’s most fragile step.

Mechanism Breakdown: ssAP (single-strand annealing protein) must:

1. Be taken up by cells (efficiency ~10–30%)
2. Find the replication fork spatiotemporal window
3. Successfully anneal with the lagging strand
4. Avoid degradation or dilution

The multiplication of these four steps results in overall success rate < 1%.

M4 - Recombination Layer: The Two-Site Coordination Challenge

Quantifies how Bxb1 recognizes and recombines two attB sites, and why longer fragments reduce efficiency.

M5 - Experimental Assistant: From Mechanism to Practical Tool

Transforms parameters from the first four layers into a lookup-table success rate predictor for experiment design.

Three Core Discoveries

Discovery 1: The Real Bottleneck Isn’t Where You Think

Most would assume low Bxb1 protein concentration causes poor efficiency. But the model reveals:

• Intracellular Bxb1 reaches 200–400 nM (μM range)
• In vitro experiments need only 10–50 nM for efficient recombination
• So why is efficiency still low in living cells?

The answer lies in M3: randomness of attB site formation is the true killer.

• Single ssAP integration probability into chromosome < 1%
• Success rate strongly depends on chromosomal position:
  • Near oriC: ~0.6%
  • Near terC: ~0.03% (20-fold difference!)

Position Effect Formula: P(d) ≈ 0.6% × e^−0.0007d, where d = distance from oriC (kb).

Discovery 2: The Mechanistic Cost of Two-to-Two

At the same chromosomal position (oriC-proximal), we compared:

• One-to-one: ~0.57% success rate
• Two-to-two: ~0.05% success rate

Why the 10-fold difference? The model provides the answer:

The fundamental mechanism complexity differs, not just a probability multiplication:

Success_one-to-one = P₀ × [modification factors]
≈ 0.001066 × φ_tot × f_gc × f_locus × distance_factor × …
≈ 0.57%

Success_two-to-two = P₀_dual × [same modification factors]
≈ (P₀ / 15–20) × φ_tot × f_gc × f_locus × distance_factor × …
≈ 0.05%

Why is P₀_dual so much smaller? Two-to-two deletion requires:

1. Two independent ssAP integration events both succeed
2. Spatial coordination between distant chromosomal sites
3. Synchronized recombination at both attB sites
4. Any single failure aborts the entire process

Intuitive understanding: It's like needing to successfully thread two needles at once while they're moving—even if you're skilled at threading one needle (one-to-one), doing two simultaneously in coordination is fundamentally harder, not just twice as hard.

Mechanistic complexity multiplier: Two-to-two is 15-20 times harder than one-to-one

Discovery 3: Skip the Bottleneck, Efficiency Soars to 80%+

When using Ready-to-Integrate Mode B (linear DNA with pre-built attB sites):

Average prediction error < 0.4%

The model not only predicts accurately but also explains why:

• Bypasses the M3 bottleneck (no need for ssAP to form attB)
• M4 recombination becomes the dominant step, but remains efficient in the 1-2 kb range
• Bxb1 concentration reaches effective range within 30 minutes

Model Validation: Data Speaks

We calibrated the model using three length anchor points (1039/1605/1613 bp):

Mean Absolute Error (MAE) = 0.37%
Maximum deviation < 1%
Covers practical range 1.0-1.6 kb

All predictions include 95% confidence intervals, with extrapolation regions (>4 kb, extreme GC%) clearly marked.

How to Read This Model?

Choose different reading paths based on your needs:

I want to use it quickly

→ Jump directly to M5 Experimental Design Assistant
Input your fragment parameters and get success rate predictions from the lookup table.

I want to understand why efficiency is low

→ Focus on M3 vs M4 Comparative Analysis
See where the bottleneck truly lies.

I want complete mechanistic understanding

→ Read sequentially M1 → M2 → M3 → M4 → M5
Trace the complete process from induction to integration.

I want to improve or adapt this model

→ Consult parameter tables and calibration guidelines for each layer.
All assumptions and limitations are explicitly stated.

What Does This Model Tell Us?

1. Low efficiency isn't due to insufficient Bxb1, but because attB formation is too difficult.
2. Chromosomal position effects are physical limitations that traditional mode cannot avoid.
3. Two-to-two deletion is inherently harder than one-to-one insertion (mechanistic complexity).
4. Mode B system bypasses the bottleneck, maintaining ≥80% efficiency for 1–2 kb fragments.
5. Simple conditions suffice: 300–400 nM Bxb1, 30 minutes reaches the concentration plateau.

Extended Implications of the Model

Can this modeling framework be transferred to other systems?

Yes! This five-layer architecture applies to any gene editing system requiring:

• Induced expression (M1–M2)
• Chromosomal site marking (M3)
• Sequence-specific recombination (M4)

Examples include:

• Other integrases: Cre-loxP, FLP-FRT, ΦC31
• CRISPR systems: Use M3 to model sgRNA delivery, M4 to model Cas9 cleavage efficiency
• Transposon systems: Tn5, Tn7, etc.

Insights for Modern Gene Editing

The model reveals that "delivery bottleneck > reaction kinetics" — a principle explaining why:

• In vitro assembled RNPs (CRISPR) outperform plasmid expression
• Electroporated DNA is more stable than chemical transformation
• Viral vectors dominate in eukaryotic cells

The underlying logic is the same: bypass random intracellular synthesis/marking processes by directly delivering functional molecules.

Ready to Dive Deeper?

This is not just a mathematical model, but a complete roadmap from gene expression to genome rewriting.

Choose your path and start exploring

Have questions about modeling? Want to see raw data? Click the sections to begin your journey!

Model1(M1): Transcriptional Dynamics of Bxb1 Expression

Model Scope

The M1 model describes the mRNA dynamics of ssAP (single-strand annealing protein) and Bxb1 integrase under various induction regimes.

• The transcriptional outputs m_g(t) provide direct input for M2 protein accumulation.
• Together, M1 and M2 determine whether proteins reach functional steady states required for M3 (attB formation) and M4 (recombination).

A key feature of M1 is its ability to resolve induction time windows, including the effect of cooling periods where transcription is reduced but not fully terminated.

Mathematical Formulation

mRNA dynamics are governed by:

$$ \frac{dm_g}{dt} = \beta_g \cdot I_g(t) - d_{m,g} m_g, \quad g \in \{ \text{Bxb1, ssAP} \} $$

where:

• m_g(t): mRNA concentration (μM)
• β_g: transcription rate (μM/min)
• d_m,g=ln⁡(2)/t_1/2: degradation constant (1/min), defined by half-life
• I_g(t): induction function, varying with experimental design

Model Parameters

For ssAP, we explicitly incorporated a cooling phase to reflect reduced transcriptional activity:

$ I(t) = \begin{cases} 1, & 0 \le t < 45 \quad \text{(full induction)} \\ \gamma, & 45 \le t < 55 \quad \text{(cooling, residual activity)} \\ 0, & t \ge 55 \quad \text{(fully off)} \end{cases} $

with γ = 0.1, representing 10% residual transcriptional efficiency during cooling on ice [4].

Analytical solutions follow canonical forms:

• Continuous induction (I = 1):
  $$ m(t) = \frac{\beta}{d_m} \left(1 - e^{-d_m t}\right) $$
• Switch-off (I = 0):
  $$ m(t) = m(t_{\text{off}}) e^{-d_m (t - t_{\text{off}})} $$
• Steady-state level:
  $$ m_{\text{steady}} = \frac{\beta}{d_m} $$

These are standard analytical solutions [5] obtained by solving the first-order ODE for mRNA dynamics.

Simulation Results

1. ssAP (45 min induction + 10 min cooling)
   • mRNA rose rapidly to a peak of ~1.43 μM at 45 min.
   • During the cooling phase, transcription was not fully halted but reduced to 10% activity, slowing decay.
   • By 120 min, mRNA levels approached baseline, illustrating a transient profile.

   

2. Bxb1 (inducible, 240 min with L-arabinose)
   • mRNA reached 95% of steady state within ~30 min, then maintained a plateau for 4 h.
   • This demonstrates long-term sustained supply consistent with Bxb1’s catalytic role.

   

3. Bxb1 (constitutive expression, promoter J23105)
   • Steady-state mRNA concentration stabilized at ~0.109 μM, significantly lower than inducible conditions.
   • Nevertheless, it provides a low but persistent baseline.

   

Conclusions and Experimental Guidance

• Conclusions:
  • M1 captures distinct temporal regimes: transient (ssAP) vs. sustained (Bxb1).
  • Protein bottlenecks arise not from insufficient magnitude but from induction cutoffs before steady-state.
  • Cooling phases must be considered as low-efficiency transcriptional windows, not as complete shutdowns.
• Guidance:

  
  Extend ssAP induction to 45 min to align with its mRNA peak, accounting for residual transcription during cooling.

Model2(M2): Translation Dynamics and Protein Accumulation

Model Scope

Building upon the transcriptional dynamics characterized in M1, the M2 model quantifies the accumulation of ssAP and Bxb1 proteins.

• The mRNA profiles from M1 serve as direct inputs for translation.
• Protein concentrations are the decisive factors for downstream steps: ssAP enabling attB formation (M3) and Bxb1 catalyzing attP × attB recombination (M4).

Thus, M2 serves as the critical bridge between gene expression and functional integration efficiency. Its role is to evaluate whether extended induction windows are sufficient to ensure proteins act within their effective functional windows.

Mathematical Formulation with Translation Delay

We start from the canonical mRNA–protein dynamics widely used in gene expression modeling.

(1) mRNA dynamics

$$ \frac{dm_g}{dt} = \beta_g I_g(t) - d_{m,g} m_g(t), \quad g \in \{\mathrm{Bxb1}, \mathrm{ssAP}\}. \quad [1][2] $$

This is a first–order linear ODE describing production at rate $\beta_g$ (controlled by induction function $I_g(t)$) and degradation at rate $d_{m,g}$. The analytical solution is obtained using the integrating factor method:

$$ m_g(t) = m_g(0)e^{-d_{m,g}t} + \int_0^t \beta_g I_g(s)e^{-d_{m,g}(t-s)}\, ds. \quad [1][2] $$

Special cases include:

• Continuous induction, $I_g=1$, $m_g(0)=0$:
  $$ m_g(t) = \frac{\beta_g}{d_{m,g}}(1 - e^{-d_{m,g}t}). \quad [1] $$
• Switch-off at $t_{\mathrm{off}}$:
  $$ m_g(t) = m_g(t_{\mathrm{off}}) e^{-d_{m,g}(t - t_{\mathrm{off}})}. \quad [1] $$
• Steady state: $$ m_{g,\mathrm{ss}} = \beta_g / d_{m,g}.\quad [1] $$

Moreover, the degradation constant is related to the mRNA half-life by $d_{m,g} = \ln 2 / t_{1/2}$, which in E. coli has been globally quantified to be on the order of a few minutes [5].

(2) Protein dynamics with translation delay

We incorporate an average translation/maturation delay $\tau$, such that fully functional proteins appear only after time $\tau$:

$$ \frac{dP_g}{dt} = k_{\mathrm{trans},g}\, m_g(t - \tau) - d^{\mathrm{eff}}_{p,g}\, P_g, \quad [3][4] $$

where the effective decay constant combines protein degradation and growth dilution:

$$ d^{\mathrm{eff}}_{p,g} \equiv k^{P}_{\mathrm{deg},g} + \mu. \quad [1][2] $$

The convention is $m_g(t - \tau) = 0$ for $t < \tau$.

(3) Analytical solutions (continuous induction example)

Assume $I_g(t) = 1$ with initial conditions $m_g(0) = P_g(0) = 0$.

• mRNA solution
  $$ m_g(t) = \frac{\beta_g}{d_{m,g}}\left(1 - e^{-d_{m,g}t}\right). \quad [1][2] $$
• Protein solution with delay
  • For $t < \tau$: $$ P_g(t) = 0. \quad [3][4] $$
  • For $t \ge \tau$: $$ P_g(t) = \frac{k_{\mathrm{trans},g}\,\beta_g}{d_{m,g}\,d^{\mathrm{eff}}_{p,g}} \left(1 - e^{-d^{\mathrm{eff}}_{p,g}(t-\tau)}\right) - \frac{k_{\mathrm{trans},g}\,\beta_g}{d_{m,g}\left(d^{\mathrm{eff}}_{p,g}-d_{m,g}\right)} \left(e^{-d_{m,g}(t-\tau)} - e^{-d^{\mathrm{eff}}_{p,g}(t-\tau)}\right). \quad [3][4] $$

This expression shows that the protein dynamics are identical to the no-delay case, but shifted by $(t - \tau)$.

(4) Generalization

For switch-off or pulse induction with residual activity \(\gamma\), one substitutes the time-dependent \(I_g(t)\) into the mRNA solution and propagates it through the protein convolution. The closed form remains valid provided all terms are shifted by \((t-\tau)\) and \(d_p\) is replaced by d_p,g^eff.

Protein removal is modeled through an effective decay constant \(d_{p,\mathrm{eff}}=d_p+\mu\), where \(\mu\) represents dilution due to cell growth and \(d_p\) denotes intrinsic protein degradation. In M2 we assume a constant \(\mu\), corresponding to balanced exponential growth, to keep the focus on transcription–translation–recombination dynamics. A burden-aware formulation \(\mu([\textit{Bxb1}])\) is introduced later as an optional extension.

Simulation Results

1. ssAP (45 min induction + 10 min cooling)

• Protein concentration peaked at ~9.5 μM around 48 min.
• Levels declined to ~6 μM by 120 min due to degradation.
• The transient dynamic matches its role as a short-lived auxiliary factor.

2. Bxb1 (inducible, 240 min with L-arabinose)

• Protein accumulated continuously, reaching ~42 μM at 240 min.
• Due to its long half-life (>16 h), Bxb1 maintained a near-plateau even after induction ceased.
• Represents long-term stable supply, ideal for recombination catalysis.

3. Bxb1 (constitutive expression, promoter J23105, 240 min)

• Accumulated more slowly, reaching ~0.3 μM at 240 min.
• Despite lower levels, still sufficient to support recombination, serving as a baseline self-sufficient regime.

Biological Interpretation

Comparison with Literature

In vitro assays typically employ serine integrases at 200–400 nM, occasionally up to 1000 nM.
Our simulations predicted μM-range intracellular concentrations, one to two orders of magnitude higher.
This confirms that absolute protein abundance is not the limiting factor in vivo.

Temporal Considerations

The critical determinant is not whether proteins reach steady-state, but whether they are present at the right time window.

• ssAP must first bind to the oligo to generate attB sites.
• Bxb1 then needs sufficient time to encounter both the newly formed attB and the pre-installed attP to catalyze recombination.
  Because this step requires two-site coordination and may be limited by conformational accessibility, short induction leaves Bxb1 with too little opportunity to capture both sites.
  Extending induction therefore increases the encounter probability within the effective reaction window, rather than merely raising overall abundance.

Role Differentiation

• ssAP: rapid accumulation and decay → transient support for attB formation.
• Bxb1: gradual but persistent accumulation → sustained catalytic supply, provided the induction window is long enough to cover the recombination step.

Dry–Wet Interaction and Design Decisions

Dry → Wet (Prediction to Guidance):

• Initial inductions (30 min ssAP, 2 h Bxb1) were predicted to leave proteins insufficiently aligned with their effective functional windows, restricting recombination efficiency.
• Based on M1–M2 integration, we recommended extending induction windows: ssAP to 45 min and Bxb1 to 4 h, ensuring transient ssAP action is captured and Bxb1 has enough time for dual-site encounters.

Wet → Dry (Feedback to Refinement):

• Experimental teams extended induction accordingly (ssAP: 30 → 45 min; Bxb1: 2 h → 4 h).
• Feedback from colony counts will recalibrate degradation constants and refine M2 protein predictions.

This cycle demonstrates how modeling can directly shape experimental protocols by moving beyond concentration magnitude and instead emphasizing the importance of temporal coordination and encounter probability.
Without this temporal adjustment, downstream recombination (M4) would systematically underestimate efficiency, despite protein abundance being non-limiting.

Conclusions and Experimental Guidance

Conclusions:

• Intracellular protein levels are already far above in vitro requirements, ruling out insufficiency as a bottleneck.
• The decisive factor is whether induction windows allow ssAP to generate attB and give Bxb1 adequate time to engage both attB and attP.

Experimental Guidance:

• Maintain ssAP induction for 45 min to align with its transient peak and attB formation.
• Extend Bxb1 induction to 4 h to maximize the probability of dual-site recombination events.
  This adjustment ensures effective coordination of M3 and M4, improving recombination efficiency without requiring higher expression strength.

Graphics for Wiki/Publication

• Figure 1: Overlaid protein dynamics curves (ssAP transient, Bxb1 inducible plateau, Bxb1 constitutive slow rise).

• Figure 2: Flow diagram of Dry–Wet interaction: Model prediction → Experimental adjustment (longer induction) → Data feedback → Model refinement.

Model3(M3) Model: Single-Cell Stochastic Process of attB Formation

Model Scope

The M3 model adopts a single-cell, site-discrete probabilistic framework to describe how synthetic oligonucleotides (oligos) are integrated into the E. coli MG1655 chromosome via RecT (ssAP), forming attB sites. Compared to M1 (mRNA dynamics) and M2 (protein dynamics), addressing not only whether a cell forms at least one attB site, but also how genomic position relative to oriC affects integration probability.

Motivation and Design Philosophy

Traditional concentration-based kinetic models fail to capture the polarization phenomenon often observed in recombineering: in the same batch, a minority of cells succeed in integration while most fail completely. This arises because each cell contains only 1–4 genome copies, making the relevant question not the average concentration of intermediates but the probability that a given cell achieves ≥1 attB integration event. This motivated a stochastic, site-limited modeling strategy.

Positioning and Function

• Upstream dependence: Takes ssAP dynamics from M2, converting concentrations to per-cell molecule counts.
• Core feature: Models attB appearance as a non-homogeneous Poisson process, naturally capturing population heterogeneity.
• Outputs: Time landmarks (t₁₀, t₅₀, t₉₀), distributions of attB copy numbers per cell, and optimal Bxb1 induction windows.
• Bridging role: Links the protein concentration layer (M2) with recombination efficiency (M4), defining both the fraction of competent cells and their timing.
• Locus dependence: Integration efficiency declines with increasing distance from oriC, consistent with replication-timing and gene-dosage effects.
• Replication-fork gating: Integration is efficient only when the replication fork passes through the target locus; otherwise probability is near zero.

Mathematical Formulation

(1) Concentration-to-molecule conversion
The ssAP concentration output from M2 is converted into per-cell molecule numbers:

$$ N_{ssAP}(t) = ssAP \times 602 $$

where 602 is the approximate conversion factor for E. coli (1 µM ≈ 602 molecules per 1 fL cell volume) [1]. This links the protein-level dynamics in M2 to discrete single-cell stochastic modeling in M3.

(2) Single-cell reaction rate function (Cycle 4 baseline)
For cell i, with genome copy number (G_i {1,2,4}) (replication-state dependent [2]), the instantaneous rate of attB formation is modeled by a Michaelis–Menten form:

$$ v_i(t) = G_i \cdot \frac{k_{cat} , N_O^{(i)}(t)}{K_m + N_O^{(i)}(t)} $$

where:

• N_O⁽ⁱ⁾(t): available oligos in cell i
• k_cat: catalytic turnover constant (min⁻¹) [3]
• K_m: effective Michaelis constant (molecules) [3]
• G_i: genome copy number (sampled as 1, 2, or 4 [2])

Oligos decay exponentially due to nucleases [4,5]:

$$ N_O^{(i)}(t) = N_O^{(i)}(0) e^{-k_{loss}t} $$

(3) Stochastic integration as a non-homogeneous Poisson process
The cumulative intensity is:

$$ \Lambda_i(t) = \int_0^t v_i(\tau) d\tau $$

and the probability that cell i has formed at least one attB site by time t is:

$$ P_i(t) = 1 - e^{-\Lambda_i(t)}. $$

Monte Carlo sampling (10³–10⁴ cells) over initial oligo distributions (log-normal, GSD ≈ 2–3 [6]), ssAP multiplicative noise (GSD ≈ 1.3–1.5 [7]), and genome copy number distributions yields the population-level fraction of cells with ≥1 attB site.

(4) Cycle 5 extension: locus dependence and replication-fork gating
Cycle 4 assumed all chromosomal loci were equally accessible. In Cycle 5, we introduce explicit locus and time-dependence:

$$ v_i(x,t) = G_i(x,t) \cdot \frac{k_{cat} N_O^{(i)}(t)}{K_m + N_O^{(i)}(t)} \cdot \chi_{\text{rep}}(x,t) $$

where:

• x: genomic locus (distance from oriC in kb)
• G_i(x,t): locus-specific copy number, higher near oriC and lower near terC [2,8]
• χ_rep(x,t): replication-window indicator (1 when the replication fork passes locus x, 0 otherwise [9])

Thus, attB formation can only occur during fork passage, and success probability decreases with locus distance from oriC.

(5) Empirical oriC–ter decay function
Fitting experimental and predicted efficiencies (hisA, metA, galK, leuD, lacZ), the success probability declines approximately exponentially with oriC distance:

$$ P(d) \approx P_0 \cdot e^{-k d} $$

where:

• d: locus distance from oriC (kb)
• P₀ ≈ 0.6%: maximal success rate at ori-proximal loci
• k ≈ 7x10^-4, kb^-1: decay constant

This fit indicates that each additional ~1000 kb from oriC reduces integration probability by ~50–70%, consistent with oriC-proximal loci (hisA, 0.57%) being an order of magnitude more efficient than ter-proximal loci (leuD, 0.05%).

(6) Numerical implementation
Analytical solutions remain tractable only in the simplified Cycle 4 case (constant inputs, no locus dependence). With locus-specific copy number and replication timing, numerical integration combined with Monte Carlo sampling is required. This allows computation of population-level trajectories, polarization histograms, and locus-dependent gradients (Figs. M3-2A–E).

Simulation Results

(1) Baseline molecular framework .
Simulations under baseline parameters (g_m = 300 molecules/cell, k_cat = 0.02 min⁻¹, K_m = 6000 molecules) reproduced the characteristic recombineering behavior. Population success curves exhibited lag–rise–plateau dynamics with a final plateau of ~30–35%, consistent with reported λ-Red/RecT efficiencies. Polarization histograms revealed bimodal outcomes: most cells remained unsuccessful, while a small fraction achieved near-complete integration. Sensitivity analyses confirmed parameter priority as oligo delivery (g_m) > catalytic turnover (k_cat) > K_m, with delivery dominating across regimes.

(2) Locus dependence .
In incorporating locus position and replication-fork gating fundamentally reshaped predictions. Ori-proximal loci (hisA) reached high success probabilities (0.57%), mid-replichore loci (galK) were intermediate (~0.19%), and ter-proximal loci (leuD) dropped to 0.05%. A simple exponential model captured this oriC–ter gradient:

$$ P(d) \approx 0.6\%\,\cdot\, e^{-0.0007d}, $$

indicating each additional ~1000 kb from oriC reduces success probability by ~50–70%. These results demonstrate that positional constraints outweigh steric hindrance: two-to-two recombination at lacZ (0.0195%) did not outperform one-to-one at leuD (0.0386%).

To explicitly incorporate replication timing, Cycle 5 simulations compared representative loci across the E. coli chromosome. As shown in Fig. M3-2E, ori-proximal genes (hisA) reached significantly higher success fractions than mid- and ter-proximal sites (galK, leuD), consistent with experimental observations. This establishes a genome-wide oriC–ter gradient as the dominant factor shaping integration efficiency.

Simulated integration success trajectories for representative chromosomal loci demonstrate a clear oriC–terC gradient. Ori-proximal loci (hisA) reach high plateaus (0.57%), mid-replichore locus (galK) reaches intermediate levels (~0.19%), and ter-proximal loci (leuD) remain very low (0.05%). The curves highlight that genomic position, mediated by replication timing and locus copy number, is the dominant determinant of success, outweighing steric hindrance effects.

Biological Interpretation

• Polarization arises from stochastic molecular heterogeneity (oligo uptake, ssAP variability, nuclease degradation), not population averages.
• S-shaped dynamics emerge naturally from the Poisson framework, reflecting the transient window where oligo availability and RecT activity overlap.
• Locus dependence dominates: ori-proximal loci are ~10–20× more favorable than ter-proximal sites due to earlier replication and sustained higher copy number.
• Steric hindrance is secondary: detectable only at high-efficiency loci; at low-efficiency loci, positional effects override.
• Design implication: integration efficiency is best predicted by combining molecular noise with replication timing.

Dry–Wet Interaction and Design Decisions

Dry → Wet (predictions):

• Optimal induction window aligns with replication fork passage; inducing too early (<15 min) is premature, while induction after the fork has passed results in negligible success at distal loci.
• oriC-proximal loci should be prioritized whenever possible. If distal loci must be targeted, extended induction or oligo stabilization is required.
• Parameter priority remains: (1) oligo delivery, (2) RecT/Bxb1 expression and activity, (3) catalytic fine-tuning.

Wet → Dry (feedback):

• Experimental colony assays confirmed the oriC–ter gradient, motivating explicit fork-gating in the model.
• Time-course qPCR and single-cell FACS measurements of integration efficiency across loci are needed to refine locus-specific parameters.

Conclusions and Experimental Guidance

• attB formation is governed by single-cell stochasticity shaped by genomic context.
• Molecular heterogeneity (oligo delivery, nuclease degradation, ssAP variability) explains polarization and S-shaped dynamics at the population level, establishing the baseline variability upon which locus-dependent replication timing imposes an additional positional gradient.
• Genomic position is the dominant determinant: ori-proximal loci achieve order-of-magnitude higher efficiencies than ter-proximal loci.
• Design rules:
  

1. Select oriC-proximal loci for maximum efficiency.
2. For distal loci, extend induction and stabilize ssDNA to overlap with replication-fork passage.
3. Induce Bxb1 when 30–50% of cells are predicted to carry attB sites (~20–30 min post-electroporation at proximal loci).

• This combined framework unifies molecular noise and chromosomal context, providing actionable guidance for rational locus selection in genome engineering.

Model4(M4) — Dual-End Tetrameric Recombination Core Model

1. Model Scope and Hierarchical Design

Overall Layering and Data Flow

The M4 model operates within a hierarchical framework where upstream models (M1–M3) provide boundary conditions: M1–M2 establish the temporal profile of Bxb1 concentration following induction, capturing the rise, plateau, and potential decay phases [1][2]; M3 quantifies chromosomal attB site availability and the position-dependent efficiency gradient along the oriC–terC axis [3][4]; M4 constitutes the mechanistic core, converting these upstream conditions into dynamic trajectories of success versus failure pathways through explicit treatment of binding, synapsis, and cleavage steps [5][6].

This modular architecture enables independent validation and parametric refinement of each layer,maintaining model interpretability and experimental tractability [7][8].The design supports two-level population mapping: from conditional success rates within the attB-bearing sub-population to whole-population outcomes.Specifically, we define:P_TX = n_attB / n_total (transformation efficiency, e.g., 1700 / 100,000 = 0.017)Example calculation: 1.147 % × PTX = 0.017 = 0.0195 % overall, confirming consistency between conditional and population-level rates.We propagate binomial SE for P_TX in confidence intervals of whole-population rates;central value 0.017 with propagated uncertainty (see Methods).We calibrate parameters using the simpler one-to-one insertion system before extending to two-to-two replacement,avoiding identifiability problems inherent in direct fitting of complex mechanisms [7][8].In our GenOMe system, one-to-one means that one pair of attB and attP sites recombines once through Bxb1 integrase to form one integration event.

Two-to-two, on the other hand, means two pairs of attB and attP sites are available on the genome and the DNA fragment, allowing two recombination events to happen at the same time — one at each site.

One-to-One System: Baseline Calibration

The one-to-one insertion mechanism proceeds through sequential steps: high-affinity recognition (k_on / k_off) → single-end homologous synapsis (k_f1 / k_synr) → catalytic strand exchange (k_cat) → product release with passive attL/attR binding [5][9]. This simpler topology establishes the quantitative relationship between sequence specificity (encoded in site-dependent k_on / k_off ratios) and overall conversion efficiency, providing a validated parameter window for the two-to-two model [10].

Critical correction on product binding: attL and attR sites do not catalytically inhibit the forward reaction. Instead, they function as non-reactive binding sinks that transiently sequester Bxb1 through reversible association, reducing the free enzyme fraction available for productive attP–attB recognition. This passive sequestration becomes significant only at high product accumulation and does not drive reverse conversion in the absence of recombination directionality factor (RDF) [11][12].

$$ \begin{aligned} \text{attP} + \text{Bxb1} &\xrightleftharpoons[k_{-1}]{k_1} C_P \\ \text{attB} + \text{Bxb1} &\xrightleftharpoons[k_{-2}]{k_2} C_B \\ C_P + C_B &\xrightleftharpoons[k_{-s}]{k_s} \text{synapsis} \xrightarrow{k_{cat}} C_L + C_B \\ C_L &\xrightleftharpoons[k_{-3}]{k_3} \text{attL} + \text{Bxb1} \\ C_R &\xrightleftharpoons[k_{-4}]{k_4} \text{attR} + \text{Bxb1} \end{aligned} $$

Two-to-Two System: Non-Length-Increasing Replacement

The two-to-two replacement mechanism embodies “dual-end synchrony”: two spatially separated attP/attB pairs first form independent single-end synaptic complexes (S_left and S_right) with Bxb1 dimers; these must then encounter each other through chromosomal looping to achieve bimolecular ring-closure (rate constant k_r2), assembling a tetrameric synaptic complex (D); finally, concerted cleavage and strand exchange complete the replacement [6][13].

Unlike one-to-one insertion—which opens a chromosomal gap and introduces new DNA—two-to-two performs direct segment replacement between pre-installed attB sites without altering total genome length. Critically, the newly inserted segment can carry orthogonal att sites at its termini, enabling unlimited replacement cycles without cumulative length increase. This architecture enables unlimited “replace-without-growth” cycles and shifts selection from additive stacking to iterative updating.

$$ \begin{aligned} \frac{d[\text{attP}]}{dt} &= -k_{\text{on}}[\text{attP}][\text{Bxb1}] + k_{\text{off}}[C_p] \\ \frac{d[\text{attB}]}{dt} &= -k_{\text{on}}[\text{attB}][\text{Bxb1}] + k_{\text{off}}[C_b] \\ \frac{d[\text{Bxb1}]}{dt} &= -k_{\text{on}}[\text{attP}][\text{Bxb1}] - k_{\text{on}}[\text{attB}][\text{Bxb1}] + k_{\text{off}}[C_p] + k_{\text{off}}[C_b] + 2k_{\text{cat}}[C_{\text{syn}}] \\ \frac{d[C_p]}{dt} &= k_{\text{on}}[\text{attP}][\text{Bxb1}] - k_{\text{off}}[C_p] - k_{\text{syn}}[C_p][C_b] + k_{-\text{syn}}[C_{\text{syn}}] \\ \frac{d[C_b]}{dt} &= k_{\text{on}}[\text{attB}][\text{Bxb1}] - k_{\text{off}}[C_b] - k_{\text{syn}}[C_p][C_b] + k_{-\text{syn}}[C_{\text{syn}}] \\ \frac{d[C_{\text{syn}}]}{dt} &= k_{\text{syn}}[C_p][C_b] - k_{-\text{syn}}[C_{\text{syn}}] - k_{\text{cat}}[C_{\text{syn}}] \\ \frac{d[\text{attL}]}{dt} &= k_{\text{cat}}[C_{\text{syn}}] - k_{\text{on}}[\text{attL}][\text{Bxb1}] + k_{\text{off}}[C_l] - k_{\text{on}}[\text{attL}][\text{Bxb1}] + k_{\text{off}}[C_r] \\ \frac{d[C_l]}{dt} &= k_{\text{on}}[\text{attL}][\text{Bxb1}] - k_{\text{off}}[C_l] \\ \frac{d[C_r]}{dt} &= k_{\text{on}}[\text{attL}][\text{Bxb1}] - k_{\text{off}}[C_r] \end{aligned} $$

Parameter relationship correction:

Notation (first-use definitions). k_f1 (single-end synapsis), k_synr⁻¹ (its reverse), k_f2 (dual-end ring-closure), k_rel (single-end disassembly), k_cat (cleavage). S_left/right: single-end synaptic intermediates; D: dual-end tetrameric complex. We normalize [S] ∈ [0,1] as the fraction of cells (or occupancy) carrying a single-end synaptic intermediate. f_ssAP (0.05–0.10): fraction of chromosomal attB sites activated by RecT/CspRecT-mediated single-strand annealing. a_site (0.28–0.40): position-dependent reduction in association rate due to nucleoid structure. g_loop (0.30–0.40): reduction in ring-closure rate due to DNA topology and confinement.

2.Locus Geometry: Literature Basis and Experimental Context

Both published data and our one-to-one analyses demonstrate a robust positional effect: integration efficiency decreases monotonically with distance from oriC toward terC [3][14]. Quantitative examples include high-efficiency loci (hisA, metA, galK; all oriC-proximal) versus low-efficiency regions (leuD, lacZ; ter-proximal) [15][16].

Our experimental two-to-two design utilized the lacZ region primarily for blue-white screening convenience, accepting the trade-off of operating in a low-efficiency zone where lacZ and leuD nearly coincide.

The 30–60 min / 0.3–1.0 μM window corresponds to the time-over-threshold (ToT) concept originating from M1–M2 models and is localized here as the standard operational (SOP) range for Mode B optimization.

This locus choice explains why the efficiency differential between our two-to-two and one-to-one implementations was smaller than initial predictions: both operated under comparable geometric constraints. Relocating one or both att sites toward oriC would substantially increase overall success probability by enhancing both local DNA accessibility and effective encounter rates [4].

3. Mathematical Formulation

I. One-to-One Baseline Mechanism

Binding Equilibria (Pre-Synaptic Complex Formation)

$$ \frac{d[C_P]}{dt} = k_{on}[\text{attP}][\text{Bxb1}] - k_{off}[C_P] $$

$$ \frac{d[C_B]}{dt} = k_{on}[\text{attB}][\text{Bxb1}] - k_{off}[C_B] $$

These expressions inherit sequence-dependent affinity from heat-map data [10]. We map in vitro specificity ratios to in-cell \(k_{on}/k_{off}\) via proportional scaling while preserving \(K_d\); absolute values are re-fit in Layer 2. Paired adjustments maintain thermodynamically reasonable \(K_d\) values (typically 10–50 nM). High-affinity recognition ensures rapid, reversible pre-synaptic complex formation, establishing the substrate pool for downstream synapsis and cleavage. Sequence-dependent differences manifest primarily at this recognition step, avoiding mis-attribution to catalytic phases.

Single-End Synapsis

$$ \frac{d[S]}{dt} = k_{f1}[C_P][C_B] - k_{synr}[S] - k_{cat}[S] $$

Synaptic complex stability governs the waiting-time distribution before potential dual-end encounters in the two-to-two system. Large \(k_{synr}\) values induce rapid cycling between pre-synaptic and synaptic states, reducing the steady-state synaptic intermediate pool. This step exhibits strong concentration dependence and sensitivity to the Bxb1:DNA stoichiometric ratio.

Product Binding as Passive Sinks

In the absence of RDF, we do not include reverse catalysis. Instead, attL/attR contribute to enzyme sequestration:

$$ \frac{d[C_L]}{dt} = k_{on}^{attL}[\text{attL}][\text{Bxb1}] - k_{off}^{attL}[C_L] $$

$$ \frac{d[C_R]}{dt} = k_{on}^{attR}[\text{attR}][\text{Bxb1}] - k_{off}^{attR}[C_R] $$

These enter the Bxb1 mass balance but leave product concentrations unchanged; we do not subtract binding terms from d[attL]/dt or d[attR]/dt. Products are not consumed by reversible binding; the effect is reduced free enzyme, not catalytic inhibition [11][12].

These terms subtract only from the free [Bxb1] pool and do not enter d[attL]/dt or d[attR]/dt, thereby preserving mass balance without introducing reverse catalysis.

II. Two-to-Two Replacement Mechanism

Dual-End Ring-Closure as the Bottleneck

Two single-end intermediates \( S_{left} \) and \( S_{right} \) must meet to form \( D \).

Defining \( J_{f2} = k_{f2} [S_{left}] [S_{right}] \):

$$ \frac{d[S_{left}]}{dt} = \dots - J_{f2} - k_{rel}[S_{left}] $$ $$ \frac{d[S_{right}]}{dt} = \dots - J_{f2} - k_{rel}[S_{right}] $$ $$ \frac{d[D]}{dt} = + J_{f2} - k_{cat}[D] $$

Under \( [S_{left}] = [S_{right}] = [S] \), dual-end encounter is second-order in [S], yielding:

$$ P_{ring} = \frac{k_{f2}[S]^2}{k_{f2}[S]^2 + k_{rel}} $$

Here [S] is a normalized single-end pool (fraction of cells or dimensionless occupancy), so \( k_{f2}[S]^2 \) and \( k_{rel} \) are effective first-order terms and \( P_{ring} \) is dimensionless. This square dependence explains the sharp, sigmoidal gain when ends are closer; \( k_{rel} \) captures the decay of single-end intermediates outside the synchrony window.

Reverse directionality (RDF) policy

gp47 (RDF) was not expressed in any condition; models therefore set k_rev = 0 and treat the forward pathway as effectively irreversible after dual-end closure. Any effects of attL/attR are modeled only as passive, reversible enzyme-sequestration terms in the Bxb1 mass balance (no catalytic inhibition) [11][12]. As RDF is excluded by design, we do not analyze reversibility further in this work.

4. In-Cell Effective Parameter Mapping

We avoid double counting geometry by applying it only to the steps it governs:

$$ I_{\text{eff}} = I_{\text{tot}} \times f_{\text{ssAP}} $$ $$ k_{\text{on,eff}} = a_{\text{site}} \cdot k_{\text{on}} \quad \text{(association only)} $$ $$ k_{f2,\text{eff}} = g_{\text{loop}} \cdot k_{f2} \quad \text{(ring-closure only)} $$

We map chromosomal position and locus accessibility to \( a_{\text{site}} \) (association only), and polymer looping/topology (length, orientation, torsion) to \( g_{\text{loop}} \) (ring-closure only); no geometry factors are applied elsewhere [3][4][17][18].

Ranges reflect our fit priors anchored to macrodomain/topology literature; see Table S2 for locus-wise posteriors. Reported I_eff (1–2 μM) reflects M2/M3-inferred f_ssAP under our induction regime and should be re-calibrated if growth rate, copy number, or crowding differ.

5. Length-Dependent Kinetics

Length penalties apply exclusively to ring-closure and disassembly; recognition parameters are length-invariant:

$$ k_{f2}(L) = k_{f2}^{\circ} \cdot \exp[-\alpha (L - L_0)] $$ $$ k_{rel}(L) = k_{rel}^{\circ} \cdot \exp[+\beta (L - L_0)] $$

where \( L_0 = 1039 \, \text{bp} \). Fitting the 1039/1605/1613-bp triad yields (α, β) with 95% CIs (see Table S3). Longer segments reduce end-to-end encounters (\( k_{f2}\downarrow \)) and increase strain-induced disassembly (\( k_{rel}\uparrow \)) [17][18][19][20][21]; shorter segments elevate \( P_{ring} \) and accelerate completion.

6. Parameter Identifiability and Calibration Hierarchy

This hierarchy prevents over-parameterization [7][8]. Minimal observation sets per layer:

• L1 — single-site heat-map ratios
• L2 — one-to-one time-course (≥2 time points pre/post 60 min)
• L3 — two-to-two conditional success at 240 min plus ≥2 length points
• L4 — locus-wise (oriC→ter) comparisons

7. Core Parameter Summary

A compact set of calibrated parameters is summarized in Table 1. Full parameter tables, initial conditions, and locus-wise geometry estimates are provided in Tables S1–S3 .

Table 1. Core Kinetic Parameters

8. Simulation Results

A. One-to-One Baseline

Pre-synaptic complex formation reaches quasi-equilibrium rapidly (\( \tau \approx 5 \, \text{min} \)); product binding becomes relevant only in late phases (>60 min) when attL/attR accumulation exceeds ~30% conversion. This baseline establishes that if two-to-two underperforms, the bottleneck lies in \( k_{f2} \) versus \( k_{rel} \) competition rather than deficient recognition.

B. M4₁₆ Results (Two-to-Two with Dynamic Bxb1)

Bxb1 dynamics: The total translated Bxb1 (I_tot) rises during 0–240 min and approaches the order of 10 μM in our simulations (model outputs informed by M2/M3). The effective reactive fraction (I_eff = I_tot × f_ssAP × f_free) remains around 1–2 μM and serves as the input concentration for downstream kinetic models. This scale is consistent with the 0.3–1.0 μM operational window defined in M5 and used for Mode B optimization.

Dual-end competition: The success pathway flux (J_f2 = k_f2 [S_left][S_right]) directly competes with failure (k_rel [S_i]), determining whether dual-end ring-closure completes within the synchrony window before single-end intermediates commit to failure pathways.

Population mapping: Calibrated to 1.147% conditional success at 240 min within \( n_{attB} = 1700 \) cells; multiplication by \( P_{TX} = 0.017 \) yields 0.0195% whole-population rate.

C. M4₂₃ Results (Steady-State Multi-Round)

Under constant Bxb1 (~4.36 μM), success is governed by \( k_{f2}/k_{rel} \) balance and fragment length. The observed ordering (1039 bp: 87.5%, 1605 bp: 80%, 1613 bp: 79%) emerges purely from length-dependent \( k_{f2}/k_{rel} \) modulation.

Bars represent mean success percentage across \( n = 3 \) biological replicates (error bars: ±1 SD). Fragment lengths: 1039 bp (reference), 1605 bp, 1613 bp. Solid line: model prediction using \( k_{f2}(L) = k_{f2}^{\circ} \exp[-\alpha (L - L_0)] \) with \( \alpha = 1.5 \times 10^{-3} \, \text{bp}^{-1} \) fitted to data. Dashed lines: 95% prediction interval. Bootstrap resamples were drawn at the replicate level with a fixed RNG seed for reproducibility; point estimates are medians, and bands denote 95% CIs.

9. Biological Interpretation

Near-Irreversibility Threshold

Without RDF, the forward pathway after dual-end closure is effectively irreversible under our conditions; we did not observe reverse flux, and products enter only via passive enzyme sequestration (no catalytic inhibition). This explains why increasing \( k_{f2} \) (or effective [S]) dramatically shifts overall success. Conversely, “waiting longer” provides diminishing returns once the synchrony window closes: single-end-committed failures become kinetically locked.

Sequence Geometry and Length Effects

DNA looping probability scales approximately as \( P_{loop} \propto L^{-3/2} \) for \( L > \) persistence length (~50 nm ≈ 150 bp) [19][20][21]. For our 1.0–1.6 kb range, this predicts substantial reduction in encounter frequency, consistent with observed efficiency drops. Shorter segments elevate \( P_{ring} \) faster than they increase torsional strain because bending energy scales as \( \kappa \cdot \theta^2 / L \), while encounter frequency scales as \( L^{-3/2} \).

Chromosomal Position (oriC Distance)

Ori-proximal loci exhibit higher supercoiling density, greater replication-associated accessibility, and reduced NAP crowding [3][4]. These effects compound into \( a_{site} \) gradients ranging from ~0.40 (hisA, 0.2 Mb from oriC) to ~0.28 (lacZ, 1.5 Mb from oriC), directly modulating \( k_{on,eff} \) (see Table S2).

Orthogonal Design Benefits

Operationally, we require \( K_{d,cross}/K_{d,self} > 100 \) and observe <1% cross-recombination in plasmid challenge assays; see Methods. This validates clean heat-map → \( k_{on} \) translation and enables stable multi-round performance by preventing inter-round cross-talk. We require (i) <1% cross-recombination in plasmid challenge assays and (ii) no off-target junctions by junction-specific PCR across two rounds.

10. Structural and Geometric Intuition

Dual-end recombination assembles two Bxb1 dimers (one per att pair) into a tetrameric synaptic complex. When sites are closely spaced (≈1–1.3 kb) and co-oriented, DNA bending strain remains small, elevating effective \( k_{f2} \). Beyond a critical spacing, DNA must loop back under torsional constraint, increasing \( k_{rel} \) probability.

This geometric framework explains the “shorter is better” principle: confinement that might hinder one-to-one insertions (steric crowding) actually facilitates two-to-two replacements by constraining the search volume for dual-end encounters. In molecular visualization, the tetramer functions as a rigid hinge—once both ends dock, cleavage proceeds with near-irreversible commitment.

11. Wet-Lab Verification Pathways

To experimentally validate M4 predictions:

1. Position gradient assay: Fix attP-plasmid and relocate genomic attB pairs along oriC → terC axis (metA → hisA → lacZ → galK → leuD). Success probability should decay monotonically with increasing terC proximity.
2. Length-dependence experiment: At constant locus (lacZ), vary replacement lengths (1.0, 1.3, 1.6 kb). Success should follow \( P_{succ} \propto \exp[-\alpha \cdot L] \) with α within the predicted confidence interval.
3. Combined optimization: Shift one attB terminus toward oriC while maintaining total length. The model predicts substantial gain in \( P_{ring} \) for moderate relocation.
4. Reversibility control: Introduce RDF (gp47) expression. The reverse branch should appear only after gp47 induction, validating the forward-only assumption under standard conditions.

12. Model Limitations and External Factors

While M4 captures essential dual-end recombination kinetics, several biological processes remain simplified:

• DNA repair competition: Mismatch repair (MMR, \( \tau \approx 2\text{–}5 \, \text{min} \)) and RecBCD nuclease can prematurely erase single-end intermediates. Inclusion would require coupling to the RecA/RecBCD network.
• Topological constraints: Supercoiling density fluctuations and NAP binding dynamics are averaged into \( a_{site} \) and \( g_{loop} \).
• Resource competition: RecT/CspRecT translation competes with ribosomal capacity, implicit through M2/M3 outputs.
• Spatial heterogeneity: Nucleoid macrodomains (Ori, Ter, Left, Right) exhibit distinct physical properties, coarse-grained to single parameters.
• Cell-to-cell heterogeneity: Oligonucleotide uptake and ssAP expression vary substantially across cells due to stochastic gene expression [22][23].

Quantitative predictions are most trustworthy for 1–2 kb segments in mid-log E. coli (\( \text{OD}_{600} = 0.4\text{–}0.6 \)) where repair activity and topological state are stable.

13. Conclusions and Actionable Guidance

Key Conclusions

1. Central bottleneck: \( k_{f2} \) vs \( k_{rel} \) competition (\( P_{ring} = \frac{k_{f2}[S]^2}{k_{f2}[S]^2 + k_{rel}} \)) determines success/failure branching; maximizing \( [S]_{eff} \) elevates \( P_{ring} \).
2. Non-length-increasing advantage: Two-to-two enables unlimited cycles without genome growth, transforming selection from additive to modular.
3. Design levers: Position (\( a_{site} \) variation) and length (α, β parameters) are most sensitive parameters.
4. Orthogonality requirement: Specificity ratio > 100:1 minimizes off-target recombination to < 1% per round in our assays.

Actionable Experimental Guidance

1. Locus selection: Achieve ≥1 oriC-proximal terminus; if constrained, compensate via length reduction (~200 bp yields substantial \( P_{ring} \) gain).
2. Length planning: Prefer ≤1.1 kb; for ≥1.6 kb, use sequential rounds (efficiency_{per_round^0.5} > efficiency_{single_long}).
3. Calibration hierarchy: Heat-map → \( k_{on}/k_{off} \) → one-to-one (\( k_{f1}, k_{cat} \)) → two-to-two (\( k_{f2}, k_{rel} \)) → population (\( P_{TX} \)).
4. Reversibility control: Omit RDF for forward-only; introduce gp47 only when designed reversibility required.
5. Multi-round QC: Execute rounds sequentially with inter-round PCR validation to confirm orthogonality and expected amplicon sizes.

Table S1. Complete Parameter Set

Length-dependence update: The early fit yielded α = 1.5 × 10⁻³ and β = 1.0 × 10⁻³ bp⁻¹. Subsequent tri-point regression across 1039/1605/1613 bp constructs refined these to α_mech = β_mech = 8.5 × 10⁻⁴ bp⁻¹ (SLOPE = 1.70 × 10⁻³ bp⁻¹), ensuring consistency with the logit-linear fit reported in M5 (Fig. 1).

Table S2. Position-Dependent Parameters

Table S3. Length-Dependent Fitting Results

Fit statistics:
α̂ = 1.5 × 10⁻³ bp⁻¹ (95% CI: 1.2–1.9 × 10⁻³)
β̂ = 1.0 × 10⁻³ bp⁻¹ (95% CI: 0.8–1.3 × 10⁻³)
RMSE = 0.53%; R² = 0.991; Cross-validation residual < ±1.2%.

M5 — Design-Factor Predictive Model for Genome Integration

Overview

After establishing mechanistic insight across M1–M4, we developed M5 as a compact, interpretable, and experimentally actionable predictive model.
While early layers involved >40 kinetic and structural parameters, sensitivity analysis showed that recombination outcomes are consistently dominated by only a small subset of biological determinants.

M5 therefore serves two roles:

1. A unifying abstraction of the mechanistic dynamics discovered in M1–M4.
2. A practical predictive engine enabling real-time design decisions inside GenOMe Navigator, supporting both
   • Mode A: initial attB installation
   • Mode B: two-to-two cassette integration.

This factor-based reduction converts complex ODE/stochastic simulations into a form that experimental users can apply directly.

Architecture Across Scales

M5 integrates three biological layers into a single predictive output:

1. Protein Layer (from M2)

Determines whether Bxb1 remains above its functional activation threshold long enough to support recombination (captured by φ_tot).

2. DNA Layer (from M3 & M4)

Encodes genomic-position effects, GC-dependent oligo stability, dual-end synapsis geometry, spacer dependence, and length-penalty behavior.

3. Population Layer

Maps per-cell recombination probability into expected colony counts, accounting for uptake heterogeneity.

Across these scales, M1–M4 consistently highlighted five dominant factors, which M5 isolates into a single multiplicative framework.

Five Design Factors

For any given experimental condition, M5 predicts final integration success as:

$$P_{\mathrm{final}}=\alpha, f_{\mathrm{pos}}(d_{\mathrm{ori}}), f_{\mathrm{dist}}(L,s), f_{\mathrm{orient}}(\theta), f_{\mathrm{GC}}(g), f_{\mathrm{prot}}(\phi_{\mathrm{tot}})$$

Where:

1. Genomic Position (f_pos)

Reflects the oriC–ter replication gradient discovered in M3.
ori-proximal loci exhibit sharply higher attB formation.

2. Distance / Geometry (f_dist)

Encodes spacer-length optima and length-dependent penalties revealed in M4.

3. Orientation (f_orient)

Penalizes unfavorable attB/attP polarity combinations that reduce synapsis probability.

4. GC Stability (f_GC)

Models reduced oligo survival and pairing at extreme GC%.

5. Protein Sufficiency (f_prot)

Derived from φ_tot, reflecting induction time and Bxb1 expression sufficiency.

Interpretation:
  M5 is not an empirical curve fit—it is a mechanistic factorization of the causal dependencies discovered in M1–M4.

Limiting-End Principle for Two-to-Two Recombination

Dual-end recombination requires both attB–attP pairs to synapse successfully.
Mechanistic analysis (M4) shows that:

• DNA looping energy is dominated by the most difficult bend.
• Synapsis is a coupled, not independent, process.
• Therefore the less favorable end overwhelmingly determines the final success.

Thus, Mode B predictions use:

$$P_{\mathrm{2to2}} \approx \min(P_1, P_2)$$

This approximation preserves biological realism while remaining computationally simple.

Mode A vs. Mode B

Mode A — Initial attB Installation

Dominated by:

• genomic position (oriC–ter gradient)
• GC content
• ssAP window
• uptake heterogeneity

Typical success rate: 0.001–0.05%.

Mode B — Production Cassette Integration

Dominated by:

• insert length penalty
• spacer geometry
• orientation
• limiting-end principle

Typical success rate: 50–90%.

M5 automatically switches between these regimes within the Navigator interface.

Model Validation Across Genomic Loci

(Integrated from Text 1 — essential for credibility)

To test whether the factorized M5 model captures real genomic behavior, we compared predictions to experimental measurements across five loci:

• hisA (ori-proximal)
• metA (ori-proximal)
• galK (mid-range)
• leuD (mid-range, single-end)
• lacZ (ter-proximal, two-to-two)

Key agreement points

1. ori-proximal loci (hisA, metA) showed highest efficiency,
   fully consistent with the f_pos factor.
2. lacZ displayed the lowest efficiency,
   explained by its ter-proximal location and the limiting-end penalty from its two-to-two architecture.
3. galK and leuD matched mid-range predictions,
   confirming the GC–distance and geometry terms.

Together, these results validate that the five-factor structure of M5 captures the dominant determinants of recombination efficiency across the chromosome.

Engineering Design Rules Derived from M5

(Integrated from Text 1 — gives actionable value)

These rules form the decision backbone of GenOMe Navigator, enabling users to diagnose limiting factors and optimize constructs before performing experiments.

Lookup-Table Workflow (Software Integration)

M5 supports an intuitive three-step workflow:

1. Select Mode (A or B)
2. Provide design parameters (distance to oriC, GC, insert length, spacer, induction)
3. Navigator computes P_final, identifies bottlenecks, and generates ranked suggestions

This connects the mechanistic model to real experimental planning.

Cross-Strain Recalibration

Because all mechanistic terms are factorized, porting M5 to new strains or integrases requires only:

• 3–5 reference measurements,
• re-scaling α,
• adjusting φ_tot if growth kinetics differ.

Geometry and GC factors remain universal.

Summary

M5 converts the mechanistic insights from M1–M4 into a compact, interpretable, and experimentally actionable model. It preserves biological causality while enabling real-time prediction and design optimization through GenOMe Navigator.

Supporting Trends from M1–M4

Below, we summarize the key quantitative trends emerging from M1–M4 that justified the final five-factor structure of M5. These heatmaps and sweeps are not supplementary decoration—they represent the mechanistic evidence that allowed us to simplify dozens of kinetic parameters into a compact, interpretable predictive model.

Each trend corresponds to one biological layer and directly supports one or more terms in M5.

Trend 1 — Protein Layer (M2): Functional Availability Plateau (φ_tot)

This analysis revealed a broad plateau where φ_tot becomes insensitive to both induction time and enzyme concentration. Specifically:

• φ_tot rapidly increases during the first 20–30 minutes,
• and saturates for 0.3–1.0 μM Bxb1,
• with almost no additional gain beyond 30 minutes.

This explains why:

1. Mode B predictions are weakly dependent on Bxb1 concentration,
2. A standard 30-min induction reliably reaches functional enzyme levels,
3. The protein-layer term in M5 can be represented by a single normalization factor f_prot.

Supports: Protein sufficiency factor (f_prot).

Trend 2 — DNA Layer (M3): Replication Timing × GC% Map

Single-cell stochastic simulations showed two strong effects:

1. oriC–ter replication gradient
   • Loci closer to oriC display sharply higher attB-installation probability.
   • Efficiency drops ~10–20× toward the terminus.
2. GC% penalty
   
   • GC extremes (<40% or >65%) reduce oligo survival and pairing probability.
   • A mild optimum appears near ~50% GC.

These two effects were highly consistent across parameter perturbations, providing the mechanistic basis for incorporating genomic position and GC stability as explicit factors in M5.

Supports: f_pos (locus effect), f_GC (GC penalty).

Trend 3 — Geometry (M4): Spacer × Insert Length Synergy

M4 simulations revealed a geometric optimum for dual-end capture:

• A broad maximum occurs at ~60–100 bp spacers,
• Very short (<40 bp) or long (>150 bp) spacers disrupt synapsis,
• Insert length introduces a monotonic penalty, becoming dominant above 2 kb.

This trend provides the structural justification for:

• Treating spacer as a geometry-specific factor,
• Treating insert length as a global multiplicative penalty, and
• Using the limiting-end principle for two-site recombination.

Supports: f_dist (length + spacer geometry).

Trend 4 — Length-Penalty Curve (M4): Dominance of Insert Size

A striking observation across M4 is that:

• Length dictates 90–95% of the variance in Mode B,
• Concentration sweeps converge onto nearly identical curves,
• Experimental validation points fall directly on the predicted curve.

This validated the use of a logit-linear length penalty inside M5 and supported the engineering claim that insert length is the primary determinant for Mode B efficiency.

Supports: f_dist (dominant length effect).

Trend 5 — Population Layer: Uptake Mean × CV Map

M3–M4 propagation into M10k-scale population simulations revealed:

• Low mean uptake limits the available attB precursor pool.
• High CV amplifies left-tail events where a large subpopulation receives almost no oligo.
• These effects are mode-dependent: strong for Mode A, marginal for Mode B.

This justified the inclusion of the population-level scaling factor α, which differentiates the two integration modes.

Supports: α (population correction for Mode A/B).

Trend 6 — Steric Kinetic Factors (M4): kf2 × krel Sensitivity

While not directly tunable during experiments, sensitivity analysis showed:

• Length penalties emerge naturally in regions where synapsis formation slows (kf2↓)
• or release becomes rate-limiting (krel↑).

This validated that the geometry/length effects in M5 do not require ad-hoc curve fitting, but arise from fundamental kinetic asymmetry.

Supports: Mechanistic plausibility of f_dist and orientation factor.

Trend 7 — CFU → Expected Positive Colonies

This mapping connects:

• single-cell probability (P_final)
  
• → population outcome (expected positives)
  

This trend allowed M5 to integrate with the software’s AI suggestion system (e.g., advising additional plates if uptake is poor).

Supports: Population interpretation of P_final.

Dry-Lab Results

Through the M1–M5 modeling framework, we established an interpretable, quantitative system that deconstructs Bxb1-mediated genome integration into three interconnected layers — protein, DNA, and population.
Both simulations and experimental data consistently revealed that integration success is not determined by a single kinetic parameter, but by the combined influence of replication timing, DNA geometry, and protein expression dynamics.

1. Locus Effect Dominates (Replication Timing)

Integration efficiency exhibited a clear hierarchy across genomic loci: hisA and metA achieved the highest success rates, galK showed intermediate performance, while leuD and lacZ were markedly lower.
This gradient mirrors the loci’s relative distance from oriC, indicating that replication fork timing is the primary determinant of integration efficiency.
Regions near oriC are duplicated early and remain more accessible to recombineering enzymes, whereas loci near ter replicate later, limiting accessibility.

Notably, lacZ is physically positioned near hisA and galK, not in the terminal region. Its low success rate—comparable to leuD—stems from its two-to-two configuration, which requires recombination across two pairs of attB/attP sites, unlike the single pair used in one-to-one designs.
This contrast reveals that integration efficiency is shaped not only by replication timing but also by geometric coordination within the DNA architecture.

2. Geometry as a Limiting Factor

In a one-to-one configuration, recombination proceeds independently at a single attB/attP pair.
By contrast, the two-to-two configuration demands synchronized recombination at both ends, introducing challenges in molecular alignment and spatial coordination.
Our model predicts that if one end completes recombination prematurely, Bxb1 molecules may become transiently trapped at that site, leading to temporal asymmetry and stalling of the opposite reaction.
This explains the lower initial efficiency observed for lacZ, which reflects geometric instability rather than locus position.

3. Turning Geometry into Advantage — The Strength of Two-to-Two Integration

While the two-to-two configuration introduces additional geometric constraints, it uniquely enables full cassette integration within a single step, a feat impossible for one-to-one architectures.
Our simulations revealed that the inefficiency of dual-end recombination is not intrinsic to the mechanism but results from unoptimized parameters such as inter-site spacing, enzyme concentration, and induction timing.
Once these were optimized in the M5 framework, the predicted success rate of two-to-two events converged with that of one-to-one integrations.

This finding redefines two-to-two recombination as a design-driven advantage rather than a structural limitation.
By quantifying geometric coordination, the model transforms a former bottleneck into a high-performance, scalable solution for production-level genome engineering — proving that dual-end integration can be both complete and efficient when designed systematically.

4. Distinction Between the First and Second Two-to-Two Rounds

The first two-to-two process (installation of dual attB sites) and the second two-to-two process (dual-end cassette integration) represent distinct mechanistic stages.
The first is probability-limited, governed by ssAP expression level and replication-fork timing, resulting in a very low success rate (~0.00195%).
The second occurs only in cells already containing dual attB sites; it is kinetics-limited, dominated by Bxb1 binding and catalytic efficiency.
Under optimized conditions (0.3–1.0 μM Bxb1 and short inter-site spacing), integration proceeds rapidly, reaching success rates exceeding 80%.

This dramatic contrast marks a shift from probabilistic to kinetic limitation, clarifying why integration efficiency can increase by four orders of magnitude between the two stages.

5. DNA Length Defines an Efficiency Window

Insert length strongly affects dual-end recombination success.
Both simulations and experiments show that fragments of 1–2 kb yield consistently high efficiencies, while inserts longer than ~3 kb result in a sharp decline.
This pattern reflects a trade-off between reaction kinetics and molecular proximity: longer fragments reduce the likelihood that both ends align simultaneously, thereby lowering the probability of complete integration. Thus, insert length emerges as a highly tunable design parameter in optimizing two-to-two architectures.

6. Protein Timing Raises the Integration Ceiling

The M1–M2 models predict optimal induction durations of ≈45 minutes for ssAP and ≈240 minutes for Bxb1. When residual transcription during the cooling phase was incorporated, the simulated protein accumulation curves matched experimental observations. These results show that fine-tuning expression timing determines the achievable “integration ceiling”: operating within the predicted windows significantly increases both attB formation and total integration success.

7. Population-Level Heterogeneity Explained by Replication Gating

The M3 model simulates attB formation as a stochastic single-cell event and successfully reproduces the experimentally observed “polarized” outcomes, where most cells fail to integrate while a minority succeed. This heterogeneity arises from replication-fork gating: integration occurs primarily during the brief window when the replication fork passes through the target locus, rendering DNA temporarily accessible to ssAP-mediated recombination. Replication timing, therefore, is the underlying driver of both population heterogeneity and the oriC–ter efficiency gradient.

8. Integration of All Layers — From Model to Design Assistant

Together, these models define three interacting principles:

The M5 model, implemented as our open-source software GenOMe Navigator, integrates these three layers into a predictive design assistant. Users can input locus position, insert length, and induction schedule to estimate integration efficiency and receive quantitative SOP recommendations, transforming experimental planning from trial-and-error into rational, data-driven design.

9. Future Refinements

While the M5 framework already reproduces experimental trends with strong fidelity,several refinements could further elevate its predictive precision:

1. Expanding locus-level replicates to enhance statistical robustness;
2. Incorporating direct replication-timing data instead of inferred positions;
3. Refining geometric loss modeling for complex two-to-two architectures;
4. Automating parameter calibration across experimental batches.

These extensions represent next-generation upgrades rather than current limitations— laying the foundation for M5 to evolve into a fully predictive genome-design platform.

10. Summary

Integration of modeling and experimental data yields the following key insights:
Integration efficiency follows an oriC-to-ter gradient defined by replication timing.
Two-to-two geometry, once optimized, achieves high-efficiency dual-end recombination comparable to one-to-one events.
The first and second two-to-two rounds are governed by distinct kinetic bottlenecks — the former by attB formation probability, the latter by enzymatic and geometric factors.
A 1–2 kb length window represents the optimal trade-off between completeness and efficiency.
Protein induction timing defines the integration ceiling.
The M5 model (GenOMe Navigator) converts these multilayer rules into a reproducible, quantitative, and design-driven framework for genome engineering. In essence, integration success emerges from the interplay of timing, geometry, and kinetics. The locus provides opportunity, the geometry sets boundaries, and protein dynamics determine attainment.
Through calibrated modeling, these elements are unified into a single predictive platform — transforming genome integration from empirical trial-and-error into a rational, engineering discipline.