YOKABIO

1. Complete Parameter Reference
2. Mathematical Properties of the Three-Parameter Decomposition
3. Numerical Implementation
4. Heatmap Generation Technical Details

1.1 Interactive Simulator Parameters

The interactive simulator uses a three-parameter decomposition approach to explore how modifications affect prey growth kinetics. Instead of directly adjusting enzyme and template concentrations, users control three scaling factors that represent elementary processes. All parameters can be adjusted via sliders with live updates to the time series visualization.

Enzyme and template concentrations are fixed at standard values (pol = 3.7 nM, rec = 32.5 nM, G = 150 nM), while kinetic parameters $k_1$ and $b$ are modified according to:

$$k_1' = k_1 \cdot \frac{r_{\mathrm{assoc}} \cdot r_{\mathrm{poly}}}{r_{\mathrm{nick}}}, \quad b' = b \cdot \frac{r_{\mathrm{assoc}}}{r_{\mathrm{nick}}}$$

Physical interpretation: These three parameters correspond to distinct mechanisms by which chemical modifications (e.g., terminal amino acids on template DNA) can alter oscillation dynamics. Each modification can affect multiple elementary processes simultaneously—for example, altering association strength while also changing turnover rate or saturation threshold. The framework enables systematic exploration of how these combined effects determine oscillation characteristics.

1.2 Fixed Kinetic Parameters (PP₁ System)

The following kinetic parameters are fixed to the optimized values obtained from the original paper's Supporting Information (Table S5 for PP₁ system). These values were determined through experimental fitting and are used as the baseline for all simulations.

1.3 Initial Conditions and Simulation Settings

1.4 Dimensionless Parameters

For bifurcation analysis and theoretical understanding, the system can be expressed in dimensionless form. The key dimensionless parameters are:

These dimensionless parameters completely characterize the system's dynamical behavior. The bifurcation structure (extinction, coexistence, oscillations) is determined by the values of $g$, $\beta$, $\lambda$, and $\delta$.

2.1 Definitions

The three-parameter decomposition framework is based on the following definitions. The effective kinetic parameters $k_1$ and $b$ that appear in the prey growth term are related to the underlying biochemical parameters by:

$$k_1=\frac{K_a^{GN}k}{K},\quad b=\frac{K_a^{GN}}{K} \tag{1}$$

where $K_a^{GN}$ is the G:N association constant, $k$ is the polymerase turnover rate, and $K$ is the Michaelis-Menten constant for the nicking process.

Any experimental condition change (terminal modifications, temperature, salt concentration, etc.) can be represented by three scaling factors:

