# nNxB: a new coarse-grained model for RNA and DNA nanotechnology

F. Tosti Guerra  Department of Physics, Sapienza University of Rome  E.
Poppleton  School of Molecular Sciences and Center for Molecular Design and
Biomimetics, The Biodesign Institute, Arizona State University  Biophysical
Engineering Group, Max Planck Institue for Medical Research  P. Šulc  School
of Molecular Sciences and Center for Molecular Design and Biomimetics, The
Biodesign Institute, Arizona State University  Department of Bioscience,
School of Natural Sciences, Technical University Munich  L. Rovigatti
Department of Physics, Sapienza University of Rome

(August 28, 2024)

###### Abstract

> The folding processes of RNA and DNA play crucial roles in biological
> systems and nanotechnology. However, studying these processes with high-
> resolution models is beyond current computational capabilities. In this
> article, we present a new coarse-grained model for investigating the folding
> dynamics of nucleic acids. Our model represents $n$ nucleotides using a
> single patchy particle and is parametrized using well-established nearest-
> neighbor models. By simplifying the system and interactions, our model
> allows for simulations at timescales and length scales that are currently
> inaccessible to more detailed models. To validate the performance of our
> model, we conducted extensive simulations of various single-stranded
> systems. We examined the thermodynamics of DNA hairpins, capturing their
> stability and structural transitions. We also investigated the folding of an
> MMTV pseudoknot, a complex RNA structure involved in viral replication.
> Furthermore, we explored the folding of a challenging RNA tile containing a
> k-type pseudoknot. The new coarse-grained model accurately captures the
> folding behavior of diverse nucleic acid structures, providing insights into
> their thermodynamics. The successful reproduction of experimental data and
> favorable comparisons with existing coarse-grained models validate the
> effectiveness of our approach.

††preprint: AIP/123-QED

##  I Introduction

RNA and DNA molecules play critical roles in biological systems and have
become increasingly important in constructing nanoscale architectures
Pitchiaya and Krishnan (2006); Dey  _et al._ (2021). Understanding the folding
processes of these molecules is essential for unraveling their functional
mechanisms and facilitating the design of innovative nanoscale structures,
with broad implications across scientific disciplines. One significant
application of understanding folding pathways is in RNA and DNA structure
prediction, since kinetic traps can prevent access to the conformation
corresponding to absolute free-energy minima. Accurately predicting secondary
and tertiary structures is crucial for elucidating their biological functions,
designing therapeutic interventions, engineering RNA/DNA-based nanomaterials
Tinoco Jr and Bustamante (1999), facilitating vaccine delivery Tang  _et al._
(2021), and improving drug delivery Li  _et al._ (2018).

Another area where knowledge of folding pathways is valuable is functional RNA
design. RNA molecules have diverse functional roles, including catalysis, gene
regulation, and molecular recognition Benenson (2009). Understanding folding
pathways enables the rational design of functional RNA structures with desired
properties, such as enhanced catalytic activity or improved binding affinity
to specific targets Du, Hu, and Ren (2023).

Finally, the use of RNA and DNA in nanotechnology is gaining traction, as they
serve as versatile building blocks for constructing a wide range of nanoscale
devices and materials Seeman (2010). These include self-assembling complex 2D
and 3D architectures using DNA origamiRothemund (2006); Douglas  _et al._
(2009); Dey  _et al._ (2021), dynamic structures like DNA motors Bath  _et
al._ (2009); Ouldridge  _et al._ (2013a); Shi  _et al._ (2023), and even DNA-
based computational systems Adleman (1994); Chen  _et al._ (2013); Wang  _et
al._ (2023). Exploring the intricate folding pathways of RNA and DNA helps
uncover fundamental principles of self-assembly, enabling the design of
increasingly complex and programmable nanostructures with exciting
possibilities for applications in various fields.

Computational models have been developed to explore the processes involved in
RNA and DNA folding. These models range from atomistic to coarse-grained
resolutions, each offering unique advantages and insights into the dynamics of
folding. The all-atom model is a commonly employed approach that explicitly
represents individual atoms and their interactions. Atomistic force fields
such as AMBER Cornell  _et al._ (1995) and CHARMM Brooks  _et al._ (1983) have
been extensively used to investigate the folding of nucleic acidsSponer  _et
al._ (2018). While these models accurately capture atomistic details and
provide valuable insights into the structural and energetic aspects of RNA and
DNA, their simulations are computationally demanding and limited in timescale.

To overcome the limitations of all-atom models, coarse-grained models have
emerged as powerful tools for studying longer lengths and timescales of
nucleic acid folding. Coarse-grained models simplify the representation of
nucleic acids by grouping multiple atoms into a single particle or bead,
allowing for more efficient simulations and exploration of extended
timescales. Several coarse-grained models have been specifically developed for
studying DNA and RNA Parisien and Major (2008); Savelyev and Papoian (2009);
Paliy, Melnik, and Shapiro (2010); Pasquali and Derreumaux (2010); Ouldridge,
Louis, and Doye (2011); Šulc  _et al._ (2012); Hinckley  _et al._ (2013);
Cragnolini, Derreumaux, and Pasquali (2013); Xia  _et al._ (2013); Denesyuk
and Thirumalai (2013); Šulc  _et al._ (2014); Korolev  _et al._ (2014);
Maciejczyk  _et al._ (2014); Maffeo  _et al._ (2014); Rovigatti  _et al._
(2015); Machado and Pantano (2015); Uusitalo  _et al._ (2015); Dans  _et al._
(2016); Ivani  _et al._ (2016); Chakraborty, Hori, and Thirumalai (2018);
Maffeo and Aksimentiev (2020); DeLuca  _et al._ (2023).

Coarse-grained models offer a reduced number of degrees of freedom compared to
atomistic approaches when representing DNA/RNA. They typically treat the
solvent and solvated ions implicitly and represent groups of atoms in the
DNA/RNA structure with effective interactions. This simplification enables the
study of larger and more extended molecular systems for longer times. Coarse-
graining techniques can be broadly categorized as either from a bottom-up or
top-down, each aiming to capture specific aspects of the system. The bottom-up
approach formally maps the statistical behavior of a more detailed model into
a coarse-grained description, while the top-down approach aims to reproduce as
many experimentally relevant properties as possible Sengar  _et al._ (2021).

Many of the models mentioned above are parametrized using nearest-neighbor
(NN) models. Initially introduced by Poland and Scheraga to investigate duplex
denaturation phase transitions Poland and Scheraga (1966), this approach has
been meticulously parameterized SantaLucia Jr and Hicks (2004); Allawi and
SantaLucia (1998a); Allawi and SantaLucia Jr (1998); Allawi and SantaLucia
(1998b, 1997); Peyret  _et al._ (1999); Allawi and SantaLucia (1998a); Turner
and Mathews (2010); Huguet  _et al._ (2010); Bae, Fang, and Zhang (2020);
Zuber  _et al._ (2022) to describe binding equilibria for oligonucleotides.

A nearest-neighbor model calculates the free energy change ($\Delta G$)
associated with forming a duplex by summing the contributions of individual
base pair interactions. These interactions are characterized by experimentally
determined enthalpy ($\Delta H_{ij}$) and entropy ($\Delta S_{ij}$) values for
each possible combination of adjacent base pairs $(i,j)$. Different base pair
combinations have different thermodynamic parameters, reflecting variations in
hydrogen bonding, stacking interactions, and other factors.

The concentration of the fully bound duplex product, denoted as $[AB]$, is
related to the thermodynamic parameters and the equilibrium constant ($K$) of
the duplex formation. The equilibrium constant is defined as the ratio of the
concentration of the product $[AB]$ to the concentrations of the individual
strands $[A]$ and $[B]$:

$K=\frac{[AB]}{[A][B]}\propto\exp{(-\beta(\Delta H-T\Delta S))},$ (1)

where $\Delta H$ and $\Delta S$ are the total enthalpy and entropy change upon
binding, and $\beta=1/k_{B}T$.

By considering nearest-neighbor interactions along the entire sequence length,
the model predicts the stability of duplexes and provides insights into their
melting temperatures, binding affinities, and overall thermodynamic
properties.

In our study, we introduce a new coarse-grained description of nucleic acids
where a single patchy particle represents $n$ nucleotides, enabling efficient
exploration of RNA and DNA folding processes. We parametrize our model using
well-established nearest-neighbor models SantaLucia Jr and Hicks (2004);
Allawi and SantaLucia (1998a); Allawi and SantaLucia Jr (1998); Allawi and
SantaLucia (1998b, 1997); Peyret  _et al._ (1999); Zuber  _et al._ (2022) to
capture essential interactions involved in folding dynamics while maintaining
computational efficiency. Our model combines the strengths of existing coarse-
grained and patchy particle models, providing new possibilities for
investigating folding phenomena at previously inaccessible time- and length-
scales. This first iteration of the model has been parametrized by considering
single-stranded systems only, and its thermodynamic performance has been
compared with available experimental and numerical data.

