INTRODUCTION

Rumino utilizes the toehold-mediated strand displacement (TMSD) reaction to detect viral RNA. An exposed toehold on a DNA substrate (S) strand binds a complementary RNA or DNA target (T) strand, initiating branch migration and displacing a pre-hybridized incumbent (I) strand functionalized with a hydrophobic molecule, driving our wettability-based readout.

The toehold domain controls the forward binding rate, and the branch-migration domain controls displacement completion via sequence-dependent effects such as mismatch placement. Because each step depends on sequence, domain length, and assay conditions, we built a modelling framework to predict reaction behaviour before synthesis.

Our modelling has three parts:

  • Law of Mass Action and Kinetics Modelling: We model TMSD with a two-step mass-action scheme (target binding then branch migration) and solve the resulting ODEs in MATLAB. The same scheme is implemented in Visual DSD under deterministic and stochastic simulations to ensure consistency of results.
  • Sequence-specific design and analysis: Using a conserved 30nt region of HPAI H5N1 2.3.4.4b, we design S/I/T strands (7nt toehold, 23nt branch migration), validate thermodynamics with NUPACK, and verify mechanisms with coarse-grained simulation via the oxDNA.org web server.
  • Amplification of strands and RNA/DNA hybrids: We evaluate enzyme-free amplification modules (strand-displacement cascade and hybridization chain reaction), and assess the feasibility of using RNA invaders in the TMSD reaction based on the literature.


LAW OF MASS ACTION AND KINETICS MODELLING

This section describes the kinetics of the toehold-mediated strand displacement (TMSD) reaction that controls the Rumino detector. Because the reaction mechanism depends only on strand design and assay conditions rather than on a specific pathogen, we built a generic model based on the Law of Mass Action. Ordinary differential equations (ODEs) implemented in MATLAB describe the time-dependent concentrations of all species, enabling parameter sweeps and sensitivity analysis. Visual DSD simulations complement the ODE model by generating the reaction network from domain-level specifications and providing both deterministic and stochastic solvers to validate kinetics and noise effects before experimental testing.

Diagram to visualize the TMSD reaction in its stages:

Stages of TMSD figure

[1]

Figure 1. Stages of TMSD: An invader (target) binds the substrate's toehold, competes with the incumbent strand through branch migration, and displaces it once a few base pairs remain. Arrowheads indicate 3' ends of strands.

*Note: τ represents the toehold domain and β represents the branch migration domain.

Reaction Scheme (τ = 7nt; β = 23nt):

We model the Rumino TMSD gate as two reversible steps: toehold binding and branch migration/incumbent release, to determine the relationship between the involved complexes and the kinetics during the reaction.

Substrate and toehold migration figure

Species: S (substrate with exposed toehold (τ) and branch migration domain (β)), I (incumbent bound on β), T (target), T-S-I (toehold-captured intermediary complex), S-T (product duplex).
Reasoning: The toehold controls association and displays rate saturation at ~6-7 nucleotides, and the branch-migration step depends on the β-sequence and energy landscape (mismatch placement can lead to significant slowdowns) [2,3].

Law of Mass Action Representation in an ODE System:

Law of mass action representation

Assumptions

  • Isothermal, well-mixed, 25°C, 1M NaCl conditions.
  • The toehold on S is exposed (validated on a sequence-specific basis via NUPACK)[4].
  • Lower G-C content in the early branch migration domain, and higher G-C content in the late branch migration domain [1].

Parameter Table

Parameter Value Parameter Value Source
kbind Rate of toehold association 1.0x107 M-1⋅s-1 [5]
kunbind Rate of toehold disassociation 10 s-1 [6]
kbm_f Rate of forward branch migration 1.2 s-1 [7]
kbm_rev Rate of reverse branch migration 0.05 M-1⋅s-1 [6]
[T]0 Initial concentration 25nM [7]
[S-I]0 Initial concentration 25nM [7]
[T-S-I]0 Initial concentration 0 -
[S-T]0 Initial concentration 0 -
[I]0 Initial concentration 0 -

