{'fraction_non_alphanumeric': 0.02433005026031689, 'fraction_numerical': 0.011077554148166512, 'mean_word_length': 4.726944175540497, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A nucleotide-level coarse-grained model of RNA

We present a new, nucleotide-level model for RNA, oxRNA, based on the coarse-graining methodology recently developed for the oxDNA model of DNA. The model is designed to reproduce structural, mechanical and thermodynamic properties of RNA, and the coarse-graining level aims to retain the relevant physics for RNA hybridization and the structure of single- and double-stranded RNA. In order to explore its strengths and weaknesses, we test the model in a range of nanotechnological and biological settings. Applications explored include the folding thermodynamics of a pseudoknot, the formation of a kissing loop complex, the structure of a hexagonal RNA nanoring, and the unzipping of a hairpin motif. We argue that the model can be used for efficient simulations of the structure of systems with thousands of base pairs, and for the assembly of systems of up to hundreds of base pairs. The source code implementing the model is released for public use.


I. INTRODUCTION
RNA (ribonucleic acid) strands are polymers that play crucial cellular roles in gene expression, translation and regulation. 1 RNA is composed of units called nucleotides, which consist of a phosphate and ribose sugar backbone to which one of four different bases is attached: adenine (A), guanine (G), cytosine (C) or uracil (U). DNA has a similar chemical structure, but thymine (T) is present in place of uracil and the backbone includes deoxyribose sugars instead of ribose sugars. Both RNA and DNA can form double-helical molecules, stabilized by hydrogen bonds between complementary Watson-Crick base pairs: AU (AT in case of DNA) and GC. For RNA, wobble base pairs (GU) can also stabilize the duplex form and are commonly found. While DNA typically forms a B-helix, double-stranded RNA is typically found in a related, but different A-helix geometry. 2 DNA's role in biological systems is to store genetic information and it is most often found in its doublehelical form. By contrast, most naturally occurring RNA molecules are single-stranded and fold into complex structures that contain double-helical segments as well as loops, junctions and numerous tertiary structure interactions that further stabilize the molecule. The lengths of RNA strands found in living cells can range from a few nucleotides to several thousands or even tens of thousands. 3 DNA nanotechnology 4 uses DNA molecules to create nanoscale structures and active devices. Examples of experimentally realized systems include DNA motors, 5,6 self-assembled nanostructures such as DNA origamis 7 and the use of DNA strands for computation. 8 Since RNA molecules are more difficult to handle and preserve than DNA in experimental settings, developing RNAbased nanodevices has been a more challenging task. However, the emerging field of RNA nanotechnology 9 offers promising applications of RNA-based nanodevices both in vitro and in vivo. The general design principles from DNA nanotechnology can be applied to RNA, but because specific RNA sequences form functional structures that are known to interact with proteins and other RNA molecules in the cell, RNA molecules would be the natural choice for nanotechnology devices operating inside cells. 10 Successful examples of experimentally realized RNA nanotechnology include self-assembling RNA nanocubes 11 and RNA tiles. 12 Recently, an RNA strand displacement reaction cascade was proposed that could be used for conditional gene silencing inside the cell. 13 Because of its importance for biological systems, RNA has been the subject of intensive study. In addition to numerous experimental studies, a range of theoretical and computational methods have been applied to the study of its properties. Currently, there are multiple tools and computational approaches that can describe RNA structure at various levels of detail. 14,15 In an important series of works [16][17][18][19][20][21][22] the thermodynamics of RNA secondary structure (i.e., the list of base pairs in the folded state of RNA) was characterized in terms of a nearest-neighbor model for calculating the free energies of various secondary-structure motifs. The nearestneighbor model is the basis of various tools for the prediction of the secondary structure. Such tools typically use dynamic programming approaches to find the secondary structure with minimal free energy. [23][24][25] Furthermore, some tools have been extended by adding simple kinetic descriptions to the nearest-neighbor thermodynamics, allowing folding transitions to be modeled. 26,27 Although these methods are typically very fast, the fundamentally discrete nature of the description and the lack of structural and mechanical detail places a limit on what they can treat.
At the highest level of detail, quantum chemistry methods have been employed to study the interactions between nucleotides. 28,29 Due to the complexity of the calculations, such methods remain limited to interactions between nearest-neighbor base pairs in vacuum. Fully atomistic molecular dynamics simulation packages such as Amber 30 or CHARMM 31 include the solvent, and use classical effective interaction potentials between atoms to represent systems with high resolution. However, the study of rare events such as the formation or breaking of individual base pairs remains a very challenging task. Simulations of several folding pathways of a short RNA hairpin 32 and tetraloops 33 provide examples of the limit of what is currently possible. At present, atomistic molecular dynamics simulations can attain timescales of the order of µs. Moreover, while the forcefields are improving, they are still under development so that different versions can generate different behavior. [34][35][36] A recently developed approach 37 combines fully atomistic representation with hierarchical Monte Carlo sampling, where different series of moves are used to move whole sections of a molecule (such as all atoms contained in one stem) at once. Such methods have for instance been used to study the effects of mutations in a sequence on the conformational freedom of a tRNA molecule and of a nanosquare composed of four tRNAs. 38 In order to access longer timescales relevant to rare events, such as the breaking of base pairs or the formation of large structures, one needs to use a more coarsegrained description. In this approach, atoms are incorporated into a reduced set of degrees of freedom that experience effective interactions. Solvent molecules are often integrated out. Such models always present a compromise between accuracy, efficiency and the level of detail, which determines their scope. Several coarse-grained RNA models have been developed in recent years. [39][40][41][42][43][44][45][46][47][48][49][50][51] Knowledge-based coarse-graining uses the information extracted from experimentally determined crystal structures to develop potentials, usually with the goal to predict the folded structure for an RNA sequence, either de novo or with some additional input of data from the user. An example of such an approach is the NAST model 39 which represents each nucleotide as a single pseudoatom in the simulation and uses a statistical potential, inferred from known structures of RNA molecules, that depends on the distances and angles between the nucleotides. This model requires the secondary structure and tertiary contacts of the final folded RNA structure as an input for the folding simulation. It has been used to study RNA structures of up to 158 nucleotides. While it was able to reproduce the structures of folded RNA, it was not parametrized to reproduce their thermodynamic properties.
Similarly, the model of Xia et al., 40 which uses 6 beads to represent a nucleotide, has interactions parametrized to reproduce known RNA structures. Xia et al. were able to predict the tertiary structure of several RNA molecules of lengths up to 122 nucleotides by using simulated annealing to attempt to find the global potential energy minimum of the structures in their coarse-grained model. The resulting structures were then refined by a simulation with a fully atomic representation. The thermodynamic properties of the coarse-grained model were not reported.
Finally, some knowledge-based methods combine together various structural motifs from database of experi-mentally determined RNA structures to predict folded RNA structure for a given sequence. The algorithms match fragments with a particular section in the RNA strand. For example, Parisien and Major 41 used such an approach to predict the secondary and tertiary structure of RNA strands with up to 50 nucleotides. The FARFAR method 42 further uses sampling with a fully atomistic representation of the respective RNA residues in order to obtain the final structure. It successfully predicted de novo folded structures for RNA sequences of size up to 20 nucleotides.
An alternative approach for model development is to use fully atomistic simulations to parametrize the effective interactions between coarse-grained representations of groups of atoms. Such an approach was adopted, for example, by Paliy et al., 43 who presented different levels of coarse-graining, using either one or three beads per nucleotide. The interactions between beads were fitted to reproduce the probability distribution of their mutual orientations and distances, calculated from from a simulation with a fully atomistic representation. The authors were then able to simulate the conformations of an RNA nanoring structure which consisted of 330 nucleotides.
The HiRe-RNA model 44,45 represents each nucleotide as 6 or 7 beads with empirically chosen interactions based on a combination of atomistic simulations and known structures. It reproduces the structure of RNA duplexes and was used to simulate the association and dissociation of small oligonucleotides (16 base pairs). The model further allows the reconstruction of a fully atomistic representation of an RNA molecule from its coarse-grained representation. The model was also used to study some transitions in RNA, although a direct link between the parameters and experimental melting temperatures has not yet been made.
The above mentioned models were parametrized to structure, either through comparison to experiment, or to atomistic simulations from which thermodynamic quantities are hard to extract. While that is useful for the structure of folded RNA complexes, it makes it hard to compare with available experimental data on RNA thermodynamics, or to simulate reactions involving multiple RNA strands. The next set of models do include explicit thermodynamic information in their parametrization.
The coarse-grained model of Ding et al. 46 uses three beads (sugar, phosphate and base) to represent each nucleotide and has been used to study the folding of various RNA structures, including tRNA and pseudoknots, of sizes up to 100 nucleotides. The parametrization of the interactions combines a knowledge-based approach with a parametrization of the interaction strengths to the free energies of base pairs taken from the nearestneighbor model. Their simulation algorithm furthermore takes into account explicitly the free-energy cost for closing a loop as predicted by the nearest-neighbor model. This added free-energy contribution does not come from the model's interactions and hence ties the use of the model to this particular simulation algorithm.
The model of Hyeon and Thirumalai 47 also uses three beads per nucleotide. Its interaction strengths are based on nearest-neighbor model parameters. The model was used to study mechanical unfolding of hairpins. Recently, the model was extended by Denesyuk and Thirumalai 48 with interactions parametrized using thermodynamic data from pseudoknot and hairpin melting experiments combined with the free energies in the nearest-neighbor model for RNA thermodynamics. 19 The new model has been used to study the thermodynamics of folding of a 34-nucleotide pseudoknot. The model also includes explicit electrostatic interactions and can also represent tertiary structure contacts such as hydrogen bonds in noncanonical base pairs. We note that Hyeon et al. also developed the SOP model for RNA, 49 which only uses one site per nucleotide, to study larger systems. The interactions in the model were set to a given energy scale and were not compared with RNA thermodynamics. The SOP model was used to study the mechanical unfolding of a 421-base ribozyme. 49 The nearest-neighbor model was also used to parametrize the lattice-based model of Cao and Chen 50 which represents the conformations of RNA as a selfavoiding walk on a 3D-lattice. This model was used to compute the heat capacities from partition functions for different mutants of the so-called 72 RNA structure, which were found to be in good agreement with experimental measurements. It was further used to study the free-energy landscape at different temperatures for a 76nucleotide P5abc RNA structure.
Finally, the lattice model of Jost and Everaers 51 is parametrized to reproduce the nearest-neighbor model thermodynamics. The parametrization was verified by studying thermodynamics of ten RNA hairpins and an ensemble of structures with varying internal loop sizes. The model was then used to study folding pathways of a 76 nucleotide long tRNA and a pseudoknot. While lattice based models allow for an efficient sampling of the possible conformations, the structural description of the RNA is necessarily limited by the requirement that it is placed on a lattice.
Most of the existing coarse-grained RNA models are aimed at the correct prediction of the most probable folded structure for a given RNA sequence. In these cases, the thermodynamics of RNA was either not considered, or was used to guide parameter choice which was then tested on a few selected systems. Of the described models, the most detailed verification of the thermodynamics was for the model of Jost and Everaers. We further note that mechanical properties have not been reported for any of these models.
Here we propose a new off-lattice coarse-grained RNA model, oxRNA, that follows the coarse-graining approach developed for the DNA model oxDNA. [52][53][54] Given that oxDNA has been successfully used to model DNA nanotechnological systems, such as motors, 55,56 tweezers, 57 kissing hairpins, 52,58 strand displacement, 59 as well as for biophysical applications including cruciforms, 60 the pulling of single 52 and double-stranded DNA 61 and the hybridization of short oligonucleotides, 62 our goal is to derive a model of similar applicability for RNA. We aim to capture basic RNA structure, mechanics and thermodynamics in as minimalistic a description as possible. We replace each RNA nucleotide by a single rigid body with multiple interaction sites. The interactions between rigid bodies are parametrized to allow an A-helix to form from two single strands and to reproduce RNA thermodynamics. The resultant model goes beyond nearest-neighbor thermodynamics because it has the ability to capture topological, mechanical and spatial effects and allows for the study of kinetic properties of various processes within a molecular dynamics simulation framework. The model uses only pairwise interactions to facilitate the use of cluster Monte Carlo algorithms for simulations. The simple representation, one rigid body per nucleotide, allows for efficient simulation of structures of sizes up to several hundred nucleotides on a single CPU as well as of rare events such as the dissociation or the formation of a double helix.
This paper is organized as follows. In Section II we present our coarse-grained model and its parametrization. In Section III we study the thermodynamic, structural and mechanical properties of the model. We then illustrate the utility of the model through applications to pseudoknot thermodynamics, hairpin unzipping and kissing hairpins in Section IV. The detailed description of the interaction potentials in our model is provided in the Supplementary Material.

II. THE RNA MODEL AND ITS PARAMETRIZATION
In this section we describe the parametrization of the oxRNA model to reproduce the RNA thermodynamics of the nearest-neighbor model, which we briefly outline in Section II A.

A. RNA thermodynamics and the nearest-neighbor model
In an extensive series of investigations, [17][18][19]22 Turner et al. parametrized a nearest-neighbor model (hereafter referred to as the NN-model) to describe the thermodynamics of RNA duplex and hairpin formation, which is widely used in RNA secondary structure prediction. 23,27,[63][64][65] The model treats RNA at the level of secondary structure, estimating enthalpic and entropic contributions to the stability from each pair of consecutive base pairs (bp) in a structure and including corrections for end effects and enclosed loops of unpaired bases. The parametrization used melting experiments of short duplexes and hairpins at 1 M [Na + ]. The results were fitted using a two-state assumption in which RNA either adopts the fully-formed structure or is completely disor- FIG. 1. A schematic representation of (a) an A-RNA helix as represented by the model and of (b) the attractive interaction in oxRNA. The nucleotides can also interact with excludedvolume interactions.
dered. The yield of the duplexes is then given by where [C] is the concentration of closed strands, and [O] is the concentration of open strands. The melting temperature of a duplex, at a given strand concentration, is defined as the temperature at which half of the duplexes present in the solution are dissociated into single strands. Similarly, the melting temperature of a hairpin is defined as the temperature at which the strand has a 50% probability of being open. The NN-model has been shown to reproduce the melting temperatures of RNA oligonucleotides with Watson-Crick base-pairing with 1.3 • C accuracy. 17 In our work, we will treat the NN-model as an accurate fit to the melting data and use its melting temperatures predictions for fitting the oxRNA model. The average accuracy (the fraction of correctly predicted base pairs in known secondary structures) of the NN-model is reported to be 73% for a database of strands that has in total 150 000 nucleotides, containing domains of up to 700 nucleotides. 18

B. The Representation
OxRNA uses a single rigid body with multiple interaction sites to represent a nucleotide. Each rigid body has a backbone, 3 ′ -stacking, 5 ′ -stacking, cross-stacking, and hydrogen-bonding interaction sites. The detailed description of the representation of a nucleotide is given in the Supplementary Material (Fig. A2). In the figures we use a schematic ellipsoid to represent the stacking and hydrogen-bonding sites as this allows the orientation of the base to be clearly seen. The potentials between the nucleotides are effective interactions that are designed to capture the overall thermodynamic and structural consequence of the base-pairing and stacking interactions, rather than directly representing the microscopic contributions such as electrostatics, dispersion, exchange repulsion and hydrophobicity.
We choose the functional forms of our coarse-grained interactions to reproduce directly experimentally measured properties of RNA. For these reasons, we label our coarse-graining approach "top-down", as opposed to a "bottom-up" approach which starts from a more detailed description of the system. We point out that any coarsegrained interaction is actually a free-energy for the real system, rather than a potential energy, and therefore it is in principle state-point dependent. So it should come to no surprise that our potential contains an explicit dependence on the temperature, although for simplicity we try to limit this as much as possible and only introduce it in one of the interaction terms (V stack , as we will discuss later). Our coarse-graining aims to retain the relevant geometric degrees of freedom in order to still correctly capture the relative entropies of different states, despite not having temperature dependence in most of the interaction potentials. 66 The potential energy of the model is where the first sum is taken over all the nucleotides that are neighbors along an RNA strand and the second sum is taken over all the non-nearest-neighbor pairs of nucleotides. All potentials are two-body potentials. There is a maximum distance beyond which all potentials are zero (with the exception of V backbone which diverges to infinity as the distance between adjacent backbone sites approaches its maximum value). The interactions are schematically shown in Fig. 1. We discuss briefly the potentials here while the detailed description is given in the Supplementary Material A 2 b. The backbone interaction, V backbone , is an isotropic FENE (finitely-extensible nonlinear elastic) potential and depends only on the distance between the backbone sites of the two adjacent nucleotides. This potential is used to mimic the covalent bonds in the RNA backbone that constrain this intramolecular distance. The nucleotides also have repulsive excluded-volume interactions V exc and V ′ exc that depend on the distance between the interaction sites, namely the backbone-backbone, stacking-stacking and stacking-backbone distances. The excluded-volume interactions ensure that strands cannot overlap, or pass through each other in a dynamical simulation.
The duplex is stabilized by hydrogen bonding (V H.B. ), stacking (V stack ) and cross-stacking (V cross st. ) interactions. These potentials are highly anisotropic and depend on the distance between the relevant interaction sites as well as the mutual orientations of the nucleotides. The anisotropic potentials are of the form where the functions f are a product of multiple terms, one of which depends on the distance between the relevant interaction sites and the remaining are angular modulation functions that are equal to one if the relevant angles between the nucleotides correspond to the minimum potential energy configuration, and smoothly go to zero as they depart from these values. The set of angles is different for each potential and includes angles between intersite vectors and orientations Ω i and Ω j of the nucleotides. The constant prefactors α ij , η ij , γ, and µ set the strength of the interactions, with α ij , η ij being dependent on the nucleotides involved. The hydrogen-bonding term V H.B. is designed to capture the duplex stabilizing interactions between Watson-Crick and wobble base pairs. The potential reaches its minimum when two complementary nucleotides (AU, GC or GU) are coplanar, directly opposite and antiparallel with respect to each other and at the right distance.
The stacking interaction V stack mimics the favorable interaction between adjacent bases, which results from a combination of hydrophobic, electrostatic and dispersion effects. It acts only between nearest-neighbor nucleotides and its strength depends on both the distance between the respective 3 ′ and 5 ′ stacking sites of the nucleotides as well as their mutual orientations. It also depends on the vector between the backbone interaction sites in a way that ensures the inclination of the bases in the duplex structure matches that for A-RNA. We note that the nucleotides can also interact via the stacking interaction when they are in the single-stranded state. To ensure the right-handedness of the RNA helix in the duplex state, the stacking interaction has an additional modulation term that is equal to one if the nucleotides adopt a right-handed conformation and goes smoothly to zero in the left-handed conformation.
Similarly to oxDNA, the interaction strength of the stacking potential has a temperature-dependent contribution (the term κ k B T in Eq. 5). This term was introduced in oxDNA in order to correctly reproduce the thermodynamics of the stacking transition. We also found that retaining this temperature dependence enables oxRNA to reproduce more accurately the widths of the melting transitions, which are discussed in more detail in Section II D.
The cross-stacking potential, V cross st. , is designed to capture the interactions between diagonally opposite bases in a duplex and has its minimum when the distance and mutual orientation between nucleotides correspond to the arrangement of a nucleotide and a 3 ′ neighbor of the directly opposite nucleotide in a A-helix. This interaction has been parametrized to capture the stabilization of an RNA duplex by a 3 ′ overhang. 19 OxRNA does not include any interaction with the 5 ′ neighbor of the directly opposite nucleotide, as 5 ′ overhangs are significantly less stabilizing than 3 ′ overhangs. 19 Finally, the coaxial stacking potential V coaxial st. represents the stacking interaction between nucleotides that are not nearest-neighbors on the same strand.
We note that, although our model does not include an explicit term for electrostatic interactions between phosphates, these interactions are effectively incorporated into the backbone repulsion. We chose to parametrize our model to the experimental data at 1 M [Na + ], where the electrostatic interaction is highly screened, making our approach reasonable. Furthermore, we are only able to capture those tertiary structure motifs that involve Watson-Crick and wobble base pairing or stacking, such as kissing hairpins or coaxial stacking of helices. In particular, oxRNA does not include Hoogsteen or sugar-edge hydrogen-bonded base pairs, or ribose zippers (interactions involving the 2 ′ -OH group on the ribose sugar). In principle these interactions could be included, but for this version of the model we chose not to as there are no systematic thermodynamic data to which we could parametrize the relevant interaction strengths.
While the strengths of the hydrogen-bonding and stacking interactions depend on the identities of the interacting nucleotides, as in oxDNA, all nucleotides in oxRNA have the same size and shape. Therefore we do not expect oxRNA to capture detailed sequencedependent structure of the A-helix.
The positions of the minima in the potential functions have been selected so that the model reproduces the structure of the A-RNA double helix, which RNA duplexes have been shown to adopt 1,2 and which we describe in more detail in Section III A. The widths of the potential functions and the strengths of the interaction potentials were parametrized to reproduce RNA thermodynamics as described in Section II D.

Algorithms
For the majority of our simulations, unless noted otherwise, we use the Virtual Move Monte Carlo algorithm (VMMC) (specifically the variant described in the Appendix of Ref. 67) to obtain the thermodynamics of oxRNA. VMMC is a cluster-move Monte Carlo algorithm that has previously been used in many applications for the oxDNA model. 66 We also implemented the oxRNA forcefield in a molecular dynamics (MD) code with the choice of a Langevin 68 or an Andersen-like 69 thermostat. For the MD simulations in this work, we used the Andersen-like thermostat.
In our simulations, we often combine the VMMC algorithm with the umbrella sampling method 70 in order to help systems overcome free-energy barriers. This technique involves splitting the configuration space of the system into regions using order parameters, and then assigning weights to each state so defined. In our case, the order parameter is typically the number of base pairs in a given system. The probability of accepting a move to a state with a different weight is then scaled by the ratio of the weights of the two respective states. By assigning larger weights to less probable states (such as states with only a few base pairs), one can achieve an efficient sampling of all values of the order parameter. When extracting the occupancy probabilities for different states from our simulation, we unbias them by rescaling by the inverse of the weight of each state. The weights are typically chosen by experience and then adapted iteratively by hand, until a selection of weights that samples efficiently all relevant states of the system is obtained.

Calculation of melting temperatures
When calculating the free energy of a system as a function of an order parameter (typically the number of base pairs between two strands or, for example, the number of base pairs in a hairpin stem), we run multiple simulations that sample all relevant states, aided by using umbrella sampling. We then calculate p i , the unbiased probabilities of the system having a particular value i of the order parameter, and take the free energy of such a state to be F i = −k B T ln (p i ) + c, where k B is the Boltzmann constant and c is a constant. We are interested only in differences in free energies between different states and hence c can be set to an arbitrary value.
When calculating the melting temperature of a duplex, we simulate a system of two complementary strands and calculate the ratio Φ of the times that the system spends in a duplex state and in a dissociated state. At the melting temperature T m , Φ = 2 for heterodimers, accounting for finite-size effects that come from simulating only two strands, as explained in Refs. 71, 72 and 73. For transitions involving single strands, such as hairpin formation, the ratio Φ is the time spent in closed states (i.e. states with at least one base pair in the stem) divided by the time spent in open states. In this case, Φ = 1 at the melting temperature because the finite size effects mentioned above are irrelevant for unimolecular assembly. We note that we consider the system to be in a duplex (hairpin) form if the complementary strands (stems) have one or more base pairs present.
Using the histogram reweighting method, the states generated during the umbrella sampling simulations can be used to calculate the melting temperatures for a model with the same functional forms of the potential but with different interaction strengths. Previously, we used this approach for the parametrization of the oxDNA model and the method is described in detail in Ref. 52. The reweighting method recalculates the ratio Φ for temperature T and a given set of interaction strengths by counting the contribution to the bonded or unbonded state of each sampled state with an additional weight factor where V 0 is the potential energy of the state generated in the simulation with the original set of parameters at temperature T 0 , and V (T ) is the potential energy of the same state recalculated for the interaction strengths for which we want to find the ratio Φ at temperature T . The reweighting method allows us to extract Φ for a given set of interaction strengths at a range of temperatures which then can be interpolated to find the melting temperature. This method assumes that the ensemble of states generated with the original set of parameters is representative of the states that the system would visit in a simulation with the new parameters. Using this approach it is possible to calculate the melting temperatures for thousands of different sets of interaction constants in an hour of CPU time without the necessity of rerunning the umbrella sampling simulation, which, by contrast, can take several days on a single CPU for each sequence to calculate the melting temperature to within 1 • C accuracy.

D. Parametrization of the model
The anisotropic potentials V stack , V H.B. and V cross st. have interaction strengths of the form η ij (1 + κ k B T ), α ij and γ respectively, where the stacking interaction strength depends also on the simulation temperature T and i and j correspond to the types of interacting nucleotides (A, C, G, U). The magnitude of the temperature dependence of the stacking, κ, and the cross stacking interaction strength γ are set to be the same for all nucleotide types.
In the first step of the fitting procedure, we parametrize the model to reproduce the melting temperatures of the averaged NN-model, for which we define the enthalpy and entropy contribution per base-pair step by averaging contributions of all possible Watson-Crick base-pair steps in the NN-model. In calculating average melting temperatures of different motifs, such as hairpins or terminal mismatches and internal mismatches, the additional entropy and enthalpy contributions for a particular motif in the NN-model were again averaged over all possible combinations of bases. In the averaged NNmodel, the melting temperature is hence independent of the particular sequence, but depends only on the lengths of the sequence and the particular secondary structure motif.
The fitting of the interaction strength parameters was done by a simulated annealing algorithm, which aims to find a set of parameters that minimizes the sum of absolute differences between the melting temperatures of a set of systems as calculated by oxRNA and as predicted by the NN-model. The algorithm uses the reweighting method outlined in the previous section on an ensemble of states generated by a VMMC simulation to calculate melting temperatures for a particular set of parameters. This algorithm for model parametrization is described in detail in Ref. 52.
First, the oxRNA model was parametrized to reproduce averaged NN-model melting temperatures of structures with only Watson-Crick base pairs. The interaction strength η ij was hence set to η avg for all base pair types i and j and α ij was set to α avg for Watson-Crick complementary nucleotides and 0 otherwise. The initial values for α avg , η avg and γ were first chosen by hand and then refined based on results of VMMC simulations in order to reproduce melting temperatures as predicted by the averaged NN-model of short duplexes of lengths 5, 6, 7, 8, 10 and 12 bp and of duplexes of lengths 5, 6, 8 bp with one overhanging nucleotide at either both 3 ′ ends or both 5 ′ ends. We set κ to be equal to 1.9756 (in the inverse of the energy unit used by our simulation code, as defined in Supplementary Material A 2), the same value as was used by the previous oxDNA model. 52,53 We found that leaving κ as a free fitting parameter did not lead to a significantly better fit to the considered sequences and motifs.
We note that for some applications, where one is more interested in the qualitative nature of the studied system or one wants to average over all possible sequences, it might be more useful to study the system with a sequence independent model. We refer to such a model as the "average-base" model meaning that η ij are set to η avg for all types of bases and α ij are set to α avg for all Watson-Crick complementary base pairs (GC and AU) and 0 otherwise. If one is interested in sequence-specific effects, then a more complex parametrization is necessary.
We used the final values of η avg , α avg as the initial values for fitting the sequence-dependent values η ij , α ij , with i and j being Watson-Crick or wobble base pairs (AU, GC, GU). The parameters were fitted to an ensemble that contained oligomers and hairpins of the above mentioned sizes, with 100 randomly generated sequences (with only Watson-Crick base pairs) for each size, and 533 further random duplexes of lengths 5 to 12 bp containing GU wobble base pairs. We excluded sequences with neighboring wobble base pairs in the fitting process as these can lead to duplexes with particularly low melting temperatures (some of the base pair steps containing wobble base pairs are actually destabilizing at room temperature 19 ). We found that our model was unable to accurately fit melting temperatures of duplexes that contain neighboring GU/UG or UG/GU wobble base pairs, probably due to the fact that we do not account for the structural changes that these induce in the duplex.
We note that if one included only Watson-Crick base pairs in the sequence-dependent fitting (as was the case for the parametrization of the oxDNA model 52 ), it would not be possible to distinguish between certain stacking interaction types. For instance, the contribution of AA and UU base stacking interactions always appear together in the AA/UU base pair step free-energy contribution in the NN-model. However, including wobble base pairs in the fitting ensemble provides additional information, for example the UU stacking contribution also appears in the AG/UU base-pair step. We therefore do not need to restrict the strength of stacking interaction to be the same for certain types of nucleotides, as was the case for the oxDNA model.
Finally, we parametrized the coaxial stacking interaction potential, V coaxial st. , which captures the stacking interaction between two bases that are not neighbors along the same strand. Experiments have measured this interaction by a comparison of the melting of a 4-base strand with its complement, or with a hairpin with a 4-base overhang with the complementary sequence adjacent to the hairpin stem. They found for both DNA and RNA that the melting temperature increases for the 4bp long strand attached to the overhang on the hairpin stem, which was attributed to the extra stabilizing interactions with the adjacent stem. 20,21,74,75 The coaxial stacking free energy has been incorporated in the NN-model by assuming that the free-energy stabilizations in these experiments are similar in strength to the actual base-pair steps with the same sequence. The NN-model hence uses the same free energy contribution for a base pair coaxially stacked on a subsequent base pair (as illustrated in Fig. 1) as it uses for a base pair step in an uninterrupted duplex. In order to parametrize these interactions for oxRNA, we performed melting simulations of a 5-base strand, which was able to associate with a complementary 5 ′ overhang on a longer duplex (which itself was stable). We fitted the interaction strength µ of the coaxial stacking interaction in our model so that it would match the prediction of the melting temperature by the averaged NN-model. We note that in our model, the contributing factors to stabilization are both the coaxial stacking interaction and the cross-stacking interaction between the 5-base strand and the hairpin.

III. PROPERTIES OF THE MODEL
In this Section, we describe the structural properties of the model and report the thermodynamics of duplexes, hairpins and other secondary structure motifs as represented by the model. We further study some of its mechanical properties, namely the persistence length of a duplex, the force-extension curve for duplex stretching, and the overstretching transition. A. Structure of the model As mentioned in Section II D, the coarse-grained interactions were selected so that the model reproduces the A-form helix that RNA duplexes have been shown to adopt at physiological conditions. 2, 76 The A-RNA structure is significantly different from B-DNA, the prevalent duplex structure found in DNA molecules. These differences are mainly caused by the sugars in A-RNA adopting a more twisted conformation (C3 ′ endo pucker) as a result of the presence of an extra OH group on the sugar. The A-RNA duplex has a reported helical twist ranging 76 from 32.7 • to 33.5 • per base pair, corresponding to a pitch of 10.7 to 11 base pairs. The rise per base pair reported by X-ray measurements 76 is about 0.28 nm. The bases are displaced from the helical axis, i.e. the helical axis does not pass through the base pair mid-points as it is approximately the case for the B-DNA helix. Finally, the bases are not perpendicular to the helical axis, but are inclined at an angle of about 15.5 • . Although the width of A-RNA is reported to be about 2.1 nm from X-ray crystal structures, 77 Reference 78 uses an effective hydrodynamic diameter of 2.8 nm for the structure. The A-RNA helix has a narrow major groove (0.47 nm) and a wide minor groove (1.08 nm).
To characterize the structure of the oxRNA duplex, we simulate a 13-base-pair duplex at 25 • C using Monte Carlo simulation. We generated 30 000 decorrelated configurations that were analyzed in the following manner. The helical axis was fitted for each saved configuration, as described in Supplementary Material A 1. The rise per base pair was measured as the distance between the projections of the midpoints of base pairs onto the helical axis. The length scale in the oxRNA model is defined so that the rise per base pair is 0.28 nm. The twist per base pair was measured as the angle between the projections of the vectors connecting bases in the base pairs onto the plane perpendicular to the helical axis. The mean turn per base pair in the model is 33.0 • , corresponding to a pitch of 10.9 base pairs. The inclination, measured as the mean angle between the vector pointing from the center of mass of a nucleotide to its base and the plane perpendicular to the helical axis, is 16.1 • .
The width of the helix is measured as twice the distance of the backbone from the axis, and includes the excluded volume interaction radius of the backbone site. The helix width in oxRNA is 2.5 nm. The major and minor grooves in oxRNA are 0.48 nm and 1.07 nm, respectively, where we measured the groove distances in a manner analogous to a method employed by the Curves+ software 79 for analyzing atomistic structures of DNA and RNA. For a selected nucleotide, we measured distances between its backbone site and points on a curve that was linearly interpolated through the backbone sites of the nucleotides on the opposite strand. The distances measured along the curve have two minima, one for each groove. The excluded volume interaction radius for each backbone site was subtracted from these measured distances.

B. Thermodynamics of the model
In this section, we examine the thermodynamics of duplexes, hairpins, bulges, and internal and terminal mismatches as represented by oxRNA. We compare the melting temperatures as predicted by oxRNA with the melting temperatures calculated from the NN-model (denoted as T m (NN)) for different sequences and different secondary structure motifs. To calculate the melting temperatures, we used the reweighting method, described in Section II C 2, with the states that were generated from VMMC simulations of melting for the average-base parametrization of oxRNA.

Duplex and hairpin melting
A comparison of the melting temperatures of the average-base parametrization of oxRNA with the thermodynamics of the averaged NN-model for structures involving only Watson-Crick base pairs is shown in Table  I. For this averaged model, the differences are roughly on the order of the accuracy of the NN-model itself.
To test the sequence-dependent parametrization of the hydrogen-bonding and stacking strength of the interactions, we calculated the melting temperatures for randomly generated ensembles of RNA duplexes, different from the ensemble used for parametrization. A histogram of the differences in the melting temperature pre-  Figure 2(a). The main histogram is for duplexes with both Watson-Crick and wobble base pairing, but not containing GU/UG or UG/GU base pair steps. For convenience, the generated ensembles of sequences also do not include any self-complementary sequences because the calculation of their melting temper-atures requires a different finite size correction. 71,72 The average difference in melting temperatures is 2.1 • C, with an average absolute deviation of 3.4 • C. The histogram in the inset of Fig. 2(a) is for sequences containing at least one GU/UG or UG/GU base pair step. The average difference in melting temperatures for the ensemble is 9.4 • C and the absolute average deviation is 9.7 • C. We note that in the NN-model, the free-energy contribution of these base pair steps is positive at 37 • C, meaning that they actually destabilize the duplex. However, in the oxRNA model, the cross-stacking, stacking and hydrogen-bonding interactions are always stabilizing interactions and the interaction strength of two hydrogenbonded nucleotides does not depend on the identity of their respective neighbors on the strand. Our coarsegrained model hence cannot capture the free-energy contributions of GU/UG and UG/GU base pair steps. One could imagine adding multi-body interactions, but for the sake of computational efficiency and maintaining the consistency of our coarse-graining methodology, we do not do so in this study. Another option might be to introduce a structural perturbation of the helix caused by the GU base pairs. The histogram in Fig. 2(b) shows the difference between melting temperatures calculated by oxRNA and those predicted by the NN-model for an ensemble of randomly generated hairpins. The average melting temperature difference was −2.8 • C with the average absolute deviation being 4.1 • C.
The transition widths for duplex and hairpin formation were calculated for the averaged model as the difference between the temperatures at which the yield is 0.8 and 0.2, respectively. This quantity is important because the widths of the transition determine the change of the duplex melting temperature with concentration. It can be shown 53 that the derivative of the melting temperature as a function of strand concentration is proportional to the width of the transition. The melting simulations of oxRNA with the average-base parametrization were compared with the width predicted by the averaged NN-model. For the duplexes of lengths 6, 8, 10, and 12 the width was on average underestimated by 0.9 • C, but was overestimated for a 5-mer by 0.5 • C. The width of the transition for the averaged NN-model decreases from 20.5 • C for a 5-mer to 9.2 • C for a 12-mer. For a set of hairpins with stems of length 6, 8, and 10 bp and with loops of lengths 5, 6, 7, 8, 10 and 15, the width of the melting transition was on average underestimated by 1.5 • C. The width of the hairpin transition decreases from about 12 • C for stems of length 6 to approximately 8 • C for stems of lengths 10 in the averaged NN-model. However, the trend of increasing width with decreasing size is always captured by the oxRNA model.
Finally, we note that we could have parametrized the sequence-dependent model only to duplex melting temperatures, as for oxDNA, which would then have led to a larger underestimate of hairpin melting temperatures, presumably because our representation of the strand is too simple to exactly capture the entropy and enthalpy of the loop formation. We hence chose to parametrize to the ensemble of duplexes and hairpins because hairpins are a prominent secondary structure motif in RNA.

Thermodynamics of secondary structure motifs
Given that our aim is to design a model that goes beyond describing hybridization of simple duplexes, it is important to assess how well it is able to reproduce the thermodynamic variation induced by common secondary structure motifs such as bulges or mismatches. To this end, we have studied the melting temperature changes induced by internal mismatches, terminal mismatches and bulges of different lengths.
To assess the effects of bulges, we consider a systems of two strands, one with 10 bases and the other with 10 complementary bases and extra bases that create a bulge motif. We considered bulges of lengths 1, 2, 3 and 6, positioned in the center. For each sequence considered we calculated the melting temperatures by reweighting a set of 6000 states that were sampled from a melting simulation using the average-base parametrization. For each bulge length, we considered 1000 randomly generated sequences with Watson-Crick base pairing in the complementary part.
We further evaluated the melting temperatures for randomly generated 10-mers which contained 1, 2 or 3 consecutive mismatches (and therefore had 9, 8 or 7 complementary Watson-Crick base pairs). The average difference and absolute average deviation for the considered bulges and mismatches are shown in Table II. The melting temperatures of duplexes with bulges are underestimated by a few degrees. However, the model presently significantly overestimates the stability of internal mismatches by the order of 10 • C or more. Even though the mismatching base pairs do not gain stabilization from hydrogen-bonding interactions (which are zero for bases that are not complementary), they still have favorable cross-stacking and stacking interactions, which have their minimum energy configuration in an A-helical configuration, which the oxRNA model can still form with the mismatches presents. We further note that our model represents each nucleotide by the same rigid body structure, so the mismatching base pairs do not cause any distortion to the duplex structure in our model. Improving the model to more accurately represent the secondary structures with mismatching nucleotides will be the subject of future work. The persistence length L p of dsRNA molecule measured in experiments is reported to be between 58 nm are the average differences in melting temperature between the sequences with a secondary structure motif and a duplex with the same sequence but with no bulge or internal mismatch as predicted by the NN-model and oxRNA respectively. Each of the motifs considered is destabilizing, resulting in a decrease of the melting temperature. The averages were taken over an ensemble of randomly generated sequences (1000 for each motif) that had 10 complementary Watson-Crick base pairs for the bulges, and 9, 8, and 7 complementary base pairs for internal mismatches of size 1, 2 and 3 bases, respectively. The bulges that we consider were of the size 1, 2, 3 and 6 bases. All the melting temperature calculations were calculated for an equal strand concentration of 3.36 × 10 −4 M. To measure the persistence length in our model, we simulated an 142-bp long double-stranded RNA with the average-base oxRNA model, and measured the correlations in the orientation of a local helical axis along the strand. The local axis vectorl i is fitted through the i-th base pair and its nearest neighbors, using the approach described in Supplementary Material A 1, but considering only the nearest neighbors. We found the results to be robust even when 2 or 3 next nearest neighbors were included in the construction of the local axis. To account for edge effects, a section of ten base pairs at each end of the duplex was ignored when accumulating averages. For an infinitely long, semi-flexible polymer in which the correlations decay exponentially with separation along the strand, the persistence length L p can be obtained from 83 where r is the rise per base pair. This correlation function is shown in Fig. 3 along with the exponential fit from which we estimated the persistence length of our model to be about 101 base pairs, somewhat lower than the 171 bp persistence length expected at this salt concentration. Our model's persistence length is hence smaller than the experimentally measured values, but still within a factor of 2. OxRNA hence captures the correct order of magnitude for the persistence length and as long as one studies structures whose duplex segments are smaller than the persistence length of the model, it should provide a physically reasonable description. We note that the persistence length is quite hard to correctly reproduce. We expect this issue to hold for other coarse-grained models as well. To our knowledge, the persistence length has not been measured yet in these models.

Force-extension properties
As a further test of the mechanical properties of the model, we measured the extension of a 99-bp RNA duplex as a function of applied force for the average-base model. We used a constant force acting on the center of mass of the first and last nucleotides in one of the two strands in the duplex. We focus on the average extension between the 11 th and 91 st nucleotide on this strand in order to avoid end effects, such as from fraying of base pairs, in our measurements.
We compare our force-extension data with an extensible worm-like chain model, 84 which provides the following expression for the projection of the end-to-end distance R along the direction of the forceẑ: where K is the extension modulus and L 0 is the contour length. This expression comes from an expansion around the fully extended state, and thus it is expected to be valid at forces high enough for the polymer to be almost fully extended. It was shown experimentally 82 that this extensible worm-like chain model describes the behavior of dsRNA prior to the overstretching transition. In particular, at 0.5 M [Na + ], the extensible worm-like chain model fit to the experimentally measured force-extension curve yielded L p = 57.7 nm and K = 615 pN.
The force-extension curve for oxRNA is shown in the inset of Fig. 3. We used data from forces between 2.4 pN and 69 pN for our fitting and allowed L 0 , K and L p to be free parameters. From the fit we obtained L 0 = 23.4 nm (84 bp), K = 296 pN and L p = 26.0 nm (93 bp). We note, however, that Eq. 9 is not a robust fit for our model: changing the fitting interval and thus including or excluding points at either high or low forces significantly changes the resulting values of the fitting parameters, even though the residual error in the fit remains small. The persistence length, for instance, can change by more than a factor of two. In the interval we selected for the fit, the L p obtained approximately corresponds to that obtained from the correlation function in Fig. 3, but given the sensitivity of the fit (which was not observed for the oxDNA force-extension curve 53 ), the errors on these extracted parameter values should be taken to be quite large. Another issue to keep in mind is that the inclination angle is also affected by the force. At 10 pN, this is roughly a 1 to 2 degree change, but close to the point where the chain starts to significantly overstretch (as discussed in the next section), the inclination has disappeared, and the bases are almost perpendicular to the axis. It is likely that this deformation is not entirely physical. However, the physical structure of stretched RNA is not experimentally known. In DNA, the structure of the extended state is a very active topic of research.

Overstretching
Both DNA and RNA duplexes are known to undergo an overstretching transition at high enough force. Recent experiments 82 for different salt concentrations find 63.6 (2.0) pN for RNA overstretching at 0.15 M [Na + ] up to 65.9 (3.3) pN at 0.5 M [Na + ]. Following the approach taken in the study of DNA overstretching with the oxDNA model, 61 we used the average-base oxRNA model to run VMMC simulations of a 99-bp RNA duplex with equal and opposite forces applied to the first and last nucleotide of one strand. In our simulations, only native base pairs were allowed to form to aid equilibration, i.e. no misbonds in the duplexes or intrastrand base pairs in the unpeeled strand. The simulations were started from a partially unpeeled state to sample states which have between 65 and 71 bp. The obtained freeenergy profiles as a function of the number of base pairs are shown in Fig. 4(a). As the force increases, the slope of the free energy profiles changes from negative (states with more base pairs are favored) to positive (it is favorable for duplexes to unpeel). By estimating the force at which the slope becomes zero, we obtained 86.2 pN as the overstretching force. We note that our model was parametrized for 1 M [Na + ], whereas the overstretching experiment was done at 0.5 M. Furthermore, by not allowing any formation of secondary structure in the unpeeled strands, we overestimate the overstretching force in the model, because these intramolecular base pairs stabilize the unpeeled state. For the oxDNA model, it was shown that allowing secondary structure decreases the overstretching force by about 3 pN. 61 We would expect the stabilization to be slightly higher for RNA, as RNA base pairs are more stable. Our model hence overestimates the value of the overstretching force by about 16-20 pN. This overestimation is partly due to the higher extensibility of the duplex in oxRNA, which is aided by the loss of inclination in the duplex when higher forces are applied, as we already discussed in the previous section.

A. Pseudoknots
Pseudoknots are a common structural motif in RNA. If a strand is represented as a circle and base pair contacts are represented as chords, then its structure is pseudoknotted if the chords cross. Most secondary structure prediction tools do not include pseudoknotted structures in their computations and thus cannot be used to study systems where they are relevant, although some progress has been made in developing efficient algorithms for this task. 85,86 Since oxRNA provides an explicitly three-dimensional representation of the RNA strands, it can be used to simulate the folding and thermodynamics of RNA structures that involve pseudoknots. In this section, we use our model to study the well known MMTV pseudoknot. 87 The sequence and the three-dimensional representation of the MMTV pseudoknot by oxRNA are shown in Fig. 5. The MMTV pseudoknot's thermodynamic properties were previously studied in experiment, 87 as well as with another coarse-grained RNA model. 48 Moreover, the MMTV pseudoknot's structure has also been investigated by NMR experiments 88 and two stems were identified in the folded structure: stem 1 with 6 base pairs and stem 2 with 5 base pairs, as schematically shown in Figure 5.
To study the thermodynamics of the system, we ran VMMC simulations of oxRNA for 3.4 × 10 11 steps at 75 • C. Umbrella sampling, using the number of base pairs in each of the pseudoknot stems as order parameters, was employed to enhance thermodynamic sampling. We also used histogram reweighting to extrapolate our results to other temperatures. The occupation probabilities of the unfolded state, a single hairpin with stem 1 or stem 2 (denoted as hairpin 1 and hairpin 2), and the pseudoknot are shown in Fig. 6(a). Our simulations also allow us to extract the heat capacity C V from the expression where we use a cubic interpolation to our simulation data for U in order to compute the derivative with respect to T . The results are shown in Fig. 6(b). The experimentally measured C V at 1 M [Na + ] has two peaks, one at 73.5 • C and the other at 95.0 • C. 87 It was hypothesized that the two peaks correspond to the transition from an unstructured strand to a hairpin structure and a second transition from a hairpin structure to the full pseudoknot. Analysis of our yields supports this claim, showing a pseudoknot to hairpin 2 transition at 67.7 • C and transition from hairpin 2 to a single strand with no bonds in stem 1 or stem 2 at 84.8 • C. The higher temperature transition coincides with a peak in the heat capacity, whilst the lower temperature transition gives rise to a shoulder. While our simulations reproduce qualitatively the behavior of the experimental system, the position of the transitions is not exactly the same as the ones measured experimentally. This is not surprising, as we have shown in Section III B that the model generally underestimates the melting temperatures of hairpins.
It is of further interest to analyze the free-energy landscape of the system (Fig. 7). Perhaps unsurprisingly, our results suggest that the minimum free-energy pathway for folding this pseudoknot from a single strand involves first forming one of the stems (forming stem 2 first is more likely as it is more stable and has a higher yield at the considered temperature) and then closing the second stem, instead of simultaneously forming both of them. We have previously seen similar pathways for a DNA pseudoknot. 66 One subtlety concerns the formation of the GU base pair between the seventh and thirty-fourth nucleotide. The NMR study 88 did not observe the presence of this GU base pair in the pseudoknot structure. However, in our simulations, we find some structures where this base pair forms (thus extending the size of stem 1 from six to seven base pairs), although it has a 5 k B T free energy penalty at 48 • C with respect to a pseudoknot state which had only six bases in stem 1. Including this additional base-pair within the definition of stem 1 had only a small effect on the calculated yields (the positions of the equal yields points changed by less than 0.3 • C) and we saw at most 0.5 k B T free-energy change for the folding pathways. We thus did not include this extra base pair in the definition of stem 1.
The experimental NMR study 88 of the structure of the MMTV pseudoknot found that the two stems of the pseudoknot are bent with respect to each other at about 112 • , and the AA mismatch between the sixth and the fourteenth nucleotides to be not stacked. As can be seen in Fig. 5, in oxRNA, this mismatch is typically stacked leading to an angle closer to 160 • . We think that this stacking of the stems reflects the overestimation of the stability of mismatches in simpler motifs (see Table II).
In summary, our model is able to describe the thermodynamics of the pseudoknot folding and predict the stabilities of the two stems, supporting the hypothesis that the peak in heat capacity at higher temperature corresponds to the pseudoknot to hairpin 2 transition. The overall secondary structure of the pseudoknot is correct in our model, which also helps to understand the tertiary structure even though we found the angle between the two stems to be larger than the one reported from experiment.

B. Kissing hairpin complex
A kissing complex is a naturally occurring motif in RNA structures 1 and consists of two hairpins that have complementary loops and can thus bind to each other. An example of such a complex as represented by oxRNA is shown in Fig. 8. The kinetics and thermodynamics of forming an RNA kissing complex with 7 bases in the loops was experimentally studied in Ref. 90 at varying salt concentrations, including 1 M [Na + ], the concentration at which our model was parametrized.
To examine the capability of oxRNA to describe kissing complexes, we studied the melting of this kissing complex using both the average-base and the sequence-dependent parametrization of oxRNA and found the transition at a point which is approximately consistent with the observed experimental behavior. The 7-base loops in the two hairpins have fully complementary sequences (5 ′ -GGAAAUG-3 ′ and its Watson-Crick complementary sequence). All melting simulations were run in a volume corresponding to an equal strand concentration of 3.36 × 10 −4 M.
For the average-base representation we found the melting temperature of the kissing hairpins to be 62.7 • C which compared to 53.6 • C for a 7-bp duplex with the same sequences as the loops. For the sequence-dependent model, we found the melting temperature of the kissing complex to be 44.8 • C, similar to 45.2 • C for this 7-bp duplex. The free-energy profiles for both average-base and sequence-dependent models at 45 • C are shown in Fig. 8.
For most sequences, we find that the kissing hairpin loop is more stable than the separate duplexes with respect to the unbound state. We find that with increasing temperature, the kissing complex loses less stability with respect to the unbound state than a duplex at the same temperature and strand concentration. This trend can be rationalized in terms of the fact that a constrained loop loses less configurational entropy upon binding than an unconstrained single strand does. Furthermore, the kissing complex also gets an extra enthalpic stabilization from a coaxial stacking interaction between the loop and the stem nucleotides. These two effects help explain why, on average, the kissing hairpins are more stable, especially at higher temperatures. However, the kissing hairpins do not satisfy the enthalpic contributions as well as a perfectly formed duplex does. Thus, at low temperature, the duplex can be more stable. Which of these effects dominates depends on the sequence, and if the melting happens before the switch of which motif is more stable, then the duplex will have a higher melting temperature, which we find for a minority of sequences at our strand concentration. For the sequence above, the melting temperatures are very close. Note that the hairpin loops are sufficiently short that kissing complex formation is not inhibited by the topological requirement that the linking number between the loops must remain zero. This contrasts with previous simulations of DNA kissing complexes between hairpins that have 20-base complementary loops. 58 The thermodynamics of this kissing complexes was studied in Ref. 90 using isothermal titration calorimetry (ITC) at 1 M [Na + ]. Up to 35 • C the authors found evidence of a transition to a kissing complex after the injection of the reactants, but did not observe a transition at 45 • C. The authors claim to measure only a 0.6 kcal/mol change in the ∆G for forming the kissing complex between 10 • C and 35 • C while observing a significant increase (19.5 kcal/mol) in ∆H along with a compensating increase in ∆S. Such behavior is not observed in our simulations, where we see a classic quasi-two-state transition in which ∆H and ∆S change slowly with temperature, similarly to a duplex association. If the yields of the kissing complexes in our simulations were extrapolated to the concentrations used in the experiment, we would predict a yield of 36% at 35 • C and 5% at 45 • C. We note that the thermodynamic parameters in Ref. 90 were derived with the assumption that the transition was fully saturated after the injection of the reactants in the ITC experiment, which is incompatible with the experimentally inferred value of ∆G. If the transitions were in fact not fully saturated, then it is possible that the experiments are consistent with a more typical quasi-two-state transition as observed in our model with a melting temperature approximately consistent with that found by us.
It is also interesting to use oxRNA to probe the structure of this kissing complex because it is a motif that has been used in RNA nanotechnology. Molecular dynamics simulations of the kissing complex using an all-atom representation (Amber) measured the angle between the helical stems at 300 K and 0.5 M monovalent salt to be approximately 120 • . 89 Based on this finding, a hexagonal ring nanostructure that can be assembled from six RNA building blocks was proposed. Each of the proposed building blocks is an RNA strand that in the folded state has a stem and two hairpin loops. The sequences in the loops are designed to bind to the complementary block to allow the assembly of a hexagonal "nanoring" via the kissing complex interaction. This computationally proposed RNA nanoring design was later experimentally realized by self-assembly. 91 The nanoring can be functionalized by including siRNA sequences either in the hairpin stems or by appending siRNA sequences onto the stems. Experiments in blood serum have shown that the nanoring protects the loop regions of the assembly blocks from single-strand RNA endonucleolytic cleavage, making the nanoring a promising tool for in vivo nanotechnology applications. 91 Simulations of oxRNA at 25 • C allowed us to measure the angle between the helical axes of the hairpin stems. We found this angle to fluctuate around an average value of 133.9 • with a standard deviation of 14.8 • . Hence, an octagon would probably be the most relaxed nanoring for oxRNA, and therefore favored by enthalpy. Smaller rings would be favored by translational entropy.
To illustrate the capabilities of our oxRNA model, we constructed the hexagonal RNA nanoring of Ref. 89 (Fig. 9) by starting a simulation with six folded hairpins blocks and introducing a harmonic potential between the complementary loop regions, which helped the kissing interactions to form. Once the nanoring was completed, the harmonic traps were removed and the assembled structure was relaxed in a molecular dynamics simulation. We found the angle between the stems of the kissing hairpins in the nanoring to fluctuate around a mean of 124.9 • . The thermal fluctuations around the mean value had a standard deviation of 14.4 • , which is similar to that of the single kissing complex. It would also be possible to use our model to study the mechanical properties of the nanoring, as well as the thermodynamics and kinetics of its self-assembly from the building blocks, but such a study is beyond the scope of this article.
A typical relaxation time for the angle between adjacent kissing complexes or the energy self-correlation function of the assembled nanoring structure corresponds to about a minute of CPU time or less on a standard 2.2 GHz processor. This example shows that structural investigations of systems of the order of hundreds of bases are well within the capabilities of the oxRNA model using a single CPU, and if multiple CPUs are used, or a GPU chip is used, 92 then structural properties and fluctuations around equilibrium can be studied for systems on the order of thousands of base pairs, as can be done for oxDNA. 66

C. Hairpin unzipping
RNA hairpin unzipping has been used in experiment to study the thermodynamics of base pairing and the mechanical properties of RNA, with the kinetics of the process also being of interest. [94][95][96] Unzipping the same sequence under different salt and temperature conditions can provide systematic data on the physical properties of RNA that, for example, could be used to improve the parametrization of thermodynamic models of RNA.
Here we consider the unzipping of the CD4 hairpin (shown schematically in Fig. 10(a)), which has a 20-bp stem and 4 bases in the loop. It was studied experimentally by pulling at different rates and measuring the unzipping force 94-96 for each trajectory. While it is possible to simulate pulling the hairpins ends at a given rate in the oxRNA model, direct comparisons with experimentally observed unzipping forces are difficult because for a coarse-grained model there is not a straightforward way to map the simulation time to the experimental one. Furthermore, to obtain an estimate of the average unzipping force for a given pulling rate, one needs to average over multiple trajectories, and generating such trajectories can be quite time consuming, especially for very slow pulling rates.
A more suitable experimental setup for comparison with a coarse-grained simulation comes from Ref. and closed state. By interpolating results for a range of applied forces, they were able to obtain the unzipping force, which was defined as the force for which the freeenergy difference was zero.
To compare to the experimental results, we ran VMMC simulations of the CD4 hairpin with a constant force of 19.7 pN applied to the first and last nucleotide at 23 • C. We then extrapolated the free energies of the open and closed state by the histogram reweighting method to forces ranging from 16.2 to 20.7 pN and to the temperatures at which the experiments were carried out (22,27,32,37 and 42 • C). We only allowed bonds between the native base pairs to avoid sampling of metastable secondary structures that would slow down our simulations. We considered the hairpin to be closed if at least one of the bonds in the stem was present. For each temperature considered, we performed a linear interpolation of the free-energy difference between closed and open state as a function of force to obtain the unzipping force for which the difference is 0. The unzipping forces we obtain are shown in Fig. 10, along with a fit. We also show for comparison the fit to the experimentally observed unzipping forces 96 at 1 M [Na + ], expressed in the form where F unzip is the unzipping force at temperature T . The values obtained in the experiment 96 were a = 69.1 pN and c = −0.164 pN/K. Fitting Eq. 11 to our simulation data, we obtained the same value for c and 68.2 pN for a. The values of the fitting parameters varied by at most 6% between the fits to unzipping forces obtained from three independent sets of generated states. Thus, oxRNA is able to reproduce the unzipping force with 5% (1 pN) accuracy and fully captures the trend with temperature. It would be of further interest to study the forceextension properties of a hairpin which contains various secondary structure motifs such as bulges and internal loops or which has regions with variable sequence strengths (such as AU rich and GC rich regions). However, such a study is beyond the scope of this article.

V. SUMMARY AND CONCLUSIONS
We have presented a new off-lattice coarse-grained model for RNA, oxRNA, which aims to capture basic thermodynamic, structural and mechanical properties of RNA structures with a minimal representation and pairwise interaction potentials. OxRNA is specifically developed to allow for efficient simulations of large structures (up to hundreds of nucleotides on a single CPU), and for reactions involving multiple strands of RNA, which are important for applications in RNA nanotechnology. Our coarse-graining strategy is closely linked to our previous successful coarse-graining of DNA with oxDNA. 53 Rather than focusing mainly on reproducing structure, as many other previous models have done, here we tried to systematically compare to a whole suite of different properties.
We employed a "top-down" coarse-graining approach, where we aim to reproduce free-energy differences between different states (such as opened and closed state of a hairpin) as measured in experiment. OxRNA represents each nucleotide (i.e., sugar, phosphate and base) as a single rigid body with multiple interaction sites. The model is capable of representing RNA structures such as hairpins and duplexes and is designed to reproduce the A-helical form of duplex RNA. We used a histogram reweighting method 52 to parametrize the model to reproduce the thermodynamics of short duplexes and hairpins. Currently, the oxRNA model captures Watson-Crick and wobble base-pairing interactions as well as various types of stacking interaction. However it was not designed to capture non-canonical interactions such as Hoogsteen or sugar-edge hydrogen-bonded base pairs, or ribose zippers. Nevertheless, it can reproduce some important ter-tiary interaction motifs, in particular coaxial stacking of helices, kissing loop interactions, and pseudoknots.
The model is currently parametrized for a salt concentration of 1 M, as this corresponds to the conditions for the melting experiments used for the nearest-neighbor model to which our model was parametrized. Explicit electrostatic interactions are not included, because they are very short-ranged at high salt and thus can be incorporated into the short-ranged excluded volume terms in the potential. This excluded volume also prevents a strand from crossing itself or other strands, forbidding topologically impossible trajectories in kinetic simulations. It is possible to use the same coarse-graining techniques to parametrize the model at lower salt concentrations. However, as the screening lengths become longer, different longer-ranged forms for the interactions may need to be used to capture the correct physics. We note however that nanotechnology experiments in vitro are usually carried out in high salt conditions. To test our model, we investigated the thermodynamics of short duplexes and hairpins with different sequence content, as well as various other secondary structures such as bulges, internal and terminal mismatches. We found that in comparison with our previous model for DNA, parametrizing RNA thermodynamics is a more challenging task. For example, experimental melting temperatures of a duplex of a given length can differ by as much as 70 • C between weak and strong sequences with Watson-Crick base pairs, as opposed to 50 • C for DNA. 52 So sequence effects are stronger in RNA. Including wobble base pairs presents further challenges, as some base pair steps that include two or more neighboring wobble base pairs have a positive contribution to the free energy of duplex. Although it is not possible to capture this effect with the current representation of our model, adding the structural effects of wobble base pairs on the double helix may provide means to improve this aspect of the model. Finally, we found that oxRNA overestimates the stability of duplexes with mismatches in internal loops. This could lead to an overestimation of the stability of misfolded structures and complicate the study of the folding of RNA strands that have multiple metastable states with mismatches. Nevertheless, even though oxRNA does not reproduce the exact melting temperatures for structures with internal mismatches and bulges, we do observe, as expected, a decrease in melting temperatures of a duplex with internal mismatches or bulges with respect to the fully complementary strands. The observed changes in melting temperatures are within the orders of magnitude as predicted by the NN-model for the same motifs and capture correctly the direction of the change.
We have tested the mechanical properties of the RNA duplex as represented by the model and found its persistence length to be about half of the reported persistence length of RNA molecules at high salt conditions. The model hence has the correct order of magnitude for the duplex persistence length and provides sufficiently accurate mechanical behavior for most applications, where individual double helical sections are likely to be much shorter than the persistence length.
In Section IV, we provided some applications of the model that illustrate its strength and potential utility. In particular, we investigated the thermodynamics of pseudoknot folding and the thermodynamics and structure of a kissing hairpin complex. We also showed that oxRNA can be used to model large nanostructures like an RNA nanoring composed of 264 nucleotides on a single CPU. The computational cost of oxRNA is very similar to oxDNA where simulations of a DNA origami motif with 12 391 nucleotides have been achieved. 66 Finally, our model is able to reproduce experimental results for the mechanical unzipping of a hairpin quite closely, to within an accuracy of 1 pN, illustrating oxRNA's potential to study mechanical properties of RNA constructs.
Although we did not describe applications with detailed dynamics (the simulations are typically quite involved and so beyond the scope of this study), we want to emphasize that oxRNA is particularly well-suited for such tasks. For example, oxDNA has been used to study the detailed dynamics of hybridization, 62 toehold-mediated strand displacement 59 and hairpin formation. 66 Studying similar processes would be very feasible for oxRNA. For example, it should be possible to use the model to obtain estimates of the rates of a strand displacement reactions as a function of length of the toehold as well as temperature. OxRNA can further be used to study the self-assembly of nanostructures such as RNA nanorings and to investigate both the thermodynamics and the kinetics of such systems. Although the model is currently only parametrized at high salt concentration, oxRNA can be also used to qualitatively study processes of biological relevance, for instance, the folding pathways of RNA strands or in vivo applications of nanotechnology.
At this point it is interesting to reflect on some similarities and differences between the coarse-graining of oxDNA and oxRNA. Although oxRNA can clearly reproduce a good number of properties of RNA, quantitative differences with experiment for the melting temperatures of certain motifs are larger than they are for oxDNA. Moreover, it was more difficult to simultaneously reproduce the thermodynamics and the correct persistence length or the force-extension curves. One reason for these differences may be that RNA itself exhibits more complex behavior than DNA, and so is harder to coarse-grain. It is intuitively obvious that the compromises made to increase tractability mean that not all properties of the underlying system can be simultaneously captured by a more simplified description, a general phenomenon that has been called "representability problems". 97,98 One consequence of representability problems is that in general, fitting too strongly to one set of input data (say structure, as is often done for other RNA models) will often lead to larger errors in other quantities, say thermodynamics. We tried to compromise between different requirements for oxRNA. However, in order to make further progress, more detailed fitting may not be enough. Instead a more complex and most likely less tractable representation of the full interactions may need to be chosen. For example, for RNA it remains to be seen whether our single nucleotide-level model can be easily extended to generate a better representation of both structure and thermodynamics, or whether, say, a more complex model is needed to achieve the next level of accuracy. Clearly, including electrostatic effects for lower salt-concentrations, or implementing tertiary structure contacts, for example non-canonical base pairing interaction (such as Hoogsteen base pairs) and hydrogen bonding between a sugar group and a base will also need an extension of the current representation.
Given the challenges and complexity of RNA modelling, it is unsurprising that oxRNA performs less well than oxDNA for the whole ensemble of motifs we considered. However, we believe that it is a non-trivial achievement to create a model that can semi-quantitatively reproduce such a wide range of the thermodynamic data. The properties of our model have been comprehensively tested on a variety of systems, ranging from secondary structure motifs to systems such as kissing complexes and hairpin unzipping. Our model is also currently the only RNA model with reported mechanical properties, which were tested by measuring its persistence length, forceextension curve and duplex overstretching.
Finally, the simulation code implementing molecular dynamics and Monte Carlo algorithms for our oxRNA model is available for download at dna.physics.ox.ac.uk. fluctuations and hence the plane has to be fitted through the endpoints of the vectors. The first and last two base pairs were not included in order to exclude end effects. The vector perpendicular to the plane was then taken to be the helical axis. The fitting of the helix is schematically illustrated in Fig. A1.

The potentials and nucleotide representation in the RNA model
The model and its potentials were introduced in Section II and here we provide a detailed description of the interaction potentials and the nucleotides. We first describe the representation of each nucleotide as a rigid body in Section A 2 a and then give the explicit expression for each of the potential terms in Section A 2 b. All the values of the potential parameters are in an internal unit system of the downloadable simulation code where 1 distance unit = 8.4Å and 1 energy unit = 41.4 pN nm = 10 k B T for T = 300 K. In molecular dynamics simulations, we set 1 mass unit to correspond to the average weight of a nucleotide, 321.4 AMU, which gives us the simulation time unit corresponding to 3.06 × 10 −12 s in SI units.

a. Representation
Each nucleotide in the oxRNA model is represented as a single rigid body with multiple interaction sites. Each nucleotide has backbone, hydrogen-bonding, cross-stacking and 3 ′ -and 5 ′ -stacking interaction sites. The position and orientation of each nucleotide is uniquely specified by its center of mass position and the perpendicular unit vectors a 3 and a 1 , where a 1 is a unit vector pointing from the center of mass to the hydrogen-bonding site and a 3 is defined in Fig. A2. In a duplex configuration, the a 3 vector would be pointing towards the 5 ′neighboring nucleotide. We further define a 2 = a 3 × a 1 . The nucleotide as represented by oxRNA is schematically shown in Fig. A2. The small colored circles indicate the position of the interaction sites, while the large circles around hydrogen bonding and backbone sites indicate the interaction radius of the excluded-volume interactions.
The interaction potentials are functions of the distances between the relevant interaction sites as well as the angles between intersite vectors and the respective orientation vectors a 3 , a 1 , p 3 ′ and p 5 ′ , where we define We also define the following vectors which are then used in the definitions of the potentials in the oxRNA model (Eq. 3): • δr backbone : the vector between backbone sites of the nucleotides. If the nucleotides are nearest neighbors, it is pointing towards the nucleotide's 3 ′neighbor's backbone site • δr HB : the vector between the hydrogen-bonding sites of the interacting nucleotides, pointing from the first nucleotide towards the second one • δr coaxial st. : the vector between the coaxial stacking sites of the interacting nucleotides, pointing from the first nucleotide towards the second one • δr stack : the vector pointing from the 3 ′ -stacking site of a nucleotide to the 5 ′ -stacking site of its 3 ′neighbor • δr back−base /δr base−back : the vector pointing from the backbone/hydrogen-bonding site of the first nucleotide to the hydrogen-bonding/backbone site of the second nucleotide We further define the following angles that are used in the potential functions: θ 3 = arccos a 1 · δr HB (A5) θ 5 = arccos a 3 · δr coaxial st. (A7) θ 5 ′ = arccos a 1 · δr stack (A9) where we use the notation b 1 , b 2 , and b 3 to define the orientation vectors of the second nucleotide participating in the interaction (the orientation vectors of the first nucleotide are denoted as a 1 , a 2 , and a 3 ). The vector q 5 ′ corresponds to the p 5 ′ vector of the 3 ′ -neighbor of the interacting nucleotide, i.e. using the same definition as in Eq. A2, but substituting b for a.

b. Potentials
The oxRNA potential consists of a sum of potential functions designed to represent different physical interactions, with some of the potentials being products of multiple potential functions. The functions that are used in the potentials are: • FENE spring (used in V backbone ): • Morse potential (used in V stack and V H.B. ): • Harmonic potential (used in V cross st. and V coaxial st. ): • Lennard-Jones potential (used in excluded volume potentials V ′ exc and V exc ): (A22) • Quadratic modulation terms (used for angular modulation of the anisotropic potentials V H.B. , V cross st. , V stack and V coaxial st. ): • Quadratic smoothing terms for truncation (used in all potentials with the exception of V backbone ) in order to make the potentials differentiable functions that are equal to 0 beyond some specific cutoff distance: The smoothed functions used in the potentials have the following form: • The radial part of the stacking and hydrogen-bonding potentials: f 1 (r, ǫ, d, r 0 , r c , r low , r high ) =          V Morse (r, ǫ, r 0 , d) − V Morse (r c , ǫ, r 0 , d)) if r low < r < r high , ǫV smooth (r, b low , r c,low ) if r c,low < r < r low , ǫV smooth (r, b high , r c,high ) if r high < r < r c,high , 0 otherwise. (A25) FIG. A2. A schematic representation of the nucleotides as represented by the oxRNA model. The red circle indicates the backbone site, the blue circle is the hydrogen-bonding site and the green circle is the coaxial stacking site. The yellow circle is the 3 ′ -stacking site, and the black circle is the 5 ′ stacking site. The unfilled circle from which the a3, a2 and a1 vectors originate is the center of mass. All distances are given in a unit system where 1 distance unit = 8.4Å. The left image shows the projection of a single nucleotide where the a2 vector is pointing towards the reader, and the image on the right shows a projection where the a3 vector is pointing towards the reader. For comparison, we also show the schematic representation of the nucleotide that is used in producing pictures of oxRNA configurations. The backbone site is represented by a sphere because of the isotropic nature of the its interactions, whereas the base is represented by an ellipsoid whose principal axes are parallel to a1, a2 and a3 respectively.

(A29)
We note that for given parameters of the main part of the smoothed functions (for example, ǫ, r 0 , d and r c for the V Morse part of f 1 ), the parameters of the smoothed cutoff regions (b low , b high , r c,low , r c,high for f 1 ) are uniquely determined by ensuring differentiability of the function at the boundaries (r low and r high for f 1 ) and thus they are not explicitly listed in the function arguments and are not provided in the tables of values of model parameters.