##  II Model description

(a)  |  (b)
---|---
Figure 1: LABEL:sub@subfig:_sketch schematic representation of the attractive
interaction in nNxB, LABEL:sub@subfig:_mapping the diameter of the DNA/RNA
helix is approximately 2 nm. Considering our beads with a diameter denoted as
$\sigma$ and the patches on the beads’ surface with a diameter of
$r_{p}=0.219\sigma$, we can determine that
$2\sigma+2\sigma(\delta_{pb}-0.5)+r_{p}\sigma=2$ nm. Therefore, we can deduce
that $\sigma=2/(1+2\delta_{pb}+r_{p})$ nm.

In our model, we introduce a polymer composed of N beads to represent the
DNA/RNA strand. Each bead is decorated with an attractive spherical patch,
having a radius of 0.219$\sigma$ and positioned $\delta_{pb}=0.65\sigma$ away
from the center of the bead, and represents $n$ nitrogenous bases that make up
the nucleic acid molecule. We incorporate several potentials to reproduce the
thermodynamic and mechanical characteristics of these biopolymers.

Figure 1 provides a visual summary of the notation and different interaction
terms used in our model, while the detailed functional forms are reported in
Appendix A. The potentials between topologically bonded beads linked through
the backbone include the Kremer-Grest force field Grest and Kremer (1986),
combining a Weeks-Chandler-Andersen (WCA) potential to model excluded volume,
and a finitely extensible nonlinear elastic (FENE) potential to mimic the
covalent bonds along the strand.

We also introduce a term to add stiffness to the structure, a three-body
potential ($V_{\rm semiflex}$) that tends to align three consecutive beads
which results in a different persistence length between single and double
strands (see II.3). To model the stacking of bases in DNA/RNA, we incorporate
a term called $V_{\rm stack}$. This term promotes the alignment of the
directions $\vec{a}_{1}$, which determine the patch positions, of consecutive
beads by using a cosine-angle potential.

For the non-bonded interactions, we used the WCA potential to account for
excluded-volume interactions. Additionally, we employ a patch-patch
interaction potential, based on the functional form proposed by Stillinger and
Weber Weber and Stillinger (1985), to model the hybridization of nucleotides.
The strength of the attraction between two beads $i$ and $j$, $\epsilon_{ij}$,
depends on the types of nucleotides described by each interacting bead and is
computed with a NN model as described below. To ensure that each bead can bind
to only one other bead, we implement a repulsive three-body interaction
($V_{\rm 3b}$) that penalizes the formation of triplets of bonded beads
Sciortino (2017). This repulsive potential compensates for the gain associated
with the formation of a second bond and can be tuned to favor bond swapping.
More, specifically, the parameter $\lambda$ in $V_{\rm 3b}$ allows to
interpolate between the limits of swapping ($\lambda=1$) and non-swapping
($\lambda\gg 1$) bonds. Finally, the $V_{\rm sticky}$ potential is modulated
by a term $V_{\rm direct}$ that takes into account the $\vec{a}_{3}$
orientations of the beads to ensure that only antiparallel strands can bind to
each other.

In the following, we run molecular dynamics simulations at constant
temperature Russo, Tartaglia, and Sciortino (2009) using the nNxB model
implemented in the oxDNA simulation engine Poppleton  _et al._ (2023). The
equations of motion are integrated with a velocity Verlet algorithm with a
time step $\Delta t=0.002$ (in internal units).

###  II.1 Parametrization

The strength of the sticky interaction between any two beads $i,j$ is
controlled by a term $\beta\varepsilon_{ij}$ that depends on the nucleotide
sequences present in each bead. This term models the pairing of the
nucleotides that compose the two beads by taking into account the Gibbs free
energy, multiplied by $\beta$, as provided by the nearest-neighbor model. In
the following we will consider systems simulated at a fixed monovalent salt
concentration of $0.5$ M, but any other condition can be considered, as long
as the NN model employed supports it.

In order to properly account for temperature variations, we need to separate
the contributions from enthalpy and entropy. This can be achieved by rewriting
$-\Delta G/k_{B}T$ as follows:

$\displaystyle-\frac{\Delta G}{k_{B}T}$ $\displaystyle=-\left(\frac{\Delta
H}{k_{B}T}-\frac{\Delta S}{k_{B}}\right)$ $\displaystyle=-\left(\frac{\Delta
H}{k_{B}T_{\rm ref}}\frac{k_{B}T_{\rm ref}}{k_{B}T}-\frac{\Delta
S}{k_{B}}\right)$ $\displaystyle=-\left(\mu\frac{\Delta H}{k_{B}T_{\rm
ref}}-\frac{\Delta S}{k_{B}}\right),$ (2)

where $\mu\equiv T_{\rm ref}/T$ is a parameter that quantifies the temperature
difference between the simulation temperature and the temperature $T_{\rm
ref}$ used to estimate the values of $\Delta H,\Delta S,\Delta G$ as reported
in NN models. Typically, as it is the case here, $T_{\rm ref}=37$ ∘C.

Figure 2: Contributions to the interaction strength between the two colored
(red) beads, illustrating the base steps involved. The blue boxes surround
base steps made of nucleotides within the beads, while the green boxes
surround the two boundary base steps (i.e. base steps comprising nucleotides
both inside and outside of the beads considered).

When calculating the $\Delta G$ between two beads, each bead should not be
considered in isolation but as part of a larger system, and the coarse-
graining procedure should also take into account the nucleotide base steps at
the boundaries of the beads. This means including the last nucleotide of the
previous bead and the first nucleotide of the following bead for each pair of
beads involved in the interaction calculation. However, in order to not
overestimate the bead-bead free-energy contribution, the interactions
involving nucleotides outside the considered beads are halved. As an example,
consider two sequences:
$5^{\prime}-(..A)(CGG)(C..)......((..G)(CCG)(U..)-3^{\prime}$ (see Figure 2).
The $\Delta H$ and $\Delta S$ for the interaction between the two middle beads
will have the following contributions:

$\displaystyle\begin{pmatrix}5^{\prime}-A)(CGG)(C-3^{\prime}\\\\[3.0pt]
3^{\prime}-U)(GCC)(G-5^{\prime}\end{pmatrix}=$
$\displaystyle\frac{1}{2}\begin{pmatrix}AC\\\\[3.0pt]
UG\end{pmatrix}+\begin{pmatrix}CG\\\\[3.0pt]
GC\end{pmatrix}+\begin{pmatrix}GG\\\\[3.0pt]
CC\end{pmatrix}+\frac{1}{2}\begin{pmatrix}GC\\\\[3.0pt] CG\end{pmatrix}.$ (3)

In order to establish a meaningful correspondence between the potential used
in our simulations and the NN parameters, the attraction strength between two
beads $i$ and $j$, $\epsilon_{ij}$, is computed by equaling the free-energy
change due to a bond with that predicted by the NN model, $\Delta H+T\Delta S$
Reinhardt and Frenkel (2016), which for the form and parameters of the bead-
bead attraction used here yields

$\beta\varepsilon_{ij}=-\frac{1}{\alpha}\left[\left(\mu\frac{\Delta
H}{k_{B}T_{37}}-\frac{\Delta
S}{k_{B}}\right)+\ln{(\rho^{\text{\sout{$\circ$}}}\sigma^{3}\tilde{V}_{b})}\right],$
(4)

where
$\rho^{\text{\sout{$\circ$}}}=[\text{\sout{$\circ$}}]\mathcal{N}_{A}=6.022\times
10^{26}\,\textrm{m}^{-3}$ is the standard number density, $\sigma$ is the
diameter of the beads and corresponds to 0.79 nm (see Figure 1b), while
$\tilde{V}_{b}=0.0019$ and $\alpha=0.89$ are computed via Monte Carlo
integration (see Appendix B for additional details).

We note that $\epsilon_{ij}$ is the only parameter affected by the type and
nature of the nucleotides: the geometry and structure of the strands are
otherwise independent of the nucleic acid modeled (DNA or RNA) and bead
sequence.

###  II.2 Bead design

The loss of resolution caused by dividing a strand into beads of $n$
nucleotides leads to a problem where native motifs (e.g. hairpins in natural
RNAs or crossovers in origami structures) may contain nucleotides from
multiple beads or a motif may occur in the middle of a bead. We do not want,
however, to make each bead correspond to a unique interaction, thereby losing
the ability to simulate unintended folds caused by competing complementary
sequences. To balance these competing demands, our bead design is based on the
native contacts observed in folded or designed structures, with the aim of
having the majority of the beads representing the same number of nucleotides,
$n$. The optimal number of nucleotides per bead is determined by the desired
level of detail and the specific system at hand, and it is set before running
the simulation. However, we allow for a small fraction of the remaining beads
to deviate by one nucleotide, either more or less, from the standard length,
$n$.

As a consequence, the beads are composed of varying numbers of nucleotides,
requiring the establishment of a rule for calculating interactions.
Specifically, we use the strongest $\Delta G$ found by pairing the shorter
bead with all possible subsequences of the same length contained within the
longer bead. In the following, we choose $n=3$ as the standard size of the
beads.

###  II.3 Parameters tuning

The free parameters of the model are the strength, denoted as $k_{s}$, and
width parameter, $\xi$, of the $V_{\rm semiflex}$ potential, as well as the
stacking strength, $\eta$, of the $V_{\rm stacking}$ potential. We optimize
the values of these quantities by comparing to the melting temperatures of DNA
hairpins predicted by the empirically-parameterized oxDNA model and to the
persistence lengths of single- and double-stranded DNA. For the thermodynamic
data, we focus on two sets of hairpins: twelve hairpins with a stem length of
six and a loop length of six, and twelve hairpins with a stem length of six
and a loop length of nine. The melting temperature is defined as the
temperature at which the yield of the closed hairpin is $0.5$, and the hairpin
is considered closed when at least 3 nucleotides are bonded for oxDNA
simulations), or one bead-bead bond is present, for the nNxB model. The
estimation of the persistence lengths is carried out by simulating 200-beads
hairpins with stems of 98 beads, and 100-beads single strands. In both cases,
we define the unit vector $\hat{v}_{i}$ connecting each pair $i$ and $i+1$,
and define $\cos(\theta_{ij})=\hat{v}_{i}\cdot\hat{v}_{j}$, where $i$ and $j$
are separated by $m$ beads. Averaging over all pairs (excluding the terminal
10 beads), we obtain an angular correlation,
$\langle\cos(\theta_{ij})\rangle$, which is a function of $m$ only. For worm-
like chains, such a correlation decays as