Table 1: Initial conditions and kinetic constants describing the relationship between complexes within the TMSD reaction.

Limitations

The kinetic constants in Table 1 are informed by reported measurements of toehold-mediated strand displacement and branch migration [3,5,6,7], with on-rates for 6-8 nt toeholds in the 10⁶-10⁷ M⁻¹·s⁻¹ range and branch migration step rates of 0.1-10 s⁻¹. Reverse rates are typically inferred from thermodynamic bias rather than measured directly, though our estimates fit realistic and observed ranges [6,8]. These values, however, are not absolute and represent a significant oversimplification of current models describing kinetics in TMSD. TMSD kinetics are strongly dependent on the conformational free energy landscape of the displacement complex, where sequence distribution and local structural stability introduce barriers or accelerants to branch migration [1]. As such, two systems of comparable toehold length and GC content may still display divergent rates.

For our design (7-nt toehold, 23-nt branch migration domain, ~50% GC, 10-50 nM strands), the constants represent literature-grounded estimates for simulation. They provide an internally consistent framework for modelling, but necessarily simplify the sequence-specific and structural nuances that influence displacement kinetics.

Solved equations in MATLAB:

Solved equations in MATLAB Solved equations in MATLAB Solved equations in MATLAB
Solved equations in MATLAB

Figure 2: Solved ODE system describing the dynamics of the TMSD reaction.

Solving the system of ODEs under the baseline parameters (7-nt toehold, 23-nt branch domain, 25nM input strands) shows that the TMSD reaction proceeds as designed. [T] and [SI] decrease in parallel as the target binds the S-I complex, forming the T-S-I complex. [T-S-I] rises rapidly, peaks within seconds, and then decays as branch migration displaces the incumbent. [ST] and [I] rise steadily with the decreasing T-S-I complex, reaching stagnatory levels within ~5-10 minutes. This system of solved equations confirms that the model captures both the present intermediate and the desired end-product.

VisualDSD Simulations

The Visual DSD simulations complement the ODE model by generating the reaction network from domain-level specifications. Using it alongside our ODE model ensures that our kinetic scheme matches the actual DNA design, and tests noise effects at realistic molecule counts before moving to experimental work [9].

Visual DSD deterministic solver

Figure 3: Time evolution of all species computed with the Visual DSD deterministic solver under baseline conditions (7-nt toehold, 23-nt branch domain, 25nM strands). Target (T) and Gate (S-I) are consumed symmetrically; the intermediate complex (T-S-I) peaks briefly; Product (S-T) and Blocker (I) accumulate to the expected stagnated concentration within minutes.



Gillespie simulation

Figure 4: Gillespie simulation describing the same system. Despite molecule-level fluctuations at nanomolar concentrations, the mean behaviour follows the deterministic prediction, confirming correct reaction mapping and mass conservation.

The deterministic Visual DSD run replicates the kinetic results predicted by the ODE model: a short-lived intermediate complex, rapid branch migration, and accumulation of [S-T] and [I] to predicted levels. The stochastic simulation displays the same overall trajectory with visible fluctuations, reflecting the inherent noise at nanomolar copy numbers. These outputs confirm that the chosen reaction pathway and rate constants produce the intended timescale and yield of the Rumino TMSD assay.



SEQUENCE-SPECIFIC DESIGN AND ANALYSIS

Goal: Define a pipeline for the sequence design of strands to be used in the TMSD reaction from sequenced genomic data in the Rumino detector. While the kinetic model shown above represents a general model for the TMSD reaction under assay conditions, this is a method to output functional disease-specific strands to be expressed and used in the wet lab through a cross-platform and simulation-based approach.

This pipeline follows from genome to designed sequence:
genome collection → conserved-region selection → sequence specificity assessment → thermodynamic verification → kinetic optimization → coarse-grained MD simulation → functional sequence sent to wet-lab for implementation.


