David H. Mathews

Umbrella sampling is widely used to estimate biomolecular free energy landscapes and relative folding stabilities. Although equilibration is a critical component of umbrella sampling workflows, the impact of equilibration length on thermodynamic predictions remains poorly understood. Here, we investigate the effect of equilibration length on relative folding free energy predictions for four RNA hairpins with loop sequences GUGAAA, CUGGGA, GUAAUA, and UUAAUU with helical stems of three base pairs. Umbrella sampling simulations were performed using an end-to-end distance reaction coordinate spanning 15-45 Å, where equilibrium simulations (windows) were spaced at roughly 1 Å intervals. In these calculations, the hairpin stem-loops were allowed to equilibrate in an end-to-end distance window and then the coordinates were transferred to the next larger end-to-end distance window to equilibrate. Two equilibration lengths, 2 ns and 100 ns per window, were followed by 600 ns of production sampling. Potential of mean force (PMF) profiles were reconstructed using the Weighted Histogram Analysis Method (WHAM) and used to calculate pairwise free energy differences with thermodynamic cycles. Increasing the equilibration length produced substantial, sequence-dependent changes in the reconstructed free energy landscapes. The 100 ns protocol generated markedly flatter PMFs for GUGAAA and UUAAUU and pronounced reshaping of the free energy landscape for GUAAUA. These changes were accompanied by reductions in hydrogen-bonding and stacking interactions, particularly within the intermediate regions of the reaction coordinate. The resulting thermodynamic predictions, as free energy change differences, were therefore highly sensitive to equilibration length. Across nearly all hairpin pairs, the 100 ns equilibration yielded substantially larger magnitude free energy change difference values than the corresponding 2 ns equilibration, with differences that greatly exceeded replica-to-replica variability. Comparison with optical melting measurements and nearest-neighbor thermodynamic predictions revealed that 2 ns equilibration times more closely agreed with experimental values than those obtained using 100 ns equilibration. These findings demonstrate that longer equilibration can systematically alter the structural ensembles sampled during umbrella sampling and amplify predicted stability differences without improving agreement with experiment. More widely, our results highlighted equilibration length as a critical and nontrivial parameter in RNA free energy calculations and demonstrate that increased equilibration does not necessarily lead to more accurate thermodynamic predictions.

RNA structure is a key determinant of RNA function and regulation. The coupling of chemical probing technologies, such as SHAPE, with deep sequencing has enabled large-scale experimental characterization of RNA structures in complex samples and under diverse conditions. Furthermore, probing data are often used to guide thermodynamics-based secondary structure prediction algorithms and have been shown to improve their accuracy. However, current algorithms treat these single-nucleotide measurements as statistically independent signals, inherently overlooking short-range dependencies in the data. Here, we show that discretized SHAPE data display context dependence within loop regions and within stem regions and we use Markov models to formally capture such dependencies. We then leverage Markov modeling in the classification of small structure motifs from their discretized SHAPE data signatures and subsequently integrate the classifying feature into the dynamic programming recursions that underlie computational RNA folding. Compared to state-of-the-art SHAPE-guided structure prediction methods, our Markov-informed framework improves prediction performance. Furthermore, we identify SHAPE signatures characteristic of highly stable hairpins, such as GAAA, GCAA, and UUCG tetraloops, and integrate these insights into the folding recursions to further improve predictions. Overall, the proposed framework provides a foundation for context-aware statistical modeling of SHAPE data, particularly in loop regions, where signal characterization has proven challenging due to high variance. This work further demonstrates that finer modeling of SHAPE data has the potential to push the limits of data-guided secondary structure prediction.

