YOKABIO

Natural and biological phenomena, despite arising from complex molecular reactions, often exhibit universal mathematical structures. In DNA concentration oscillator systems, periodic oscillations emerge from predator-prey-like interactions, and a clear mathematical structure underlies this behavior. By examining the system through mathematical models, we can understand essential mechanisms that are difficult to observe experimentally, as well as the system's response to changes in conditions.

This system is described by the following system of coupled ordinary differential equations:

$$\frac{dN}{dt} = \underbrace{k_1 \cdot \mathrm{pol} \cdot G \cdot \frac{N}{1+b \cdot G \cdot N}}_{\text{Prey growth}} - \underbrace{k_2 \cdot \mathrm{pol} \cdot N \cdot P}_{\text{Consumption by predation}} - \underbrace{\mathrm{rec} \cdot k_N \cdot \frac{N}{1+\frac{P}{K_{m,P}}}}_{\text{Degradation}}$$ $$\frac{dP}{dt} = \underbrace{k_2 \cdot \mathrm{pol} \cdot N \cdot P}_{\text{Predator growth}} - \underbrace{\mathrm{rec} \cdot k_P \cdot \frac{P}{1+\frac{P}{K_{m,P}}}}_{\text{Degradation}}$$

Here, $N$ represents the concentration of prey DNA, and $P$ represents the concentration of predator DNA. The parameters $k_1, k_2$ are reaction rate constants, $b$ is a constant representing saturation effects, $\mathrm{pol}, \mathrm{rec}$ are the concentrations of polymerase and recombinase, respectively, $G$ is the template DNA concentration, and $K_{m,P}$ is the Michaelis constant. The Michaelis-Menten-like saturation appearing in the prey growth term arises from competition between template binding and the nicking process. For detailed derivation, refer to the Supporting Information of the original paper.

Focus on the Prey Growth Term

The reactions in the PP oscillator system can be broadly classified into three categories: prey (N) growth, predator (P) increase, and degradation (inactivation/loss of components). This analysis focuses on the prey growth term. The reason is mechanistically clear. The terminal amino acid modification most directly affects the chemical processes of template G and N binding (G:N association) and the turnover of the synthesis reaction proceeding from this scaffold, both of which correspond to internal mechanisms of the prey growth term.

In contrast, predator increase and degradation processes are not primarily linked to terminal modification effects, and allowing these to vary simultaneously introduces excessive degrees of freedom, making it difficult to identify causal relationships. By focusing on the single causal chain of modification → prey growth elementary processes → oscillation characteristics, the logical flow of the argument becomes clear.

This restriction has another advantage. The prey growth term can effectively be described by three elementary factors: "strength of G:N association," "synthesis turnover rate," and "saturation characteristics derived from the nicking process." If we express the effects of terminal amino acids as multipliers of these three factors, the apparent parameters ($k_1$ and $b$) can be uniquely reconstructed from the upstream elementary factors. Rather than having users independently vary $k_1$ and $b$, by providing only physical measures such as "how much the modification changed association or turnover," the model automatically maintains consistency while determining behavior.

In particular, when only association changes, the effective scale of $G$ changes (the oscillation emergence threshold shifts in parallel), and when turnover decreases, the period tends to lengthen (the dimensionless parameter $\beta$ is rescaled). The phenomenon of period elongation observed experimentally in the order G → G-amine → G-lysine can potentially be explained without parameter fitting by determining which elementary factors change in which direction.

Below, we introduce a framework that decomposes the prey growth term into three fundamental parameters (hereafter referred to as three-parameter decomposition) and clarify its mathematical properties.

Three-Parameter Decomposition Representation of the Prey Growth Term

Following the derivation in the Supporting Information of the original paper, the flux of prey growth can be written as:

$$\varphi_{N\triangleright N} = k \cdot \mathrm{pol}\cdot\frac{[G{:}N]}{K+[G{:}N]},\qquad [G{:}N]=K_a^{GN} G_{\mathrm{free}} N_{\mathrm{free}}$$

Here, $[G{:}N]$ is the concentration of the template G and prey N complex, $K_a^{GN}$ is its association constant, $k$ is the turnover rate of the synthesis reaction, and $K$ is the saturation constant derived from the nicking process. Under the usual low-occupancy approximation ($K_a^{GN} \cdot G \cdot N \ll 1$), we can assume $G_{\mathrm{free}}\simeq G$ and $N_{\mathrm{free}}\simeq N$, giving:

$$\varphi_{N\triangleright N} = k \cdot \mathrm{pol}\cdot\frac{K_a^{GN}GN}{K+K_a^{GN}GN} = k_1 \cdot \mathrm{pol}\cdot\frac{GN}{1+b \cdot GN}$$