Genome collection and selection of conserved regions (BV-BRC, Rumiseq software)

Goal: Retrieve complete genomes for the target sequence, ensuring that the final sequence is highly conserved and non-mutating.

Method: Download sequences from BV-BRC; filter by host, geography and recency (outbreak window). Complete multiple sequence alignment and track consensus identity at each position to identify high-confidence regions [10].

Acceptance: At least 30-nt windows with ≥90% conservation across all nucleotides.

Output: Ranked list of conserved regions with segment coordinates and sequence makeup.

H5N1 HPAI (2.3.4.4b) notes: Parameters used included samples sequenced in 2023-25, North America, Avian, Chicken, H5, segment 4 (HA gene). Though highly conserved regions were located within the recorded parameters, we were limited to the region within the genome that distinguishes the HPAI from the LPAI variant.


Specificity of sequence assessment (BLASTn)

Goal: Exclude candidates with host or environmental near-matches.

Method: BLASTn each candidate against human/avian/swine, microbiome and environmental genomic databases. Use megablast for close homology and, if needed, more sensitive blastn settings to catch near-matches [11].

Acceptance:

  • No off-target alignments with ≥80% identity over ≥24/30nt to non-target genomes.
  • No high-identity hits within the τ domain (7nt) and nearby β domain to host transcripts.

Output: BLAST report with pass/fail results and retained sequences.

H5N1 HPAI (2.3.4.4b) notes: No off-target binding was present.


Secondary structure analysis (NUPACK Utilities)

Goal: Ensure the τ domain is exposed and strands are free of structures that would block capture.

Method: NUPACK Utilities (single-complex ensemble) at assay temperature/salt for each strand: predict MFE structures and base-pairing probabilities; flag hairpins/self-dimers spanning τ or introducing strong structure in early β [12].

Acceptance: Unpaired τ domain; no unexpected, stable hairpins.

Output: Equilibrium probability and free energy display alongside a visual of the specified structure.

Gillespie simulation

Figure 5: Secondary structure analysis completed for the TMSD gate (S-I) in an HPAI H5N1 (2.3.4.3b) conserved sequence.


Thermodynamic analysis of TMSD (NUPACK Design)

Goal: Generate the design for the complete TMSD reaction with specified complementarity and domain parameters that form the intended complexes. Assess if the reaction is thermodynamically feasible.

Method: NUPACK Design (3 target tubes: reactants, products, global crosstalk; 2 maximum complex strands in each tube) at assay temperature/salt for each strand: assess if the experimental design is within the stop condition [12].

Acceptance: Presence of defects (failure in the reaction design) is below the stop condition of 0.05.

Output: Augmented objective function and multi-tube ensemble defect of the design trial, indicating susceptibility of the design to structural failure.

Structural stability design

Figure 6: Structural stability design analysis completed for the TMSD reaction.


Kinetic optimization (NUPACK Analysis; see Software page)

Goal: Adjust substrate-incumbent stability to bias displacement kinetics without losing specificity.

Method: Introduce single-base mismatches on intervals (“interval nucleotide switch”) and re-run NUPACK Test-Tube Analysis at assay conditions to quantify changes in yields and free energies [12].

Acceptance: S-T predicted yield ≥95 % with leak ≤5 %; improved forward bias; no new off-targets.

Output: Consistent analysis plots between unswitched and switched designs. Design and utilities pages show differing structural stability and ΔG°, displayed alongside the specified variant.

Secondary structure analysis

Figure 7: Secondary structure analysis completed for the TMSD gate (S-I) in an HPAI H5N1 (2.3.4.3b) conserved sequence with 7nt-interval switching.



Structural stability design analysis

Figure 8: Structural stability design analysis completed for the TMSD reaction with 7nt-interval switching.


Coarse-grained MD simulation (oxDNA)