$\langle\cos(\theta_{ij})\rangle=e^{-m/l_{p}}.$ (5)

We apply Eq. (5) to extract the persistence length from our simulations. The
potential we use tends to favor the presence of kinks in double strands, which
decreases their persistence length. Therefore, we perform the analysis
described above with and without kinks, which are defined as $i,j$ neighboring
pairs for which $\cos(\theta_{ij})<0.7$. While the proper persistence length
requires that all angles are taken into account, excluding kinks provides an
estimate of the rigidity of the double-stranded parts of the structures that
do not bend much, which is the most common state of double strands in complex
structures (i.e. origami).

In our approach, we divide the DNA sequence into groups of three nucleotides
to form beads. Through extensive analysis, we found that the closest results
to the oxDNA predictions and the estimated persistence length for single- and
double-stranded DNA are obtained when using $\beta k_{s}=4$, $\xi=0.07$ and
$\beta\eta=6$. The resulting melting temperatures are shown in Table 1,
comparing very well to those extracted from oxDNA simulations. Indeed, the two
data sets differ by an average of just about two kelvins.

Figure 3 illustrates that the persistence length of the double strand is
approximately 28 beads with kinks and 49 without kinks, equivalent to around
84 and 147 base pairs (bp). In contrast, the derived persistence length of
single strands is approximately 3 beads, corresponding to 9 bp. These values
align well with estimates for real DNA, where the persistence length of
double-stranded DNA is approximately 50 nmWang  _et al._ (1997) (equivalent to
around 147 bp), and that of single-stranded DNA is approximately $(1.98\pm
0.72)$ nm Roth  _et al._ (2018) (equivalent to approximately $(5.8\pm 2.1)$
bp), considering a base pair spacing of 0.34 nm.

Since our model accommodates representing a variable number of nucleotides per
bead, to achieve a consistent persistence length in base pairs (bp),
regardless of the number of nucleotides represented by the bead, it is
possible to keep the parameters $\xi$ and $\eta$ constant while only adjusting
the parameter $k_{s}$ to achieve the desired persistence length.

6stem6loop | oxDNA | nNxB | $\Delta$
---|---|---|---
DNAh1 | 66 | 67 | 1
DNAh2 | 72 | 70 | -2
DNAh3 | 64 | 64 | 0
DNAh4 | 61 | 59 | -2
DNAh5 | 64 | 67 | 3
DNAh6 | 64 | 62 | -2
DNAh7 | 68 | 64 | -4
DNAh8 | 71 | 69 | -2
DNAh9 | 63 | 65 | 2
DNAh10 | 60 | 59 | -1
DNAh11 | 65 | 64 | -1
DNAh12 | 65 | 68 | 3
|  |  | $\langle|\Delta|\rangle$
|  |  | 1.9

6stem9loop | oxDNA | nNxB | $\Delta$
---|---|---|---
DNAh1 | 62 | 65 | 3
DNAh2 | 66 | 67 | 1
DNAh3 | 60 | 61 | 1
DNAh4 | 57 | 56 | -1
DNAh5 | 60 | 64 | 4
DNAh6 | 59 | 60 | 1
DNAh7 | 61 | 62 | 1
DNAh8 | 65 | 66 | 1
DNAh9 | 59 | 62 | 3
DNAh10 | 56 | 56 | 0
DNAh11 | 59 | 62 | 3
DNAh12 | 60 | 65 | 5
|  |  | $\langle|\Delta|\rangle$
|  |  | 2.0

Table 1: Melting temperatures in °C. $|\Delta|$ is the absolute value of the
difference between the oxDNA and nNxB predictions, and
$\langle|\Delta|\rangle$ is its average. Figure 3: Angular correlation as a
function of the chemical distance $m$ for double-stranded and single-stranded
DNA molecules. Points are simulation data, lines are exponential fits
performed in the $m\in[0,5]$ interval. Note that, as for real DNA and RNA,
single strands do not behave as worm-like chains and therefore the angular
correlation decays only approximately as an exponential.

##  III Results

###  III.1 Melting curves of hairpins

(a)
---
(b)
Figure 4: Comparison of the melting curves predicted by oxDNA (points) and
nNxB (lines) for selected DNA hairpins with a stem length of six and a loop
length of six LABEL:sub@subfig:_loop6, and selected DNA hairpins with a stem
length of six and a loop length of nine LABEL:sub@subfig:_loop9. The resulting
melting temperatures are listed in Table 1.

We start by comparing the melting curves of DNA hairpins simulated with oxDNA
and nNxB. A hairpin is considered closed when at least half of the nucleotides
in the stem, or equivalently half of the beads, is bonded. A representative
set of these results, taken from the dataset used to tune the model parameters
(see Table 1), is shown in Fig. 4. For all the hairpins considered, the
agreement is almost quantitative for a large temperature range, with oxDNA
hairpins exhibiting a somewhat narrower transition, i.e. slightly steeper
melting curves. At high temperatures the nNxB seems to fall off more slowly
compared to oxDNA: we ascribe this behavior to the fact that the bonding
volume of bead-bead interactions is rather large, so that there is always a
sizeable chance that two beads are considered to be bonded for purely
geometrical reasons even for temperatures at which the attraction strength is
small compared to the thermal energy.

We use hairpin systems to also estimate the performance difference between
oxDNA and nNxB. Specifically, we run single-core molecular dynamics
simulations of hairpins at the melting temperature and count the number of
transitions between the two states (open and closed). Dividing the resulting
number by the simulation wall time yields a factor of $\approx 100$, which we
use to estimate the sampling speed-up of the new model compared to oxDNA.
While the precise figure should (and will) depend on the system at hand, it is
reasonable to expect the speed-up to be of the same order of magnitude,
independently on the specifc conditions.

###  III.2 Pseudoknot

A pseudoknot is a complex folding pattern that can occur in RNA molecules. It
involves the formation of additional base pairs between distant regions of the
RNA sequence, resulting in a knot-like structure. Pseudoknots are known to
have crucial functions in various biological processes such as regulation of
gene expression and ribosome function Peselis and Serganov (2014), and viral
replicationNeupane  _et al._ (2021). They have also been proposed to have a
notable influence on RNA folding pathways Kucharík  _et al._ (2016). The
thermodynamics of pseudoknot base-pairing remains poorly experimental
characterized due to their enormous diversity and complexity. Furthermore, the
fact that they are composed of non-nested base pairs means that they cannot be
computed using the efficient dynamic-programming approaches which dominate in
the RNA secondary structure prediction fieldEddy (2004). Therefore, exact and
heuristic coarse-graining methods have mainly focused on pseudoknot-free
secondary structures.