where

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

G terminal amino acid modification and changes in reaction buffer (salt concentration, temperature, dyes, SSB, etc.) can alter the association strength $K_a^{GN}$, effective synthesis turnover rate $k$, and saturation constant $K$ derived from nicking. We represent any condition change by the following three multipliers:

$$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 represent values after the change.

The physical meaning of each multiplier is as follows. The association multiplier $r_{\mathrm{assoc}}$ is an index representing how much the binding free energy of the G:N partial duplex changes due to modification, and can be written as:

$$r_{\mathrm{assoc}} = \frac{K_a'}{K_a} = \exp\left(-\frac{\Delta\Delta G_{\mathrm{assoc}}}{RT}\right)$$

(where $\Delta\Delta G_{\mathrm{assoc}}=\Delta G'_{\mathrm{assoc}}-\Delta G_{\mathrm{assoc}}$). This is the standard Boltzmann factor form, condensing the effects of terminal amino acid charge, aromaticity, hydrophobicity, linker length, as well as temperature, salt, dyes, SSB, etc., on strengthening or weakening the G and N association equilibrium. When the template hairpin opening fraction $f_{\mathrm{open}}$ is relevant, it effectively becomes $K_a \to f_{\mathrm{open}}K_a$, which can be absorbed as $r_{\mathrm{assoc}} \leftarrow r_{\mathrm{assoc}} \cdot f_{\mathrm{open}}$.

The turnover multiplier $r_{\mathrm{poly}}$ is the change in the effective turnover rate $k$ derived from polymerase extension and subsequent processing. Here, $k$ represents the effective constant for the upper limit rate per polymerase molecule and is sensitive to temperature, enzyme activity, local steric hindrance, charge environment (changes in surrounding electric field due to terminal modification), etc. $r_{\mathrm{poly}}<1$ means delay in the extension cycle (decrease in upper limit rate), while $r_{\mathrm{poly}}>1$ means acceleration.

The saturation multiplier $r_{\mathrm{nick}}$ represents the change in the Michaelis-like saturation constant $K$ (half-saturation threshold) arising from the nicking process. $K$ is a threshold determining "at what level of $[G{:}N]$ accumulation does the rate plateau," and varies with the effective affinity of the nicking enzyme, product stability, and local structural changes (mechanical properties and dissociation characteristics of the duplex accompanying terminal modification). $r_{\mathrm{nick}}>1$ indicates that saturation becomes more difficult (plateau only at larger $[G{:}N]$), while $r_{\mathrm{nick}}<1$ indicates easier saturation.

Using these three multipliers, the apparent parameters in equation (1) are transformed as follows:

$$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}$$

Hereafter, we refer to this representation as the three-parameter decomposition.

Application to Experimental Observations

The three-parameter decomposition framework provides a powerful tool for interpreting experimental results. As described in our Experiment page, we observed that G-amine modification (adding an amino group to the template 3' terminus) causes oscillation period elongation compared to unmodified G template.

Initially, we hypothesized that the additional cationic charge from amine protonation would enhance G:N hybridization, corresponding to $r_{\mathrm{assoc}} > 1$ (increased association). If this were true, our simulator predicts that the period should shorten (try adjusting the association slider above 1.0). However, the experimental result showed the opposite—the period became longer.

Within the three-parameter framework, this observation is consistent with $r_{\mathrm{assoc}} < 1$ (decreased G:N association). The amine modification may actually weaken the G:N binding, possibly due to altered local structure or electrostatic environment that interferes with hybridization. You can verify this hypothesis using the interactive simulator below—decreasing the association scaling parameter ($r_{\mathrm{assoc}} < 1$) indeed lengthens the oscillation period, matching our experimental observations.

This example demonstrates how the three-parameter decomposition enables us to distinguish between competing mechanistic hypotheses by connecting experimental observations to specific elementary processes. For the complete experimental discussion, including results with Boc-protected lysine modification, please see the Experiment page.

Note: For comprehensive technical details, complete mathematical derivations, implementation specifications, and extended analysis, please refer to the Technical Appendix.

1. Interactive Simulator

We have prepared an interactive simulator that implements the mathematical model described above. This simulator allows you to observe the system's behavior in real-time when parameters are applied to the PP oscillator equations. The left graph displays the time-series concentration changes of prey (N) and predator (P), while the right phase portrait shows the trajectory (limit cycle) in N-P space.

Please try it out. By adjusting three parameters using sliders—template concentration (G), polymerase concentration (pol), and recombinase concentration (rec)—you can intuitively understand how oscillation period and amplitude change. For example, increasing G lengthens the oscillation period, while excessive increases in rec cause the oscillations to dampen.

For simplicity, this demonstration simulator limits the adjustable parameters to three. Other parameters that cannot be adjusted ($k_1, k_2, k_N, k_P, b, K_{m,P}$, etc.) are set to standard values optimized for the PP₁ system through fitting in the original paper's Supporting Information (Table S5). Specifically, the values are: $k_1 = 0.002$ [nM⁻¹min⁻¹], $k_2 = 0.0031$ [nM⁻¹min⁻¹], $k_N = 0.021$ [min⁻¹], $k_P = 0.0047$ [min⁻¹], $b = 4.8 \times 10^{-5}$ [nM⁻²], $K_{m,P} = 34$ [nM].

For more detailed technical specifications, complete parameter descriptions, and implementation details, please refer to the Technical Appendix. Through this simulator, you can intuitively understand the dynamical properties of DNA concentration oscillator systems and their response to parameter changes.

2. Parameter Space Heatmap Analysis

To systematically investigate the influence of the scaling parameters ($r_{\mathrm{assoc}}, r_{\mathrm{poly}}, r_{\mathrm{nick}}$) introduced in the three-parameter decomposition, we performed heatmap analysis of the 3D parameter space. The objective is to visualize how changes in each elementary process due to terminal amino acid modifications or reaction conditions affect oscillation period and amplitude, and to establish correspondence with experimental results.

Specifically, we created 3D heatmaps by assigning two parameters to the X-Y plane and allocating the remaining parameter to the time axis (T-axis). At each grid point, a complete oscillation simulation was executed, and the oscillation period in the final 40% interval was extracted as the evaluation metric. This yields a "topographic map of oscillation period" in the three-parameter combination space.

From the same dataset, we generated three videos with different X-Y-T axis assignments. Each video represents the same information from a different perspective, facilitating visual understanding:

• Video 1 (X: $r_{\mathrm{assoc}}$, Y: $r_{\mathrm{poly}}$, T: $r_{\mathrm{nick}}$) — Observe the interaction between association and turnover while viewing saturation effects as temporal progression
• Video 2 (X: $r_{\mathrm{nick}}$, Y: $r_{\mathrm{assoc}}$, T: $r_{\mathrm{poly}}$) — Observe the interaction between saturation and association while viewing turnover effects as temporal progression
• Video 3 (X: $r_{\mathrm{poly}}$, Y: $r_{\mathrm{nick}}$, T: $r_{\mathrm{assoc}}$) — Observe the interaction between turnover and saturation while viewing association effects as temporal progression

Simulation Conditions

Each heatmap video was generated using the following computational conditions:

An independent oscillation simulation was executed at each grid point, and periods were automatically extracted from the resulting time-series data. Period detection employs a peak-finding algorithm, calculating the representative period by averaging time intervals between adjacent local maxima. In regions where oscillations are not observed (damping or divergence), values are treated as NaN (Not a Number) and displayed in gray on the heatmap.

Heatmap Videos

Below, we present heatmap videos created from three perspectives. In each video, the period distribution for combinations of two parameters assigned to the horizontal (X) and vertical (Y) axes is expressed by color, and as the time axis (T) progresses, the remaining parameter changes. Colors represent a gradient from blue (short period) to red (long period), using the Turbo colormap.

Interpretation of Results

From these heatmap videos, the following insights can be obtained:

• Dominance of turnover rate: Decreasing $r_{\mathrm{poly}}$ lengthens the period throughout the entire parameter space. This corresponds to the effect appearing as $\beta' = \beta / r_{\mathrm{poly}}$ in equation (4) and becomes a candidate major factor explaining the period elongation observed with G-lysine modification.
• Balance between association and saturation: Since the ratio of $r_{\mathrm{assoc}}$ to $r_{\mathrm{nick}}$ determines $b$, the balance between these two strongly governs the shape of oscillation patterns. When both change by the same scaling factor (gauge transformation), the period distribution remains nearly invariant.
• Variation of oscillation emergence region: Increasing association strength tends to expand the parameter region where oscillations emerge, while decreasing it suppresses oscillations. This corresponds to the phenomenon where amino acid modifications experimentally alter oscillation characteristics significantly.
• Nonlinear response: Period contour lines trace complex curves rather than straight lines, indicating that changes in elementary processes are reflected in oscillation characteristics in a non-trivial manner. This nonlinearity reflects the richness of the system's bifurcation structure.

In the future, once experimental data under different modification conditions can be successfully fitted to this model, plotting the resulting period data on these heatmaps will enable reverse estimation of how each modification affects the three elementary processes. By systematically combining experiments with varied reaction conditions (temperature, salt concentration, etc.) and performing parameter fitting, more precise separation and identification of each multiplier effect is expected.