oxDNA is a coarse-grained molecular dynamics (MD) framework used to simulate the structural behaviour of nucleic acids. Each nucleotide is represented as a single unit that can stack, pair, and move according to simplified physical rules derived from Watson-Crick interactions. This allows large-scale systems to be simulated efficiently while preserving realistic structural features such as strand alignment, fraying, and duplex stability [13].

For the Rumino system, oxDNA was used to verify that strand binding and displacement occur as designed. Simulations confirmed correct toehold recognition, stable product formation, and the absence of unintended structural intermediates. The resulting RMSD and energy profiles show a stable equilibrium, consistent with the reaction dynamics predicted by our kinetic model.

Goal: Confirm that the designed strands bind, branch-migrate and displace as intended under assay conditions.

Method: Simulate the TMSD system in oxDNA @, 25 °C; equilibrate (VMMC), then run MD; analyze bond occupancy, RMSF, duplex angle and energy profiles with oxDNA analysis tools [13].

Acceptance: Stable S-T product; I fully displaced; no long-lived misfolds; energies plateau.

Output: Final configurations, trajectory snapshots and analysis plots archived with the sequence and displayed alongside the simulated system.

RMSD analysis of the HPAI H5N1 TMSD

Figure 9: RMSD analysis of the HPAI H5N1 TMSD gate during oxDNA simulation. The root-mean-square deviation (RMSD) of the complex fluctuates within an acceptable range and stabilizes below 3nm (red line), indicating structural convergence under assay conditions.



Inter-strand distance measurements

Figure 10: Inter-strand distance measurements across the toehold and branch-migration regions during oxDNA simulation. Distances remain tightly clustered between 1-2nm, demonstrating stable base-pairing and alignment throughout the simulation.



Total energy plot

Figure 11: Total energy plot of the HPAI H5N1 TMSD gate during oxDNA simulation. The system rapidly equilibrates and maintains a consistent energy plateau, confirming that the simulated structure has reached a stable state under the chosen conditions.

Using this pipeline, we converted the HPAI H5N1 clade 2.3.4.4b genome into validated TMSD strands ready for synthesis and wet lab experimentation. Conservation mapping, BLAST screening, NUPACK verification, kinetic optimization, and oxDNA simulation ensured that each sequence is specific, accessible and predicted to function under our assay conditions. This section links the general kinetic model to disease-specific implementation and provides the wet lab with strands that have already passed sequence, thermodynamic and structural checks, reducing iteration and increasing reliability of the Rumino detector.



AMPLIFICATION OF STRANDS

A major challenge for the translation of the TMSD reaction into a practical diagnostic readout is amplification of input signals. In the Rumino detector, one target strand displaces a hydrophobic-modified incumbent strand to enable wettability-based readout. For sensitivity to reach clinically relevant limits, a single target input must trigger the release of multiple hydrophobic strands.


Strand Displacement Cascade (SDC)

A theoretical construct to achieve this is the strand displacement cascade, where a long “carrier” (C) strand is hybridized with multiple short, identical, non-hydrophobic incumbent strands [8]. Upon binding of the target strand to an exposed toehold, the invader displaces all incumbents in sequence. Each displaced short strand can then function as a new target in upstream reactions, triggering multiple displacement events. This design effectively amplifies one input into “x” outputs, where “x” is the number of pre-hybridized incumbent strands on the carrier. The kinetics of such cascades are shown as sequential branch migration steps and can be approximated as:


Total energy plot

Where [I] is the concentration of freely displaced incumbent strands, kbm the effective branch migration rate constant, and [T-C] the concentration of target-carrier complexes. Literature suggests that branch migration is nearly unbiased, with a per-step rate of ~10⁵ s⁻¹; thus, the overall cascade rate scales linearly with the number of base pairs displaced [14].


Hybridization chain reaction (HCR)