(a)  |  (b)
---|---
Figure 5: LABEL:sub@subfig:_smartdiv A schematic representation of the
secondary structure of the MMTV pseudoknot, showing the bead division used to
simulate the strand with nNxB. Nucleotides inside the same box belong to the
same bead, while nucleotides outside the box belong to the bead that follows.
LABEL:sub@subfig:_nNxBpse nNxB equilibrium yields for the MMTV pseudoknot.
Filled symbols represent the results of simulations without swapping, while
empty symbols represent those with swapping. The transition temperatures
determined experimentally THEIMER and GIEDROC (2000) and numerically, with
oxRNA Šulc  _et al._ (2014), are indicated by violet dashed lines and cyan
dash-dotted lines, respectively.

We test our model by exploring the melting curve of the well-known MMTV
pseudoknot THEIMER and GIEDROC (2000). We divide the sequence in beads as
depicted in Figure 5a and we compare our predictions with experimental
calorimetry measurements THEIMER and GIEDROC (2000) and oxRNA simulations Šulc
_et al._ (2014), in Figure 5b. The yields of the two pseudoknot-precursor
motifs (hairpin 1 and 2, highlighted in blue and red in Fig. 5a) display the
same qualitative trends observed in oxRNA (reported in Ref. Šulc  _et al._
(2014)), although the peaks for hairpin 1 and hairpin 2 are higher and lower
in nNxB compared to oxRNA, respectively. Summing up the yields of the two
intermediate hairpins we obtain a curve (gold points in the figure) whose
intersection with the single strand and pseudoknot data identifies the two
transition temperatures, 83 ∘ C and 92 ∘C. Compared to the experimental
temperatures (74 ∘C and 95 ∘C), these values are at least as good as the ones
obtained with oxRNA simulations (68 ∘C and 85 ∘C) which, given the lower
resolution and much higher computational efficiency of the nNxB compare to the
latter, is an excellent result.

We have also run simulations without the swap mechanism, which when enabled
serves the purpose of speeding up equilibration and sampling, by setting
$\lambda=10$. As shown in Figure 5b, the results obtained with and without the
swap (empty and full symbols, respectively) overlap perfectly, demonstrating
that the thermodynamics is not affected by the value of $\lambda$.

###  III.3 Tile

Drawing inspiration from the idea of molecular tiles, which are individual
units that self-assemble to create intricate patterns or structuresWinfree
_et al._ (1998), we simulated the assembly of a single-stranded RNA tile
containing a k-type pseudoknotPoppleton  _et al._ (2020).

Unfortunately, no comparison with oxRNA is possible, as observing the folding
of such a structure, composed of 132 nucleotides, is currently out of reach
from the computational point of view. As before, we split the sequence into
beads of average size 3, but optimize the division by hand so that nucleotides
that are supposed to be paired in the native structure belong to the same
beads.

---

(a)  (b)

(c)

Figure 6: LABEL:sub@subfig:_Snap The RNA tile composed of 132 nucleotides
Poppleton  _et al._ (2020) we simulate, as represented by oxRNA.
LABEL:sub@subfig:_LastConf A nNxB configuration simulated at $28^{\circ}$C.
The highlighted (green) patches form kissing loops.LABEL:sub@subfig:_tileMelt
Melting curve of the RNA tile simulated with the nNxB model.

We simulate a single tile for different values of $\mu$ (which corresponds to
different temperatures), evaluating the fraction of bonds between any two
beads, and between beads that should be bonded in the native structure. The
results are shown in Fig. 6 and should be compared with the dashed horizontal
line that signals the fraction of bonds that are bonded in the native
structure ($\approx 0.7$). The nNxB model predicts a melting temperature (i.e.
a temperature at which half of the native contacts are formed) around 65∘, and
that below $\approx 50^{\circ}$ nearly all bonds are formed.

Figure 7: Circular representation of the RNA tile structure and of its
sequence, alongside the adopted division into beads and the bonds that connect
them at low temperature. Blue lines indicate native bonds, while red lines
indicate the most probable bonds that connect beads that are not fully
complementary but are nevertheless observed in simulations.

At low temperatures, on top of the native contacts, which are always formed,
we observe additional bonds that link beads that are not fully compatible with
each other. Figure 7 shows the tile sequence, the splitting into beads and the
native contacts (blue lines), as well as the most probable misbonds observed
(red lines). It is clear that most misbonds happen between beads that are part
of hairpin loops, thereby forming so-called kissing loops (highlighted in
figure 6b). Kissing loops complexes are known to play an important role in
RNA-RNA interactions, both in the biological and nanotechnology contexts
Paillart  _et al._ (1996); Bindewald  _et al._ (2008); Liu  _et al._ (2020).
Although is hard to estimate the stability of these motifs within the tile
with oxRNA or similar coarse-grained models, the lack of flexibility of the
tile arms should disfavour these particular kissing loops. Therefore, it is
possible that the nNxB model overestimates their stability.

##  IV Discussion and conclusions

This study presents a novel coarse-grained model aimed at simulating folding
processes of DNA and RNA nanostructures, nNxB. By representing $n$ nucleotides
with a single patchy particle, we have achieved a balance between
computational efficiency and the ability to capture the interactions that
govern folding dynamics, as demonstrated by simulations of DNA hairpins, an
RNA pseudoknot and an RNA tile. The model is parametrized using well-
established nearest-neighbor models, and can offer insights into the stability
and thermodynamics of nucleic acid structures, while opening avenues for the
exploration of larger systems and longer timescales. We showed that the
thermodynamic performance of the model is comparable with those of the best
nucleotide-level coarse-grained models (oxDNA and oxRNA), while being two
orders of magnitude faster. While the model focuses on thermodynamics, some
structural properties of nucleic acids, such as the different persistence
length between single- and double-stranded molecules, are also retained.
Simulations of the RNA tile, whose melting curve can be straightforwardly
obtained via unbiased simulations, show that the stability of some specific
motifs, such as kissing loops, may be overestimated.

Possible future applications of the model are the exploration of folding
pathways, optimization of single-stranded motifs and origami designs Han  _et
al._ (2017); Qi  _et al._ (2018), vaccine design Leppek  _et al._ (2022); Kim
_et al._ (2023), and viral RNA folding and packaging Božič  _et al._ (2018);
Vaupotič  _et al._ (2023). The current version of the model, presented here,
has been parametrized using single-stranded motifs only, but we plan to extend
it to multi-stranded systems in the future to make it possible to simulate,
for instance, the folding and optimization of DNA origami Gambietz, Stenke,
and Saccà (2023); Cumberworth, Frenkel, and Reinhardt (2022).

###### Acknowledgements.

We thank Francesco Sciortino for the fruitful discussions. PŠ acknowledges
support by the National Science Foundation under Grant DMR-2239518.

##  Appendix A The functional forms of the interaction potential

The potentials between topologically bonded beads (beads linked through the
backbone) are:

  * •

The Kremer-Grest force field Grest and Kremer (1986) $V_{KG}$, which is the
sum of a WCA potential and a FENE potential (see Figure 8). This spring-like
potential guarantees excluded volume (WCA component) and attraction (FENE
component) that mimic the covalent bonds along the strand. In particular,
defining $r$ as the distance between bonded beads, the WCA potential is

$V_{WCA}(r)=\begin{cases}4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]+\varepsilon,&\,r\leq
2^{1/6}\sigma\\\ 0,&\,r>2^{1/6}\sigma\end{cases}$ (6)

where $\varepsilon$ and $\sigma$ are the units of energy and length
respectively (both set to 1 in our simulations), while the FENE potential is

$V_{FENE}(r)=-\frac{1}{2}Kd_{0}^{2}\ln\left[1-\left(\frac{r}{d_{0}}\right)^{2}\right],$
(7)

where $d_{0}=1.5\sigma$ and $K=30\varepsilon/\sigma^{2}$.

  * •

A three-body potential that tends to align three consecutive beads (see Figure
9):

$V_{\rm
semiflex}(\theta)=-k\left[e^{-\left(\frac{1-\cos{\theta}}{\xi}\right)^{2}}-1\right].$
(8)

where $\theta$ is the angle defined by a triplet of bonded particles.

  * •

A term $V_{\rm stack}$ that models the stacking of base in DNA/RNA, acting on
the $\vec{a}_{1}$ directions of consecutive beads:

$V_{\rm
stack}(\vec{a}_{1}^{(i)},\vec{a}_{1}^{(i+1)})=\eta(1-\vec{a}_{1}^{(i)}\cdot\vec{a}_{1}^{(i+1)}).$
(9)