RNAs regulate gene expression and cellular processes, often relying on specific conformations for function. RNA folding is hierarchical and sequence-dependent, with nearest-neighbor thermodynamic models commonly used to predict secondary structure. Current models were developed using optical melting experiments in 1 M NaCl, which does not represent the cellular environment. To address this, we developed a new model in Advanced Dulbecco’s Modified Eagle Medium (Adv. DMEM), which mimics mammalian extracellular ionic composition. This in vivo -like model provides RNA folding parameters for helical base stacks and loop motifs. Optical melting experiments revealed helical stacks, particularly tandem G-U pairs, are less stabilizing in Adv. DMEM. Loop parameters were generally destabilizing but highly dependent on both sequence and loop type, with internal loops displaying idiosyncratic behavior. Structure prediction benchmarking revealed minimal differences overall, except for tRNAs, which showed improved prediction reliability and enhanced cloverleaf stability. Notably, tRNAs lack internal loops, suggesting further studies in Adv. DMEM could refine secondary structure predictions. This in vivo -like parameter set is included in the RNAstructure software package. Grounding these parameters in a physiologically relevant environment, we improve the biological relevance of RNA secondary structure predictions and establish a foundation for studying RNA folding under in vivo conditions. Graphical

RNA design, also known as RNA inverse folding, aims to find RNA sequences that fold into a target secondary structure. However, recent work has shown that some target structures are provably undesignable, where no RNA sequence can fold into it as the minimum free energy (MFE) structure. In this paper, we go beyond this binary, MFE-based designability and explore a soft, probability-based designability that upperbounds the Boltzmann probability of any design and quantifies how easily or likely any design might possibly fold into the target structure. We introduce a theory of ensemble approximation and a probability decomposition framework for bounding the folding probabilities of RNA structures and motifs in an explainable way. We further develop a linear-time dynamic programming algorithm that efficiently searches over exponentially many decompositions. Combining ensemble approximation with dynamic decomposition search, our method efficiently identifies the optimal motif decomposition that yields the tightest probabilistic bound for a given structure. Our framework is applicable to any factorizable energy model or scoring function that decomposes onto loops. Applying our work, LinearDecompose, to both native and artificial RNA structures in the ArchiveII and Eterna100 datasets, we obtained much tighter probability bounds than baselines. Our work also provides anatomical tools for analyzing RNA structures and pinpointing the sources of design difficulty at the motif level. Source code and data are available at https://github.com/shanry/RNA-Undesign.

During translation initiation, the 40S small ribosomal subunit is recruited to the mRNA 5′ cap and scans the 5′ untranslated region (UTR) to locate the start codon. While the mechanism of 40S translocation remains elusive, the RNA helicase eIF4A has long been suspected as the primary molecular motor driving 40S scanning. In this study, we utilized GFP reporter mRNAs to investigate the impact of 5′ UTR length on translational efficiency. We found that an 8-fold variation in the length of unstructured 5′ UTRs did not lead to substantial changes in translation efficiency in wheat germ extract (WGE) and human HEK293T cell lysate. By contrast, the presence of a stable stem-loop in the middle of the 5′ UTR significantly reduced cap-dependent translation. These results suggest that mRNA scanning is not rate-limiting when the UTR is devoid of secondary structure. Inhibition of eIF4A by hippuristanol in cell-free protein synthesis systems yielded an equivalent decrease in translation for mRNAs with short and long unstructured 5′ UTRs, indicating that eIF4A may be dispensable for 40S scanning. Our data suggest that helicase-independent one-dimensional diffusion may be the primary mechanism enabling 40S movement along the 5′ UTR during initiation.