An alternative enzyme-free amplification method is the HCR, where stable hairpins are designed such that an invader opens the first hairpin, exposing a new single-stranded domain that triggers the opening of the next hairpin [15]. This results in a chain of strand displacement events and polymerization of long nicked double helices. HCR has been shown to achieve exponential amplification in vitro, with kinetics dependent on the stability of the hairpins and exposure of the toehold domains [15].

While HCR provides amplification, its integration into Rumino would require functionalization of each released strand with hydrophobic modifications, which introduces significant design constraints. For this reason, cascades are conceptually simpler for direct translation to wettability readouts.



RNA/DNA HYBRIDS IN TMSD

Recent experimental and computational studies [1,16,17] confirm that RNA/DNA hybrid strand displacement is feasible and often exhibits distinct kinetics compared to DNA/DNA systems.

Experimental kinetics.

Systematic measurements show that RNA invaders displacing DNA incumbents generally react more slowly than DNA invaders displacing DNA, due in part to the required conformational transition from B-form DNA to A-form RNA/DNA helices [16]. Importantly, kinetics are highly sequence dependent: for purine-rich invaders, RNA→DNA displacement can be more than 100-fold faster than DNA→RNA displacement, while the trend reverses in pyrimidine-rich contexts [1]. Measured rate constants for 6-nt toeholds fall in the 10⁵-10⁶ M⁻¹·s⁻¹ range, slower than the average DNA/DNA reactions, but still within diagnostic utility.

Coarse-grained modeling.

oxNA, a hybrid coarse-grained model built on oxDNA and oxRNA, was recently released and demonstrated that DNA-RNA hybrids adopt structural properties intermediate between A-form and B-form helices [17]. Importantly, free energy profiles confirmed that hybrid strand displacement is thermodynamically favourable, although base-pair asymmetry introduces sequence-dependent variability. For example, dA-rU pairs are less stable than dT-rA, and dC-rG shows reduced stability compared to dG-rC.

Implications for Rumino.

These findings confirm that designing DNA substrate/incumbent complexes with RNA targets is feasible. However, attention must be paid to:

  • Toehold placement (favouring 5' over 3') [16]
  • Toehold length (≥6nt for saturation) [16]
  • Lower G-C content in the early branch migration domain, as mismatches are more kinetically consequential near the toehold for DNA/RNA hybrids [1]

By incorporating these considerations, RNA targets such as avian influenza H5N1 can be reliably captured and amplified through DNA-based detectors, ensuring compatibility of the Rumino system with clinically relevant RNA biomarkers.



CONCLUSION

The Rumino framework combines sequence design, kinetic modelling, and structural simulation into a unified approach for predicting the behaviour of toehold-mediated strand displacement. By integrating ODE-based mass-action modelling, Visual DSD simulations, and oxDNA verification, we confirmed that the system behaves as expected under specified assay conditions and that each designed strand forms stable, functional complexes within a reasonable timescale. Future work will focus on incorporating sequence-dependent energies to improve the accuracy of rate predictions and extend the framework to multi-gated and RNA-targeted systems. These refinements will enhance the predictive power of Rumino and support its continued development as a generalizable platform for amplification-free, sequence-specific detection of viral targets.



REFERENCES

[1] Ratajczyk EJ, Bath J, Šulc P, Doye JPK, Louis AA, Turberfield AJ. 2025. Controlling DNA–RNA strand displacement kinetics with base distribution. Proc Natl Acad Sci U S A. 122(23):e2416988122. https://doi.org/10.1073/pnas.2416988122

[2] Šulc P, Ouldridge TE, Romano F, Doye JP, Louis AA. 2015. Modelling toehold-mediated RNA strand displacement. Biophys J. 108(5):1238–1247. https://doi.org/10.1016/j.bpj.2015.01.023

[3] Machinek R, Ouldridge T, Haley N, Bath J, Turberfield AJ. 2014. Programmable energy landscapes for kinetic control of DNA strand displacement. Nat Commun. 5:5324. https://doi.org/10.1038/ncomms6324