Figure 8: The different potentials acting between bonded beads. Figure 9: Plot
depicting a potential designed to simulate the distinct behaviors of DNA/RNA
in single- and double-stranded configurations (black curve). A widely used
cosine-angle potential to introduce structural rigidity is included for
comparison (red curve). The advantage of our potential is that when there are
fewer constraints like in the single-stranded case, the polymer is more likely
to rotate freely.

Non-bonded beads interact through the following potentials:

  * •

The same WCA potential, Eq. (6), is used to model the excluded-volume
interaction;

  * •

An attractive potential that models the hybridization of nucleotides, $V_{\rm
sticky}\cdot V_{\rm direct}$. Each patch has a position that is given by the
position of the bead it is attached to plus $\delta_{pb}\vec{a}_{1}$, where
$\delta_{pb}=0.65\sigma$, and $V_{\rm sticky}(r_{pp})$ acts between any two
patches, Its functional form was proposed by Stillinger and Weber Weber and
Stillinger (1985) (see Figure 10) and reads

$V_{\rm
sticky}^{ij}(r_{pp})=A\varepsilon_{ij}\left[B\left(\frac{\sigma_{s}}{r_{pp}}\right)^{4}-1\right]e^{\sigma_{s}/(r_{pp}-r_{c})}$
(10)

where $i$ and $j$ are the bead types and $r_{pp}$ is the distance between the
patches and

$\begin{split}B&=\frac{1}{1+4\left(1-r_{s}\right)^{2}},\\\
A&=-\frac{1}{B-1}\frac{1}{e^{1/1-r_{s}}},\\\
r_{c}&=\sigma_{s}r_{s}.\end{split}$ (11)

The coefficient $\varepsilon_{ij}$ modulates the strength of the sticky
attraction and depends on what nucleotides are inside the considered beads
(see Section II.1). We set $\varepsilon_{ij}\neq 0$ only if the chemical
distance between $i$ and $j$ is larger than two, since a loop of just 6
nucleotides does not occur in DNA/RNA strands. This potential depends
essentially on two parameters: $\sigma_{s}$, which defines the minimum of the
potential (i.e. the radius of the patch) and is set to $0.21875\sigma$, and
$r_{s}$, which defines the steepness of the potential between the minimum and
$r_{c}$ (set to $0.35$), after which $V_{\rm sticky}=0$. Thus, two particles
are bonded if the relative distance between their patches is less than
$r_{c}$. The directionality of DNA/RNA is enforced by multiplying $V_{\rm
sticky}$ by a term acting on the $\vec{a}_{3}$ directions of beads:

$V_{\rm
direct}(\vec{a}_{3}^{(i)},\vec{a}_{3}^{(j)})=\frac{\vec{a}_{3}^{(i)}\cdot\vec{a}_{3}^{(j)}}{2}.$
(12)

This term makes sure that only antiparallel strands can bind to each other.

  * •

$V_{\rm 3b}$, To ensure that each bead interacts only with one other, we
implement a repulsive three-body interaction $V_{\rm 3b}$ which penalizes the
formation of triplets of bonded beads Sciortino (2017). In particular, $V_{\rm
3b}$ is designed to almost exactly compensate the gain associated with the
formation of a second bond, originating an almost flat energy hypersurface
that favors bond swapping even in the presence of pair attraction energies
much larger than the thermal energy. This repulsive potential is defined as

$V_{\rm
3b}=\lambda\sum_{ijk}\min(\varepsilon_{ij},\varepsilon_{ik})V_{3}(r_{ij})V_{3}(r_{ik})$
(13)

where the sum runs over all triplets of bonded particles (bead $i$ bonded both
with $k$ and $j$), whereas we took the minimum between $\varepsilon_{ij}$ and
$\varepsilon_{ik}$ to favor the removal of the more loosely bonded bead.
$r_{ij}$ is the distance between bead $i$ and $j$. The value of the parameter
$\lambda$ allows to interpolate between the limits of swapping ($\lambda=1$)
and non-swapping ($\lambda\gg 1$) bonds. The pair potential $V_{3}(r)$ is
defined in terms of the normalized $V_{\rm sticky}(r)$ as

$V_{3}(r)=\begin{cases}1,&\qquad r\leq\sigma_{s}\\\ -\frac{V_{\rm
sticky}^{ij}(r)}{\varepsilon_{ij}},&\qquad\sigma_{s}\leq r\leq
r_{c}\end{cases}$ (14)

where $\sigma_{s}$ is the distance at which $V_{sticky}^{ij}(r)$ has a
minimum.

Figure 10: The different potentials acting between non-bonded beads.

##  Appendix B Mapping

To find a suitable mapping between the potential used in our simulations and
experimental data we followed the procedure of Ref.Reinhardt and Frenkel
(2016): we consider the hybridization of two single-stranded portions A and B
to give a hybridized (double-stranded) molecule AB:

$A+B\rightleftharpoons AB$ (15)

for which the equilibrium constant can be written as

$K=\frac{[AB]/[\text{\sout{$\circ$}}]}{([A]/[\text{\sout{$\circ$}}])([B]/[\text{\sout{$\circ$}}])}=\frac{\rho_{AB}\rho^{\text{\sout{$\circ$}}}}{\rho_{A}\rho_{B}}=\exp{(-\beta\Delta
G^{\text{\sout{$\circ$}}})},$ (16)

where $[\text{\sout{$\circ$}}]=1\,\textrm{mol dm}^{-3}$ is the standard state
concentration,
$\rho^{\text{\sout{$\circ$}}}=[\text{\sout{$\circ$}}]\mathcal{N}_{A}=6.022\times
10^{26}\,\textrm{m}^{-3}$ is the standard number density and $\Delta
G^{\text{\sout{$\circ$}}}$ is the standard Gibbs energy for the transformation
given in eq.15 where 50% of the monomers have hybridized. In our model, this
hybridization free energy is obtained from a nearest-neighbor thermodynamic
model. Assuming that the solution of monomers and dimers is ideal, which is
reasonable provided the concentration of each species is small, we can write
the canonical partition function of each species $x$ (where $x$ can be A, B,
or AB) as:

$Q_{x}=\frac{V^{N_{x}}}{N_{x}!\Lambda_{x}^{3N_{x}}}q_{x}^{N_{x}}$ (17)

where $V$ is the volume of the container, $N_{x}$ is the number of particles
of species $x$, $\Lambda_{x}$ is the de Broglie thermal wavelength of species
$x$ and $q_{x}$ is the internal partition function of species $x$. Each of the
chemical potentials can straightforwardly be calculated from the canonical
partition function:

$\mu_{x}=-k_{B}T\frac{\partial\ln{Q_{x}}}{\partial
N_{x}}=k_{B}T\ln{(\rho_{x}\Lambda_{x}^{3}/q_{x})}.$ (18)

At the equilibrium, $\mu_{A}+\mu_{B}=\mu_{AB}$, which we can solve as:

$\frac{\rho_{AB}}{\rho_{A}\rho_{B}}=\frac{q_{AB}\Lambda_{A}^{3}\Lambda_{B}^{3}}{q_{A}q_{B}\Lambda_{AB}^{3}}$
(19)

We assume that the internal state of the monomeric units that bind is not
affected by binding. We express this by setting the internal partition
functions of the two monomers equal to unity, $q_{A}=q_{B}=1$. Moreover, since
the AB molecule is described classically as a dimer of the A and B particles,
the de Broglie thermal wavelengths cancel out, since the momenta of the two
monomeric units in the dimer are uncoupled, The rotational partition function
of the dimer is thus subsumed into the translational degrees of freedom of the
constituting monomers, given that we integrate over the potential energy over
all possible states; we include this contribution in the internal partition
function $q_{AB}$, which will therefore have dimensions of volume. Hence, we
have:

$\frac{\rho_{AB}}{\rho_{A}\rho_{B}}=q_{AB}.$ (20)

and by comparing this equation with eq.16 we can write:

$q_{AB}\rho^{\text{\sout{$\circ$}}}=\exp{(-\beta\Delta
G^{\text{\sout{$\circ$}}})}.$ (21)

Knowing that:

$q_{AB}=4\pi\int r^{2}dr\int d\\{\vec{a}\\}\int d\\{\vec{b}\\}e^{-\beta V}$
(22)

where $V$ is the total interaction potential between the two beads.

We evaluate the right-hand side of Eq. 22 by means of a Monte Carlo
integration, finding that, for $\varepsilon\gtrsim 3$, $q_{AB}\simeq
V_{b}e^{\alpha\beta\varepsilon}$, where $V_{b}=0.0019\,\sigma^{3}$ and
$\alpha=0.89$ are fitting parameters. Substituting this relation in eq.21 we
find

$\tilde{V}_{b}e^{\alpha\beta\varepsilon_{ij}}\sigma^{3}\rho^{\text{\sout{$\circ$}}}=e^{-\beta\Delta
G^{\text{\sout{$\circ$}}}},$ (23)