$$r_{\mathrm{assoc}}\equiv \frac{K_a'}{K_a},\quad r_{\mathrm{poly}}\equiv \frac{k'}{k},\quad r_{\mathrm{nick}}\equiv \frac{K'}{K} \tag{2}$$

where prime symbols denote values after the modification.

These scaling factors transform the apparent kinetic parameters according to:

$$k_1' = k_1\cdot\frac{r_{\mathrm{assoc}} \cdot r_{\mathrm{poly}}}{r_{\mathrm{nick}}},\qquad b' = b\cdot\frac{r_{\mathrm{assoc}}}{r_{\mathrm{nick}}} \tag{3}$$

2.2 Representation Theorem

Theorem (Representation Theorem): If a condition change modifies the values of $(K_a^{GN}, k, K)$ while preserving the functional form of prey growth $\varphi_{N\triangleright N}=k \cdot \mathrm{pol}\cdot\frac{K_a^{GN}GN}{K+K_a^{GN}GN}$, its effect can always be completely expressed by three multipliers $(r_{\mathrm{assoc}}, r_{\mathrm{poly}}, r_{\mathrm{nick}})$ in the form of equation (3).

Proof: From definition (2), $K_a'=r_{\mathrm{assoc}}K_a$, $k'=r_{\mathrm{poly}}k$, and $K'=r_{\mathrm{nick}}K$. Substituting these into the definition of equation (1), $k_1=K_a k/K$ and $b=K_a/K$, we obtain:

$$k_1'=\frac{K_a'k'}{K'}=\frac{(r_{\mathrm{assoc}}K_a)(r_{\mathrm{poly}}k)}{r_{\mathrm{nick}}K} = k_1\cdot\frac{r_{\mathrm{assoc}}r_{\mathrm{poly}}}{r_{\mathrm{nick}}}$$ $$b'=\frac{K_a'}{K'}=\frac{r_{\mathrm{assoc}}K_a}{r_{\mathrm{nick}}K} = b\cdot\frac{r_{\mathrm{assoc}}}{r_{\mathrm{nick}}}$$

Thus, equation (3) holds and $\varphi_{N\triangleright N}$ is maintained in the same functional form. □

2.3 Constructibility Corollary

Corollary (Constructibility): Suppose a change in apparent parameters $(k_1, b) \to (k_1', b')$ is observed. If it comes from the same function family of equation (1) (low-occupancy approximation and MM-like saturation assumptions are preserved), there always exists a triplet $(r_{\mathrm{assoc}}, r_{\mathrm{poly}}, r_{\mathrm{nick}})$ satisfying equation (3).

Construction: For example, taking (as gauge fixing) $r_{\mathrm{nick}}=1$,

$$r_{\mathrm{assoc}}=\frac{b'}{b},\qquad r_{\mathrm{poly}}=\frac{k_1'}{k_1}\cdot\frac{b}{b'}$$

satisfies equation (3). Therefore, any $(k_1', b')$ can be represented by the three-parameter decomposition. □

2.4 Gauge Freedom and Physical Interpretation

An important note: while the observables $(k_1, b)$ have 2 degrees of freedom, the three-parameter decomposition has 3 parameters, so there exists a one-dimensional gauge freedom. Specifically, scaling $(r_{\mathrm{assoc}}, r_{\mathrm{nick}})$ simultaneously by the same factor $\lambda$ leaves equation (3) invariant:

$$(r_{\mathrm{assoc}}, r_{\mathrm{poly}}, r_{\mathrm{nick}}) \to (\lambda r_{\mathrm{assoc}}, r_{\mathrm{poly}}, \lambda r_{\mathrm{nick}})$$

Therefore, it is mathematically impossible to uniquely determine the three multipliers from observational data alone. Nevertheless, the reason for adopting the three-parameter decomposition is clarity of physical interpretation. By separating the three elementary processes of association $K_a$, turnover $k$, and saturation $K$, we can discuss mechanistically "what changed." Mathematically, it is possible to reduce to a two-parameter decomposition by setting $r_{\mathrm{nick}}=1$, but for interpreting terminal modifications and reaction condition changes, explicitly preserving the three elementary processes is far more transparent.

2.5 Propagation to the Entire System: Transformation Rules for Dimensionless Parameters $(g, \beta)$

In the two-variable model (reduced prey-predator system), the system's behavior is characterized by two dimensionless parameters:

$$g=\frac{k_1G}{k_2K_{m,P}},\quad \beta=\frac{b \cdot k_2K_{m,P}^2}{k_1}$$

Here, we assume that only the prey growth term parameters $(k_1, b)$ change, while predation and degradation term parameters $(k_2, K_{m,P})$ remain fixed. Substituting the three-parameter decomposition (equation (3)), we obtain:

$$g' = g\cdot\frac{r_{\mathrm{assoc}} \cdot r_{\mathrm{poly}}}{r_{\mathrm{nick}}},\qquad \beta' = \beta\cdot\frac{1}{r_{\mathrm{poly}}} \tag{4}$$

From these transformation rules, we can immediately see how each of the three elementary processes affects oscillation characteristics. When only association changes ($r_{\mathrm{assoc}}\neq1$, $r_{\mathrm{poly}}=r_{\mathrm{nick}}=1$), $\beta$ remains invariant and only $g$ scales. This corresponds to a parallel shift in effective template concentration $G$, shifting the threshold for oscillation emergence. When only turnover changes ($r_{\mathrm{poly}}\neq1$), $\beta$ is rescaled inversely, changing the bifurcation structure itself. When only saturation changes ($r_{\mathrm{nick}}\neq1$), $\beta$ remains invariant and $g$ scales in the opposite direction.

Oscillation characteristics such as period and amplitude are determined by the position in these two dimensionless quantities $(g, \beta)$. Therefore, the three-parameter decomposition provides a complete framework for tracking changes in oscillation characteristics at the elementary process level.

Furthermore, the apparent invariant $\kappa = k_1/b$ transforms as:

$$\kappa'=\frac{k_1'}{b'}=\frac{k_1}{b}\cdot r_{\mathrm{poly}}$$

In other words, $\kappa$ remains invariant even if association or saturation changes, and only turnover moves $\kappa$. By measuring $\kappa'/\kappa$ in experiments or simulations, we can distinguish what changed at the elementary process level.

3.1 Governing Equations

The system is described by the coupled ordinary differential equations:

$$\frac{dN}{dt} = k_1 \cdot \mathrm{pol} \cdot G \cdot \frac{N}{1+b \cdot G \cdot N} - k_2 \cdot \mathrm{pol} \cdot N \cdot P - \mathrm{rec} \cdot k_N \cdot \frac{N}{1+\frac{P}{K_{m,P}}}$$ $$\frac{dP}{dt} = k_2 \cdot \mathrm{pol} \cdot N \cdot P - \mathrm{rec} \cdot k_P \cdot \frac{P}{1+\frac{P}{K_{m,P}}}$$

3.2 Fourth-Order Runge-Kutta Method (RK4)

We employ the classical RK4 method for time integration, which provides fourth-order accuracy with excellent stability properties for this system. For each time step $\Delta t$:

$$\mathbf{y}_{n+1} = \mathbf{y}_n + \frac{\Delta t}{6}\left(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4\right)$$

where $\mathbf{y} = (N, P)^T$ and the slope vectors are:

$$\begin{aligned} \mathbf{k}_1 &= f(\mathbf{y}_n, t_n) \\ \mathbf{k}_2 &= f(\mathbf{y}_n + \frac{\Delta t}{2}\mathbf{k}_1, t_n + \frac{\Delta t}{2}) \\ \mathbf{k}_3 &= f(\mathbf{y}_n + \frac{\Delta t}{2}\mathbf{k}_2, t_n + \frac{\Delta t}{2}) \\ \mathbf{k}_4 &= f(\mathbf{y}_n + \Delta t\mathbf{k}_3, t_n + \Delta t) \end{aligned}$$

3.3 Rust/WebAssembly Implementation

The core integration loop is implemented in Rust and compiled to WebAssembly for high performance in the browser. Key implementation details:

3.4 Numerical Accuracy and Stability

The RK4 method with Δt = 0.5 min provides excellent accuracy for this system:

• Local truncation error: $O(\Delta t^5)$
• Global error: $O(\Delta t^4)$
• Stability region: RK4 is stable for complex eigenvalues with $|\lambda \Delta t| \lesssim 2.78$, well-satisfied for typical oscillation periods (~300 min)

For the heatmap generation, we use Δt = 0.1 min to ensure high accuracy across the entire parameter space, including regions with shorter periods or transient dynamics.

3.5 Period Detection Algorithm

To extract the oscillation period from time series data, we employ a peak-detection algorithm:

1. Identify local maxima: Find indices where N[i] > N[i-1] and N[i] > N[i+1]
2. Filter tail region: Use only peaks in the final 40% of the simulation (to exclude transient behavior)
3. Compute inter-peak intervals: Δt_i = t[peak[i+1]] - t[peak[i]]
4. Average period: T = mean(Δt_i)

If fewer than 2 peaks are detected, the period is reported as NaN (indicating non-oscillatory behavior). This algorithm is robust for sustained oscillations but can misidentify damped or chaotic dynamics.

4.1 3D Parameter Space Sampling

The heatmap videos visualize oscillation period as a function of three fold-change ratios:

Total computational cost: 340 × 340 × 150 = 17,340,000 individual simulations, each running for 3000 min with Δt = 0.1 min (30,000 integration steps).

4.2 Parallel Computation with Web Workers

To achieve reasonable computation time, the workload is distributed across multiple Web Workers:

4.3 Memory Management and IndexedDB Streaming

Large 3D parameter sweeps (340×340×150 = 17.3M simulations) require significant memory capacity. Each simulation involves RK4 integration with intermediate arrays (k₁, k₂, k₃, k₄) and time-series storage (30,000 data points per simulation at Δt = 0.1 min). To prevent browser crashes, we employ a streaming architecture:

• On-demand cell generation: Cells are generated as an iterator, never allocating the full 17.3M array
• Frame-by-frame accumulation: As simulations complete, results are accumulated into 2D grids
• IndexedDB storage: Completed frames are immediately written to disk (IndexedDB) and cleared from RAM
• Video rendering: Frames are loaded from IndexedDB one at a time during video encoding

This reduces peak memory usage significantly, enabling large-scale parameter sweeps on standard hardware.

4.4 Turbo Colormap

Period values are mapped to colors using the Turbo colormap (Google AI), which provides:

• Perceptual uniformity (equal perceived color differences for equal data differences)
• High contrast and distinguishability
• Colorblind-friendly design

The mapping is defined by polynomial approximations:

$$\begin{aligned} r(t) &= 34.61 + t(1172.33 - t(...))\quad \text{(clamped to [0, 255])} \\ g(t) &= 23.31 + t(557.33 + t(...)) \\ b(t) &= 27.2 + t(3211.1 - t(...)) \end{aligned}$$

where $t \in [0,1]$ is the normalized period value. Blue represents short periods, transitioning through cyan, green, yellow, to red for long periods.

4.5 Video Encoding

After all frames are computed, they are rendered to a <canvas> element and encoded to video. Two encoding options are available:

Both formats support configurable video duration (1–60 seconds) and produce videos that can be played directly in the browser or downloaded for offline viewing.

1. Montagne, K., Plasson, R., Sakai, Y., Fujii, T., & Rondelez, Y. (2011). Programming an in vitro DNA oscillator using a molecular networking strategy. Molecular Systems Biology, 7(1), 466.
2. Supplementary Information: Section S3 (Mathematical Model) and Table S5 (Optimized Parameters for PP₁ System).
3. Rust Programming Language: https://www.rust-lang.org/
4. WebAssembly: https://webassembly.org/
5. Turbo Colormap (Google AI): https://ai.googleblog.com/2019/08/turbo-improved-rainbow-colormap-for.html
6. MDN Web Docs - IndexedDB API: https://developer.mozilla.org/en-US/docs/Web/API/IndexedDB_API
7. MDN Web Docs - Web Workers API: https://developer.mozilla.org/en-US/docs/Web/API/Web_Workers_API