[4] Akay A, Reddy HN, Galloway R, Kozyra J, Jackson AW. 2024. Predicting DNA toehold-mediated strand displacement rate constants using a DNA-BERT transformer deep learning model. Heliyon. 10(7). https://doi.org/10.1016/j.heliyon.2024.e28443

[5] Zhang DY, Winfree E. 2009. Control of DNA strand displacement kinetics using toehold exchange. J Am Chem Soc. 131(47). https://doi.org/10.1021/ja906987s

[6] Srinivas N, Ouldridge TE, Sulc P, Schaeffer JM, Yurke B, Louis AA, Doye JP, Winfree E. 2013. On the biophysics and kinetics of toehold-mediated DNA strand displacement. Nucleic Acids Res. 41(22):10641–10658. https://doi.org/10.1093/nar/gkt801

[7] Broadwater DWB Jr, Cook AW, Kim HD. 2021. First passage time study of DNA strand displacement. Biophys J. 120(12):2400–2412. https://doi.org/10.1016/j.bpj.2021.01.043

[8] Qian L, Winfree E. 2011. Scaling up digital circuit computation with DNA strand displacement cascades. Science. 332(6034):1196–1201. https://doi.org/10.1126/science.1200520

[9] Lakin MR, Youssef S, Polo F, Emmott S, Phillips A. 2011. Visual DSD: a design and analysis tool for DNA strand displacement systems. Bioinformatics. 27(22):3211–3213. https://doi.org/10.1093/bioinformatics/btr543

[10] Olson RD, Assaf R, Brettin T, Conrad N, Cucinell C, Davis JJ, Dempsey DM, Dickerman A, Dietrich EM, Kenyon RW, et al. 2023. Introducing the Bacterial and Viral Bioinformatics Resource Center (BV-BRC): a resource combining PATRIC, IRD and ViPR. Nucleic Acids Res. 51(D1):D678–D689. https://doi.org/10.1093/nar/gkac1003

[11] Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL. 2009. BLAST+: architecture and applications. BMC Bioinformatics. 10:421. https://doi.org/10.1186/1471-2105-10-421

[12] Zadeh JN, Steenberg CD, Bois JS, Wolfe BR, Pierce MB, Khan AR, Dirks RM, Pierce NA. 2011. NUPACK: analysis and design of nucleic acid systems. J Comput Chem. 32(1):170–173. https://doi.org/10.1002/jcc.21596

[13] Poppleton E, Romero R, Mallya A, Rovigatti L, Šulc P. 2021. OxDNA.org: a public webserver for coarse-grained simulations of DNA and RNA nanostructures. Nucleic Acids Res. 49(W1):W491–W498. https://doi.org/10.1093/nar/gkab324

[14] Srinivas N, Ouldridge TE, Sulc P, Schaeffer JM, Yurke B, Louis AA, Doye JP, Winfree E. 2013. On the biophysics and kinetics of toehold-mediated DNA strand displacement. Nucleic Acids Res. 41(22):10641–10658. https://doi.org/10.1093/nar/gkt801

[15] Dirks RM, Pierce NA. 2004. Triggered amplification by hybridization chain reaction. Proc Natl Acad Sci U S A. 101(43):15275–15278. https://doi.org/10.1073/pnas.0407024101

[16] Liu H, Hong F, Smith F, Goertz J, Ouldridge T, Stevens MM, Yan H, Šulc P. 2021. Kinetics of RNA and RNA:DNA hybrid strand displacement. ACS Synth Biol. 10(11):3066–3073. https://doi.org/10.1021/acssynbio.1c00336

[17] Ratajczyk EJ, Šulc P, Turberfield AJ, Doye JPK, Louis AA. 2024. Coarse-grained modeling of DNA–RNA hybrids. J Chem Phys. 160:115101. https://doi.org/10.1063/5.0199558