All-atom simulations of RNA using molecular dynamics have the promise of modeling conformational preferences, folding thermodynamics, conformational change kinetics, and binding affinities of small molecule therapeutics. These simulations rely on a force field, a set of equations and parameters that model the potential energy as a function of conformation using classical mechanics. One popular force field for RNA is Amber OL3, with the most recent iteration derived in 1999 and with subsequent updates to backbone dihedral parameters. The Amber force field, while frequently used, is known to have limitations; for example, it does not properly stabilize native structures against alternative structures. Here, we provide a new approach to fitting the non-bonded parameters for the force field, specifically atom-centered point charges for electrostatics and the Lennard-Jones parameters. The parameters are fit to quantum mechanics (QM) interaction energies calculated with symmetry-adapted perturbation theory (SAPT), including embedded point charges to represent the electrostatic field from solvent and adjacent nucleotides. In this pilot study with a limited set of fitting data, we use the Amber ff99 equations and atom types unchanged. With the revised parameters, we observe improvement in the stability of native structures relative to alternative structures. Native tetraloop conformations, which unfold with the Amber OL3 force field, are stable on the microsecond timescale with our new force field parameters. We also see improvement in the conformational preferences of tetramers. Crucially, A-form helices are still well-modeled, but we observe additional flexibility in an internal loop that is not consistent with NMR data. Overall, we provide evidence that this new approach to fitting RNA force field parameters to SAPT interaction energies with native-structure context represented as embedded point charges is promising. It offers a flexible solution for revising the equations in future work or for extension to other molecules that interact with RNA, such as proteins and small molecules. We call this new set of force field parameters Amber RNA.ROC26.

Nearest neighbor parameters are widely used in software for estimating the conformational stability of an RNA sequence folding into a specific structure. Folding stability for RNA with canonical nucleotides A, C, G, and U has been widely studied, but the same is not true for most modified nucleotides. In this work, we present a comprehensive set of nearest neighbor parameters for estimating the folding stability of RNA including pseudouridine in helical or loop contexts. These parameters are derived from 210 optical melting experiments involving helices with pseudouridine-A and pseudouridine-G pairs and with pseudouridine in loop motifs. The experiments include sequences with pseudouridine and U in the same strand, including U-A and U-G pairs, allowing us to consider the folding stability of sequences with both U and pseudouridine. On average, pseudouridine stabilizes RNA folding compared to U in an analogous motif, although this effect is sequence-context dependent. These parameters improve the modeling of folding stability for RNA secondary structures containing pseudouridine. We demonstrate that these parameters successfully model the secondary structure change for Saccharomyces cerevisiae U2 snRNA when two additional inducible pseudouridines are present. These parameters are freely available and incorporated into the RNAstructure software package.

Nearest neighbor analysis is commonly used to estimate RNA folding stabilities. In this contribution, we report a set of RNA folding nearest neighbor parameters for estimating free energy change for RNA sequences including 1-methyl-pseudouridine. Development of mRNA vaccines has identified 1-methyl-pseudouridine as a key nucleobase modification for suppressing innate immune responses. However, the contributions of these modifications to RNA folding stability were unclear. Our new parameters provide helical terms for 1-methyl-pseudouridine-adenine and 1-methyl-pseudouridine-guanine base pairs. The parameters also estimate loop stabilities for loops with 1-methyl-pseudouridine or a combination of 1-methyl-pseudouridine and uridine. These parameters are derived using 208 optical melting experiments and tested against an additional 16 optical melting experiments. On average, we find that substitution of uridine with 1-methyl-pseudouridine stabilizes RNA folding, with the extent of stabilization depending on adjacent sequence. The estimation of tRNA folding ensembles for tRNA sequences with 1-methyl-pseudouridine was significantly improved using the new nearest neighbor parameters. The new nearest neighbor parameters are provided as part of the RNAstructure software package. With these parameters, the secondary structures of natural sequences with 1-methyl-pseudouridine and mRNA therapeutics fully substituted with 1-methyl-pseudouridine can be modeled.

RNA design aims to find a sequence that can fold into a target secondary structure. It can create artificial RNA molecules for specific functions, with wide applications in medicine. It is computationally challenging due to two levels of combinatorial explosion: the exponentially large design space and the exponentially many competing structures per design. Popular methods such as local search cannot keep up with these combinatorial explosions. We instead employ two techniques from machine learning, continuous optimization and Monte-Carlo sampling. We start from a distribution over all valid sequences, and use gradient descent to improve the expectation of an arbitrary objective function. We define novel coupled-variable distributions to model the correlation between nucleotides. We then use sampling to approximate the objective, estimate the gradient, and select the final candidate. Our work consistently outperforms state-of-the-art methods in key metrics including Boltzmann probability and ensemble defect, especially on long and hard-to-design structures.