so that

$\displaystyle\ln{(\rho^{\text{\sout{$\circ$}}}\sigma^{3}\tilde{V}_{b})}+\alpha\beta\varepsilon_{ij}$
$\displaystyle=-\beta\Delta G^{\text{\sout{$\circ$}}}$
$\displaystyle\beta\varepsilon_{ij}$
$\displaystyle=-\frac{1}{\alpha}\left[\beta\Delta
G^{\text{\sout{$\circ$}}}+\ln{(\rho^{\text{\sout{$\circ$}}}\sigma^{3}\tilde{V}_{b})}\right]$
(24)

Since $\Delta G=\Delta H-T\Delta S$, we can rewrite the previous equation as:

$\beta\varepsilon_{ij}=-\frac{1}{\alpha}\left[\frac{1}{k_{B}T}(\Delta
H-T\Delta
S)+\ln{(\rho^{\text{\sout{$\circ$}}}\sigma^{3}\tilde{V}_{b})}\right],$ (25)

Using eq.II.1, we obtain:

$\beta\varepsilon_{ij}=-\frac{1}{\alpha}\left[\left(\mu\frac{\Delta
H}{k_{B}T_{37}}-\frac{\Delta
S}{k_{B}}\right)+\ln{(\rho^{\text{\sout{$\circ$}}}\sigma^{3}\tilde{V}_{b})}\right].$
(26)

## References

  * Pitchiaya and Krishnan (2006) S. Pitchiaya and Y. Krishnan, “First blueprint, now bricks: Dna as construction material on the nanoscale,” Chemical Society Reviews 35, 1111–1121 (2006).
  * Dey  _et al._ (2021) S. Dey, C. Fan, K. V. Gothelf, J. Li, C. Lin, L. Liu, N. Liu, M. A. Nijenhuis, B. Saccà, F. C. Simmel, _et al._ , “Dna origami,” Nature Reviews Methods Primers 1, 13 (2021).
  * Tinoco Jr and Bustamante (1999) I. Tinoco Jr and C. Bustamante, “How rna folds,” Journal of molecular biology 293, 271–281 (1999).
  * Tang  _et al._ (2021) Z. Tang, X. Zhang, Y. Shu, M. Guo, H. Zhang, and W. Tao, “Insights from nanotechnology in covid-19 treatment,” Nano today 36, 101019 (2021).
  * Li  _et al._ (2018) S. Li, Q. Jiang, S. Liu, Y. Zhang, Y. Tian, C. Song, J. Wang, Y. Zou, G. J. Anderson, J.-Y. Han, _et al._ , “A dna nanorobot functions as a cancer therapeutic in response to a molecular trigger in vivo,” Nature biotechnology 36, 258–264 (2018).
  * Benenson (2009) Y. Benenson, “Rna-based computation in live cells,” Current opinion in biotechnology 20, 471–478 (2009).
  * Du, Hu, and Ren (2023) C. Du, P. Hu, and L. Ren, “Nucleic acid-based scaffold systems and application in enzyme cascade catalysis,” Applied Microbiology and Biotechnology 107, 9–23 (2023).
  * Seeman (2010) N. C. Seeman, “Nanomaterials based on dna,” Annual review of biochemistry 79, 65–87 (2010).
  * Rothemund (2006) P. W. Rothemund, “Folding dna to create nanoscale shapes and patterns,” Nature 440, 297–302 (2006).
  * Douglas  _et al._ (2009) S. M. Douglas, H. Dietz, T. Liedl, B. Högberg, F. Graf, and W. M. Shih, “Self-assembly of dna into nanoscale three-dimensional shapes,” Nature 459, 414–418 (2009).
  * Bath  _et al._ (2009) J. Bath, S. J. Green, K. E. Allen, and A. J. Turberfield, “Mechanism for a directional, processive, and reversible dna motor,” Small 5, 1513–1516 (2009).
  * Ouldridge  _et al._ (2013a) T. E. Ouldridge, R. L. Hoare, A. A. Louis, J. P. Doye, J. Bath, and A. J. Turberfield, “Optimizing dna nanotechnology through coarse-grained modeling: a two-footed dna walker,” ACS nano 7, 2479–2490 (2013a).
  * Shi  _et al._ (2023) X. Shi, A.-K. Pumm, C. Maffeo, F. Kohler, E. Feigl, W. Zhao, D. Verschueren, R. Golestanian, A. Aksimentiev, H. Dietz, _et al._ , “A dna turbine powered by a transmembrane potential across a nanopore,” Nature Nanotechnology , 1–7 (2023).
  * Adleman (1994) L. M. Adleman, “Molecular computation of solutions to combinatorial problems,” science 266, 1021–1024 (1994).
  * Chen  _et al._ (2013) Y.-J. Chen, N. Dalchau, N. Srinivas, A. Phillips, L. Cardelli, D. Soloveichik, and G. Seelig, “Programmable chemical controllers made from dna,” Nature nanotechnology 8, 755–762 (2013).
  * Wang  _et al._ (2023) B. Wang, S. S. Wang, C. Chalk, A. D. Ellington, and D. Soloveichik, “Parallel molecular computation on digital data stored in dna,” Proceedings of the National Academy of Sciences 120, e2217330120 (2023).
  * Cornell  _et al._ (1995) W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, and P. A. Kollman, “A second generation force field for the simulation of proteins, nucleic acids, and organic molecules,” Journal of the American Chemical Society 117, 5179–5197 (1995).
  * Brooks  _et al._ (1983) B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. a. Swaminathan, and M. Karplus, “Charmm: a program for macromolecular energy, minimization, and dynamics calculations,” Journal of computational chemistry 4, 187–217 (1983).
  * Sponer  _et al._ (2018) J. Sponer, G. Bussi, M. Krepl, P. Banáš, S. Bottaro, R. A. Cunha, A. Gil-Ley, G. Pinamonti, S. Poblete, P. Jurecka, _et al._ , “Rna structural dynamics as captured by molecular simulations: a comprehensive overview,” Chemical reviews 118, 4177–4338 (2018).
  * Parisien and Major (2008) M. Parisien and F. Major, “The mc-fold and mc-sym pipeline infers rna structure from sequence data,” Nature 452, 51–55 (2008).
  * Savelyev and Papoian (2009) A. Savelyev and G. A. Papoian, “Molecular renormalization group coarse-graining of polymer chains: application to double-stranded dna,” Biophysical Journal 96, 4044–4052 (2009).
  * Paliy, Melnik, and Shapiro (2010) M. Paliy, R. Melnik, and B. A. Shapiro, “Coarse-graining rna nanostructures for molecular dynamics simulations,” Physical biology 7, 036001 (2010).
  * Pasquali and Derreumaux (2010) S. Pasquali and P. Derreumaux, “Hire-rna: a high resolution coarse-grained energy model for rna,” The journal of physical chemistry B 114, 11957–11966 (2010).
  * Ouldridge, Louis, and Doye (2011) T. E. Ouldridge, A. A. Louis, and J. P. Doye, “Structural, mechanical, and thermodynamic properties of a coarse-grained dna model,” The Journal of chemical physics 134, 02B627 (2011).
  * Šulc  _et al._ (2012) P. Šulc, F. Romano, T. E. Ouldridge, L. Rovigatti, J. P. Doye, and A. A. Louis, “Sequence-dependent thermodynamics of a coarse-grained dna model,” The Journal of chemical physics 137, 135101 (2012).
  * Hinckley  _et al._ (2013) D. M. Hinckley, G. S. Freeman, J. K. Whitmer, and J. J. De Pablo, “An experimentally-informed coarse-grained 3-site-per-nucleotide model of dna: Structure, thermodynamics, and dynamics of hybridization,” The Journal of chemical physics 139, 10B604_1 (2013).
  * Cragnolini, Derreumaux, and Pasquali (2013) T. Cragnolini, P. Derreumaux, and S. Pasquali, “Coarse-grained simulations of rna and dna duplexes,” The Journal of Physical Chemistry B 117, 8047–8060 (2013).
  * Xia  _et al._ (2013) Z. Xia, D. R. Bell, Y. Shi, and P. Ren, “Rna 3d structure prediction by using a coarse-grained model and experimental data,” The Journal of Physical Chemistry B 117, 3135–3144 (2013).
  * Denesyuk and Thirumalai (2013) N. A. Denesyuk and D. Thirumalai, “Coarse-grained model for predicting rna folding thermodynamics,” The Journal of Physical Chemistry B 117, 4901–4911 (2013).
  * Šulc  _et al._ (2014) P. Šulc, F. Romano, T. E. Ouldridge, J. P. Doye, and A. A. Louis, “A nucleotide-level coarse-grained model of rna,” The Journal of chemical physics 140, 06B614_1 (2014).
  * Korolev  _et al._ (2014) N. Korolev, D. Luo, A. P. Lyubartsev, and L. Nordenskiöld, “A coarse-grained dna model parameterized from atomistic simulations by inverse monte carlo,” Polymers 6, 1655–1675 (2014).
  * Maciejczyk  _et al._ (2014) M. Maciejczyk, A. Spasic, A. Liwo, and H. A. Scheraga, “Dna duplex formation with a coarse-grained model,” Journal of chemical theory and computation 10, 5020–5035 (2014).
  * Maffeo  _et al._ (2014) C. Maffeo, T. T. Ngo, T. Ha, and A. Aksimentiev, “A coarse-grained model of unstructured single-stranded dna derived from atomistic simulation and single-molecule experiment,” Journal of chemical theory and computation 10, 2891–2896 (2014).
  * Rovigatti  _et al._ (2015) L. Rovigatti, P. Šulc, I. Z. Reguly, and F. Romano, “A comparison between parallelization approaches in molecular dynamics simulations on gpus,” Journal of computational chemistry 36, 1–8 (2015).
  * Machado and Pantano (2015) M. R. Machado and S. Pantano, “Exploring laci–dna dynamics by multiscale simulations using the sirah force field,” Journal of chemical theory and computation 11, 5012–5023 (2015).
  * Uusitalo  _et al._ (2015) J. J. Uusitalo, H. I. Ingoólfsson, P. Akhshi, D. P. Tieleman, and S. J. Marrink, “Martini coarse-grained force field: extension to dna,” Journal of chemical theory and computation 11, 3932–3945 (2015).
  * Dans  _et al._ (2016) P. D. Dans, J. Walther, H. Gómez, and M. Orozco, “Multiscale simulation of dna,” Current opinion in structural biology 37, 29–45 (2016).
  * Ivani  _et al._ (2016) I. Ivani, P. D. Dans, A. Noy, A. Pérez, I. Faustino, A. Hospital, J. Walther, P. Andrio, R. Goñi, A. Balaceanu, _et al._ , “Parmbsc1: a refined force field for dna simulations,” Nature methods 13, 55–58 (2016).
  * Chakraborty, Hori, and Thirumalai (2018) D. Chakraborty, N. Hori, and D. Thirumalai, “Sequence-dependent three interaction site model for single-and double-stranded dna,” Journal of chemical theory and computation 14, 3763–3779 (2018).
  * Maffeo and Aksimentiev (2020) C. Maffeo and A. Aksimentiev, “Mrdna: a multi-resolution model for predicting the structure and dynamics of dna systems,” Nucleic acids research 48, 5135–5146 (2020).
  * DeLuca  _et al._ (2023) M. DeLuca, T. Ye, M. Poirier, Y. Ke, C. Castro, and G. Arya, “Mechanism of dna origami folding elucidated by mesoscopic simulations,” bioRxiv (2023), 10.1101/2023.06.20.545758.
  * Sengar  _et al._ (2021) A. Sengar, T. E. Ouldridge, O. Henrich, L. Rovigatti, and P. Šulc, “A primer on the oxdna model of dna: when to use it, how to simulate it and how to interpret the results,” Frontiers in Molecular Biosciences 8, 693710 (2021).
  * Poland and Scheraga (1966) D. Poland and H. A. Scheraga, “Occurrence of a phase transition in nucleic acid models,” The Journal of chemical physics 45, 1464–1469 (1966).
  * SantaLucia Jr and Hicks (2004) J. SantaLucia Jr and D. Hicks, “The thermodynamics of dna structural motifs,” Annu. Rev. Biophys. Biomol. Struct. 33, 415–440 (2004).
  * Allawi and SantaLucia (1998a) H. T. Allawi and J. SantaLucia, “Nearest-neighbor thermodynamics of internal a. c mismatches in dna: Sequence dependence and ph effects,” Biochemistry 37, 9435–9444 (1998a).
  * Allawi and SantaLucia Jr (1998) H. T. Allawi and J. SantaLucia Jr, “Thermodynamics of internal c· t mismatches in dna,” Nucleic acids research 26, 2694–2701 (1998).
  * Allawi and SantaLucia (1998b) H. T. Allawi and J. SantaLucia, “Nearest neighbor thermodynamic parameters for internal g$\odot$ a mismatches in dna,” Biochemistry 37, 2170–2179 (1998b).
  * Allawi and SantaLucia (1997) H. T. Allawi and J. SantaLucia, “Thermodynamics and nmr of internal g$\odot$ t mismatches in dna,” Biochemistry 36, 10581–10594 (1997).
  * Peyret  _et al._ (1999) N. Peyret, P. A. Seneviratne, H. T. Allawi, and J. SantaLucia, “Nearest-neighbor thermodynamics and nmr of dna sequences with internal a$\odot$ a, c$\odot$ c, g$\odot$ g, and t$\odot$ t mismatches,” Biochemistry 38, 3468–3477 (1999).
  * Turner and Mathews (2010) D. H. Turner and D. H. Mathews, “Nndb: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure,” Nucleic acids research 38, D280–D282 (2010).
  * Huguet  _et al._ (2010) J. M. Huguet, C. V. Bizarro, N. Forns, S. B. Smith, C. Bustamante, and F. Ritort, “Single-molecule derivation of salt dependent base-pair free energies in dna,” Proceedings of the National Academy of Sciences 107, 15431–15436 (2010).
  * Bae, Fang, and Zhang (2020) J. H. Bae, J. Z. Fang, and D. Y. Zhang, “High-throughput methods for measuring dna thermodynamics,” Nucleic acids research 48, e89–e89 (2020).
  * Zuber  _et al._ (2022) J. Zuber, S. J. Schroeder, H. Sun, D. H. Turner, and D. H. Mathews, “Nearest neighbor rules for rna helix folding thermodynamics: improved end effects,” Nucleic Acids Research 50, 5251–5262 (2022).
  * Grest and Kremer (1986) G. S. Grest and K. Kremer, “Molecular dynamics simulation for polymers in the presence of a heat bath,” Physical Review A 33, 3628 (1986).
  * Weber and Stillinger (1985) T. A. Weber and F. H. Stillinger, “Local order and structural transitions in amorphous metal-metalloid alloys,” Physical Review B 31, 1954 (1985).
  * Sciortino (2017) F. Sciortino, “Three-body potential for simulating bond swaps in molecular dynamics,” The European Physical Journal E 40, 1–4 (2017).
  * Russo, Tartaglia, and Sciortino (2009) J. Russo, P. Tartaglia, and F. Sciortino, “Reversible gels of patchy particles: role of the valence,” The Journal of Chemical Physics 131, 014504 (2009).
  * Poppleton  _et al._ (2023) E. Poppleton, M. Matthies, D. Mandal, F. Romano, P. Šulc, and L. Rovigatti, “oxdna: coarse-grained simulations of nucleic acids made simple,” Journal of Open Source Software 8, 4693 (2023).
  * Reinhardt and Frenkel (2016) A. Reinhardt and D. Frenkel, “Dna brick self-assembly with an off-lattice potential,” Soft Matter 12, 6253–6260 (2016).
  * Wang  _et al._ (1997) M. D. Wang, H. Yin, R. Landick, J. Gelles, and S. M. Block, “Stretching dna with optical tweezers,” Biophysical journal 72, 1335–1346 (1997).
  * Roth  _et al._ (2018) E. Roth, A. Glick Azaria, O. Girshevitz, A. Bitler, and Y. Garini, “Measuring the conformation and persistence length of single-stranded dna using a dna origami structure,” Nano letters 18, 6703–6709 (2018).
  * Peselis and Serganov (2014) A. Peselis and A. Serganov, “Structure and function of pseudoknots involved in gene expression control,” Wiley Interdisciplinary Reviews: RNA 5, 803–822 (2014).
  * Neupane  _et al._ (2021) K. Neupane, M. Zhao, A. Lyons, S. Munshi, S. M. Ileperuma, D. B. Ritchie, N. Q. Hoffer, A. Narayan, and M. T. Woodside, “Structural dynamics of single sars-cov-2 pseudoknot molecules reveal topologically distinct conformers,” Nature Communications 12, 4749 (2021).
  * Kucharík  _et al._ (2016) M. Kucharík, I. L. Hofacker, P. F. Stadler, and J. Qin, “Pseudoknots in rna folding landscapes,” Bioinformatics 32, 187–194 (2016).
  * Eddy (2004) S. R. Eddy, “How do rna folding algorithms work?” Nature biotechnology 22, 1457–1458 (2004).
  * THEIMER and GIEDROC (2000) C. A. THEIMER and D. P. GIEDROC, “Contribution of the intercalated adenosine at the helical junction to the stability of the gag-pro frameshifting pseudoknot from mouse mammary tumor virus,” RNA 6, 409–421 (2000).
  * Winfree  _et al._ (1998) E. Winfree, F. Liu, L. A. Wenzler, and N. C. Seeman, “Design and self-assembly of two-dimensional dna crystals,” Nature 394, 539–544 (1998).
  * Poppleton  _et al._ (2020) E. Poppleton, J. Bohlin, M. Matthies, S. Sharma, F. Zhang, and P. Šulc, “Design, optimization and analysis of large dna and rna nanostructures through interactive visualization, editing and molecular simulation,” Nucleic acids research 48, e72–e72 (2020).
  * Paillart  _et al._ (1996) J.-C. Paillart, E. Skripkin, B. Ehresmann, C. Ehresmann, and R. Marquet, “A loop-loop" kissing" complex is the essential part of the dimer linkage of genomic hiv-1 rna.” Proceedings of the National Academy of Sciences 93, 5572–5577 (1996).
  * Bindewald  _et al._ (2008) E. Bindewald, R. Hayes, Y. G. Yingling, W. Kasprzak, and B. A. Shapiro, “Rnajunction: a database of rna junctions and kissing loops for three-dimensional structural analysis and nanodesign,” Nucleic acids research 36, D392–D397 (2008).
  * Liu  _et al._ (2020) D. Liu, C. W. Geary, G. Chen, Y. Shao, M. Li, C. Mao, E. S. Andersen, J. A. Piccirilli, P. W. Rothemund, and Y. Weizmann, “Branched kissing loops for the construction of diverse rna homooligomeric nanostructures,” Nature Chemistry 12, 249–259 (2020).
  * Han  _et al._ (2017) D. Han, X. Qi, C. Myhrvold, B. Wang, M. Dai, S. Jiang, M. Bates, Y. Liu, B. An, F. Zhang, _et al._ , “Single-stranded dna and rna origami,” Science 358, eaao2648 (2017).
  * Qi  _et al._ (2018) X. Qi, F. Zhang, Z. Su, S. Jiang, D. Han, B. Ding, Y. Liu, W. Chiu, P. Yin, and H. Yan, “Programming molecular topologies from single-stranded nucleic acids,” Nature communications 9, 4579 (2018).
  * Leppek  _et al._ (2022) K. Leppek, G. W. Byeon, W. Kladwang, H. K. Wayment-Steele, C. H. Kerr, A. F. Xu, D. S. Kim, V. V. Topkar, C. Choe, D. Rothschild, _et al._ , “Combinatorial optimization of mrna structure, stability, and translation for rna-based therapeutics,” Nature communications 13, 1536 (2022).
  * Kim  _et al._ (2023) Y.-A. Kim, K. Mousavi, A. Yazdi, M. Zwierzyna, M. Cardinali, D. Fox, T. Peel, J. Coller, K. Aggarwal, and G. Maruggi, “Computational design of mrna vaccines,” Vaccine (2023), https://doi.org/10.1016/j.vaccine.2023.07.024.
  * Božič  _et al._ (2018) A. L. Božič, C. Micheletti, R. Podgornik, and L. Tubiana, “Compactness of viral genomes: effect of disperse and localized random mutations,” Journal of Physics: Condensed Matter 30, 084006 (2018).
  * Vaupotič  _et al._ (2023) D. Vaupotič, A. Rosa, L. Tubiana, and A. Božič, “Scaling properties of RNA as a randomly branching polymer,” The Journal of Chemical Physics 158, 234901 (2023).
  * Gambietz, Stenke, and Saccà (2023) S. Gambietz, L. J. Stenke, and B. Saccà, “Sequence-dependent folding of monolayered dna origami domains,” Nanoscale 15, 13120–13132 (2023).
  * Cumberworth, Frenkel, and Reinhardt (2022) A. Cumberworth, D. Frenkel, and A. Reinhardt, “Simulations of dna-origami self-assembly reveal design-dependent nucleation barriers,” Nano Letters 22, 6916–6922 (2022), pMID: 36037484.
  * Veneziano  _et al._ (2016) R. Veneziano, S. Ratanalert, K. Zhang, F. Zhang, H. Yan, W. Chiu, and M. Bathe, “Designer nanoscale dna assemblies programmed from the top down,” Science 352, 1534–1534 (2016).
  * Mathews  _et al._ (2004) D. H. Mathews, M. D. Disney, J. L. Childs, S. J. Schroeder, M. Zuker, and D. H. Turner, “Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of rna secondary structure,” Proceedings of the National Academy of Sciences 101, 7287–7292 (2004).
  * Doye  _et al._ (2013) J. P. Doye, T. E. Ouldridge, A. A. Louis, F. Romano, P. Šulc, C. Matek, B. E. Snodin, L. Rovigatti, J. S. Schreck, R. M. Harrison, _et al._ , “Coarse-graining dna for simulations of dna nanotechnology,” Physical Chemistry Chemical Physics 15, 20395–20414 (2013).
  * Hyeon, Dima, and Thirumalai (2006) C. Hyeon, R. I. Dima, and D. Thirumalai, “Size, shape, and flexibility of rna structures,” The Journal of chemical physics 125, 194905 (2006).
  * Abels  _et al._ (2005) J. Abels, F. Moreno-Herrero, T. Van der Heijden, C. Dekker, and N. H. Dekker, “Single-molecule measurements of the persistence length of double-stranded rna,” Biophysical journal 88, 2737–2744 (2005).
  * Studnicka  _et al._ (1978) G. M. Studnicka, G. M. Rahn, I. W. Cummings, and W. A. Salser, “Computer method for predicting the secondary structure of single-stranded rna,” Nucleic Acids Research 5, 3365–3388 (1978).
  * Antczak  _et al._ (2018) M. Antczak, M. Popenda, T. Zok, M. Zurkowski, R. W. Adamiak, and M. Szachniuk, “New algorithms to represent complex pseudoknotted rna structures in dot-bracket notation,” Bioinformatics 34, 1304–1312 (2018).
  * Bon  _et al._ (2008) M. Bon, G. Vernizzi, H. Orland, and A. Zee, “Topological classification of rna structures,” Journal of molecular biology 379, 900–911 (2008).
  * Hochrein  _et al._ (2013) L. M. Hochrein, M. Schwarzkopf, M. Shahgholi, P. Yin, and N. A. Pierce, “Conditional dicer substrate formation via shape and sequence transduction with small conditional rnas,” Journal of the American Chemical Society 135, 17322–17330 (2013).
  * Ouldridge  _et al._ (2013b) T. E. Ouldridge, P. Šulc, F. Romano, J. P. Doye, and A. A. Louis, “Dna hybridization kinetics: zippering, internal displacement and sequence dependence,” Nucleic acids research 41, 8886–8895 (2013b).
  * Srinivas  _et al._ (2013) N. Srinivas, T. E. Ouldridge, P. Šulc, J. M. Schaeffer, B. Yurke, A. A. Louis, J. P. Doye, and E. Winfree, “On the biophysics and kinetics of toehold-mediated dna strand displacement,” Nucleic acids research 41, 10641–10658 (2013).
  * Irmisch, Ouldridge, and Seidel (2020) P. Irmisch, T. E. Ouldridge, and R. Seidel, “Modeling dna-strand displacement reactions in the presence of base-pair mismatches,” Journal of the American Chemical Society 142, 11451–11463 (2020).
  * Romano  _et al._ (2013) F. Romano, D. Chakraborty, J. P. Doye, T. E. Ouldridge, and A. A. Louis, “Coarse-grained simulations of dna overstretching,” The Journal of chemical physics 138, 02B610 (2013).
  * Wang and Pettitt (2014) Q. Wang and B. M. Pettitt, “Modeling dna thermodynamics under torsional stress,” Biophysical journal 106, 1182–1193 (2014).
  * Sutthibutpong  _et al._ (2016) T. Sutthibutpong, C. Matek, C. Benham, G. G. Slade, A. Noy, C. Laughton, J. P. K. Doye, A. A. Louis, and S. A. Harris, “Long-range correlations in the mechanics of small dna circles under topological stress revealed by multi-scale simulation,” Nucleic acids research 44, 9121–9130 (2016).
  * Rovigatti, Bomboi, and Sciortino (2014) L. Rovigatti, F. Bomboi, and F. Sciortino, “Accurate phase diagram of tetravalent dna nanostars,” The Journal of Chemical Physics 140, 154903 (2014).
  * Engel  _et al._ (2020) M. C. Engel, F. Romano, A. A. Louis, and J. P. Doye, “Measuring internal forces in single-stranded dna: Application to a dna force clamp,” Journal of chemical theory and computation 16, 7764–7775 (2020).