# Simulation-Guided Conformational Space Exploration to Assess Reactive Conformations of a Ribozyme

Sélène Forget, Marie Juillé, Elise Duboué-Dijon, Guillaume Stirnemann

► To cite this version:

Sélène Forget, Marie Juillé, Elise Duboué-Dijon, Guillaume Stirnemann. Simulation-Guided Conformational Space Exploration to Assess Reactive Conformations of a Ribozyme. Journal of Chemical Theory and Computation, 2024, 20 (14), pp.6263-6277. 10.1021/acs.jctc.4c00294. hal-04737774

HAL Id: hal-04737774

https://hal.science/hal-04737774v1

Submitted on 15 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Simulation-Guided Conformational Space
Exploration to Assess Reactive Conformations
of a Ribozyme

Sélène Forget,*,† Marie Juillé,†,‡ Elise Duboué-Dijon,‡ and Guillaume Stirnemann*,†

†PASTEUR, Département de chimie, École Normale Supérieure, PSL University, Sorbonne Université, CNRS, 24 rue Lhomond, 75005 Paris, France

‡Université Paris Cité, CNRS, Laboratoire de Biochimie Théorique, 13 rue Pierre et Marie Curie, 75005, Paris, France

E-mail: selene.forget@ens.psl.eu; guillaume.stirnemann@ens.psl.eu

Abstract

Self-splicing ribozymes are small RNA enzymes that catalyze their own cleavage
through a transphosphoesterification reaction. While this process is involved in some
specific steps of viral RNA replication and splicing, it is also of importance in the
context of the (putative) first autocatalytic RNA-based systems that could have pre-
ceded the emergence of modern life. The uncatalyzed phosphoesterbond formation is
thermodynamically very unfavorable, and many experimental studies have focused on
the understanding of the molecular features of catalysis in these ribozymes. However,
chemical reaction paths are short-lived and not easily characterized by experimental
approaches, so that molecular simulation approaches appear as an ideal tool to unveil
the molecular details of the reaction. Here, we focus on the model hairpin ribozyme.
We show that identifying a relevant initial conformation for reactivity studies can be
highly challenging due to limitations in both available X-ray and the force field accu-
racy, together with the necessity of extensive sampling. This is frequently overlooked
in mixed quantum/classical studies that predominantly concentrate on the chemical
reaction itself. These challenges stem from limitations in both available experimental
structures (which are chemically altered to prevent self-cleavage) and in the accuracy
of force fields, together with the necessity for comprehensive sampling. We show that
molecular dynamics simulations, combined with extensive conformational phase space
exploration with Hamiltonian replica-exchange simulations, enable to characterize the
relevant conformational basins of the minimal hairpin ribozyme in the ligated state,
prior to self-cleavage. We find that what is usually considered as a canonical reac-
tive conformation with active site geometries and hydrogen-bond patterns that are
optimal for the addition-elimination reaction with general acid/general base catalysis,
is metastable and only marginally populated. The thermodynamically stable confor-
mation appears to be consistent with the expectations of a mechanism that does not
require the direct participation of ribozyme residues to the reaction. While these ob-
servations may suffer from forcefield inaccuracies, all investigated forcefields lead the
same conclusions upon proper sampling, contrasting with previous investigations
on shorter timescales suggesting that at least one reparametrization of the Amber99
forcefield allowed to stabilized aligned active site conformations. Our study demon-
strates that identifying the most pertinent reactant state conformation holds equal
importance alongside the accurate determination of the thermodynamics and kinetics
of the chemical steps of the reaction.

Introduction

Ribonucleic acid (RNA) is one of the important and essential players of the life machinery.

While the original, textbook picture was that RNA serves as an intermediate in the production of proteins from DNA (a role today ascribed to so-called messenger RNAs), many more non-coding functions have been evidenced over the years.¹⁻⁵ Remarkably, it was discovered
in the 80s that some RNAs, called ribozymes, are capable of catalysis,⁶,⁷ a function which
is traditionally that of a special class of proteins (the enzymes) in living systems. However,
the chemistry enabled by ribozymes is more limited as compared to protein enzymes.⁸⁻¹¹
Typically, they catalyze either the ligation of two RNA fragments, or their cleavage (the op-
posite reaction). While the former is the key reaction for the formation of the (putative) first
autocatalytic systems based on RNA,⁸,¹² the latter remains a relevant mechanism in current
organisms, as it is involved in some specific steps of viral RNA replication and splicing.¹³

The chemical reaction catalyzed by these ribozymes (either in its forward or backward
direction) is a transphosphoesterification between a 2'OH hydroxyl group and a 3'-5' phospho-
ester bond. ¹⁴ In the absence of catalysts, this reaction exhibits an extremely high free-energy
barrier in the 30-40 kcal/mol range, ¹⁵,¹⁶ precluding the efficient polymerization of RNA in
abiotic conditions. ¹⁷ How do ribozymes succeed in catalyzing these reactions to such an
extent that they become kinetically accessible? In modern enzymes catalyzing the same re-
action, divalent cations are often involved, but in ribozymes, it is still unclear whether these
have a direct catalytic effect or whether they act as structural scaffolds. ¹⁴,¹⁸ In fact, some
ribozymes, including the hairpin ribozyme that will be the focus of the current work, do not
require divalent cations to function. ¹⁹ This poses the question of the molecular mechanism
by which a few catalytic residues of the ribozyme allow to significantly lower the reaction
free-energy barrier.

Just like for protein enzymes, experimental studies bring invaluable insights onto the
molecular details of catalysis.²⁰ For example, available structures at several steps along the
catalytic cycle, complemented by transition state or intermediate mimics, can offer impor-
tant clues about the corresponding molecular mechanisms.²¹ Mutagenesis is instrumental in
characterizing the essential residues in the active site or beyond.²² But the chemical reaction
itself is an ultrafast process that is not easily tracked in real time in the experiments. Calcu-
lation and simulation approaches, in particular molecular dynamics simulations, can prove
useful and bring relevant molecular details on the reactions, complementing, rationalizing,
and/or predicting experimental observations.

Mixed quantum-classical (QM/MM) approaches allow for the study of reactivity within the framework of atomistic simulations.²³,²⁴ If we take the example of the hairpin ribozyme, a lot of effort and attention have been devoted to the modelling of the QM part, the accuracy of the employed quantum approach, or the enhanced sampling approach employed to probe the reaction pathway.²³,²⁵⁻³¹ At the same time, an extensive exploration of the conformational phase space in the reactive state is not routinely addressed. Although this might not be an issue for many systems, for which the structure of the reactive conformation is unambiguous, ribozymes pose significant challenges, for at least three reasons. First, being typically involved in self-splicing reactions, chemical modifications of the nucleophilic 2'OH are required for structural studies of the ribozyme in a ligated state.²¹,²⁸,³² Consequently, modification-free catalytic site structures are not experimentally accessible. Second, RNAs are known to be more flexible than proteins, and a given sequence sometimes even present completely different native conformations depending on the conditions.³² Finally, while molecular dynamics can be helpful in characterizing the relevant conformations, results are heavily dependent of molecular parameters (force field) that have been considerably improved for RNA in the last years, but probably do not reach the level of reliability of those for proteins yet. There is no universal force field that would be optimal for all RNA systems.³³

Here, we focus on the hairpin (HpR) ribozyme, one of the most extensively studied small self-cleaving ribozymes,¹⁹ which shares a common catalytic scaffold with other ribozymes.³⁴ Despite many experimental and computational studies, no clear consensus has been reached so far about the self-cleavage mechanism.³⁵,³⁶ From the simulation perspective, we were initially intrigued by the report of distinct reaction pathways with very different reactant conformations,³⁶⁻⁴⁰ of important force field effects,²⁸ and by the variety of strategies that are followed to generate initial reactive conformations starting from the crystal (methylated structure) (see Figure 1).³⁶⁻⁴⁰ We pose apparently simple but crucial questions: what is/are the relevant active site conformation(s) for self-splicing? How different are they from the
crystal structure? Do several state-of-the-art force fields lead to the same results? Do
several protocols for the equilibration and the preparation of the reactive state conformation
eventually converge to the same, minimum free-energy state?

<figure><img src="image_2.png" /><figcaption>Figure 1: Different views of the minimal hairpin ribozyme (HpR) studied in this work. (A) Sequence and secondary structure of the minimal HpR, with the four helices H1 to H4, the two loops A and B, and the cleavage site indicated with a pink arrow in loop A. Residues in bold correspond are conserved among the hairpin ribozyme family.<sup>19</sup> (B) Crystallographic structure (PDB 2OUE) with the same nomenclature as in panel (A), and (C) zoom on the active site geometry of the O2'-methylated structure (other missing hydrogens were added following the procedure detailed in the Methods section).</figcaption></figure>

We perform extensive exploration of the conformational phase space using molecular
dynamics complemented with enhanced sampling techniques, validated across several force
fields, which has been systematically overlooked in previous simulation work on the system.
We show that the so-called "reactive" active site conformation, that is inspired by the features
of the HpR mimic transition state (TS) structure<sup>41</sup> and that is required with a general
acid/base mechanism where the proton transfers are assisted by neighboring nucleotides,<sup>35,36</sup>
appears as metastable. This conformation, that is very different from the crystal structure of
the reactant state, is usually generated by restraints<sup>28,30</sup> guided by chemical intuition or by
analogy with artificial mimics of transition state structures.<sup>41</sup> We show that it evolves on long
(microsecond and above) timescales toward another conformation that can serve as an initial
stage of a mechanism involving a phosphate-assisted proton transfer through a non-bridging
oxygen, which was not captured in previous simulation work restricted to shorter sampling.
Our work illustrates the difficulties in determining what is a "good" starting conformation
for reactivity studies involving ribozymes, and stresses the crucial importance of defining a
reliable strategy to characterize the conformational ensemble prior to any direct study of the
chemical step.

Materials and methods

Simulation and analysis softwares

All simulations, and some of the system preparation, were done using Gromacs 2019.4.42 patched with PLUMED 2.7.43,44 We used the VMD software45 combined with the the Collective Variables Dashboard46 for visualization, and use scripts benefiting from the python open source MDAnalysis 47,48 and Barnaba 49 libraries to analyze the data. Input files required to reproduce these simulations are available on online Zenodo repository, as referenced at the end of the article.

System preparation

In most of our study, we started from the precatalytic state structure of the minimal hairpin
ribozyme, which has been resolved by crystallography at 2.05 Å at pH 6.0 (PDB ID 2OUE).⁵⁰
As mentioned in the introduction, the precatalytic state must be neutralized by the presence
of a methyl group on the nucleophile O2' oxygen atom. The group was manually removed
from the coordinate file. Structures were then prepared using Gromacs 2019.4 pdb2gmx
tool.⁴² Following previous studies,²⁸,³⁰ A38 was protonated at the N1 position and G8 was
left neutral. Other bases were taken in their standard protonation state at neutral pH.
Crystallographic anions (SO₄²⁻ and NCO) were removed; crystallographic water molecules
were preserved. The system was solvated in a 0.2M KCl solution using Gromacs tools, either
with TIP3P or OPC water model, adding up to 32000 solvent molecules. In total, 156 K⁺ and 120 Cl⁻ were added in the bulk by replacing water molecules.

In order to study the effect of the initial ribozyme conformation on the subsequent exploration during the MD simulation, we also obtained from the authors of Ref.²⁸ a trajectory in the so-called "reactive" active site conformation from which we extracted one starting frame $L_1^{0,pucker}$, that was then resolvated following the procedure indicated above. We finally followed a strategy inspired by Veenis et al.,³⁰ where the aligned reactive conformation was generated using restraints on some active site hydrogen-bonds (H-bonds) based on chemical intuition and to mimic the expected transition state structure (see more details below).

## Force fields

The four RNA force fields tested in this study are all derivatives of the AMBER ff99 force field,⁵¹ corrected by the Barcelona $\alpha/\gamma$⁵² and the Olomouc $\chi_{OL3}$ reparametrizations.⁵³ We refer to this version as $\chi_{OL3}$. We also tested the addition of the $\epsilon\zeta_{OL1}$ corrections,⁵⁴ which were suggested to be essential to maintain the active site in a reactive, aligned conformation.²⁸ Finally, we also tried the gHBfix21 additional constraints from Frohlking et al.⁵⁵ meant to improve the description of large diversity of RNA motifs. For the protonated Adenine (RAP) residues, we used the parameters derived from Mlynsky & al.²⁷ after updating the atom types (the employed topology is given in the Supplementary Information).

We tested two water models, OPC and TIP3P, together with the Amber14\_ff99bsc0$\chi_{OL3}\epsilon\zeta_{OL1}$ and Amber14\_ff99bsc0$\chi_{OL3}$ force fields (OPC and TIP3P), while $\chi_{OL3}$-gHBfix21 was used with OPC (see Table 1). K⁺ and Cl⁻ were described with Joung and Cheatham ion parameters.⁵⁶

We also tried the latest version from the DES-Amber force field family,⁵⁷⁻⁵⁹ DES-Amber 3.20.⁵⁹ These force fields are all reparametrizations of the Amber14\_ff99bsc0$\chi_{OL3}$ force field, with improvements validated against experimental data. Released in 2022,⁵⁹ the DES-Amber 3.20 force field includes all the features of its predecessor DES-Amber SF1.0 version.⁵⁷ These
include the refinement of the van der Waals parameters and the charges of nitrogens and
polar hydrogens on the bases of the residues. Additionally it includes the correction of the
global charge of nucleobases from 1.0 to 0.9, re-distribution of the charges in between O2',
C2' and HO2' atoms, modification of certain Lennard-Jones parameters and the fine tuning
of α, β, γ, δ, ε and ζ dihedrals. DES-Amber 3.20 Gromacs files were obtained directly from
DESRES. These parameter files did not include the RAP modified residue, which was built
by combining the DES-Amber and Mlynsky et al.²⁷ charge modifications (see Supplementary
Information). DES-Amber simulations were run with TIP4P water model.⁶⁰

Relaxation and equilibration protocol

Each system is minimized prior to the MD simulation as follows. In most cases, restraints of
1000 kJ/mol/nm² are applied to all backbone atoms during the first three steps of the pro-
cedure: 5000 minimization steps followed by consecutive 500 ps NVT and NPT equilibration
at T = 298.15 K and P = 1 bar. Temperature was controlled using a V-rescale tempera-
ture coupling algorithm⁶¹ and pressure using the Parrinello-Rahman coupling scheme.⁶²,⁶³
The restraints are progressively released during three consecutive relaxations (with restraint
forces of 1000, 100, 10 kJ/mol/nm²) in order to let the nucleic backbone relax. Finally, the
entire system is freely simulated over 100 ps under NVT conditions, and then for 100 ps
under NPT conditions. The detailed equilibration protocol, along with the corresponding
parameter files of each equilibration step, can be found online on the online repository zen-
odo.11033971. The electrostatic interactions are computed using the particle-mesh Ewald
method,⁶⁴ with a Verlet cutoff scheme of 10.0 Å employed for Lennard-Jones interactions.
The LINCS algorithm⁶⁵ is applied to all bonds containing hydrogen atoms, and the simula-
tion timestep was set to 2 ps in all cases. The grid dimensions are set at 0.12 nm with cubic
interpolation.

The above protocol was applied directly starting when from the crystallographic struc-
ture. We also performed simulations starting from two conformations for which the active
site geometry had been previously constrained to adopt an "aligned" geometry. First, we kindly obtained from the authors of Ref. 28 a trajectory that was obtained as follows: they restrained the C1'-C2'-C3'-C4' dihedral of the A-1 ribose for 500 ns of simulation using Amber14 ff99bsc0χOL3 parametrization. Within this restrained simulation, they picked a structure which matched the L₁ basin characteristics described below, that was later propagated (and stable) for 500 ns with Amber14 ff99bsc0χOL3εζOL1. We took the first frame of their trajectory for subsequent equilibration and production, and we call this starting configuration L₁⁰,pucker.

We tried another strategy adapted the procedure of Veenis & al. 30 by constraining 4 distances (G8:H1···A-1:O2', G8:H2···G+1:pro-Sp, A38:H1···G+1:O5', A38:H6···G+1:pro-Rp, see Figure 2) using a harmonic potential :

$$U_i^{\text{harmonic}}(x_i) = k_i(x_i - a_i)^2 \qquad (1)$$

with $a_i = 2$ Å and $k_i = 2000$ kJ/mol/nm² for all 4 $x_i$ distances. For these specific runs, the restraints were maintained during minimization and equilibration for 1 µs. The equilibration procedure was also adapted to be closer to the spirit of Veenis & al. 30 (inputs available on zenodo.11033971).

## Production phase

We later propagated MD trajectories for several hundreds of nanoseconds (typically 500 ns to 1 µs) in the NPT ensemble using the same simulation parameters as described above. For each system, simulations were replicated thrice. Another set of simulations were performed using a Hamiltonian replica exchange scheme, here Replica Exchange with Solute Tempering (REST2). 66,67 We used 24 replicas rescaling factors λ ranging from 1 to 0.667, attempting exchanges every 1 ps. The average exchange rate is close to 15%. Simulations were propagated for 200 ns to 1 µs, and the analysis of the unscaled replica was done after an equilibration
<figure><img src="image_3.png" /><figcaption>Figure 2: View of the 4 H-bonds constrained following Veenis & al.:³⁰ A-1:O2'···G8:N1, G+1:pro-Sp···G8:N2, G+1:O5'···A38:N1 and G+1:pro-Rp···A38:N6 H-bonds.</figcaption></figure>

period of 100 ns to 200 ns. The convergence time was estimated using block analysis of the catalytic site conformations population ratios and the eRMSD⁶⁸ time evolution of the unscaled replica (see Figure S1.). A summary of all propagated trajectories of this study is given in Table 1.

## Descriptors for the conformational basins

Three main conformational basins were explored in our simulations as well as in previous works.²⁵⁻³⁰ In order to classify the obtained structures, we used three descriptors: the orientation of the phosphate of G+1 residue with respect to the nucleophilic O2', the puckering pseudo-rotation angle of the A-1 residue, and the H-bond involving O2' as acceptor. Figure 3 summarizes the combination required for a structure to be labeled as belonging to the crystallographic basin ($L_C$), the conformation usually referred to as "reactive" by Mlynsky et al.,²⁸ which we call the $L_1$ basin, and another conformation mainly characterized by different phosphate orientation and A-1 sugar puckering, $L_2$. The criteria for $L_C$ were chosen following
<table><caption>Table 1: Overview of the different simulations performed for this study, with the detailed type of sampling, the parametrisation, the initial structure and the simulations length. Here $\chi_{OL3}\epsilon\zeta_{OL1}$ stands for Amber14 ff99bsc0$\chi_{OL3}\epsilon\zeta_{OL1}$ and $\chi_{OL3}$ for Amber14 ff99bsc0$\chi_{OL3}$. $L_1^{0,pucker}$ was extracted from Mlynsky et al.²⁸ trajectory (which was generated following a protocol which involved restraints on the pucker of A-1 residue, see Methods). $L_1^{0,HB}$ was generated for the purpose of this study, following Veenis et al³⁰ procedure in which four HBs are restrained (see Methods).</caption><thead><tr><th>Sampling</th><th>Force Field</th><th>Water Type</th><th>Initial Structure</th><th>Simulation [Converged] Times</th></tr></thead><tbody><tr><td>MD</td><td>$\chi_{OL3}\epsilon\zeta_{OL1}$</td><td>TIP3P</td><td>PDB 2OUE</td><td>500 ns (×3)</td></tr><tr><td>REST2</td><td>$\chi_{OL3}\epsilon\zeta_{OL1}$</td><td>TIP3P</td><td>PDB 2OUE</td><td>500 [100-500] ns</td></tr><tr><td>MD</td><td>$\chi_{OL3}$</td><td>TIP3P</td><td>PDB 2OUE</td><td>500 ns (×3)</td></tr><tr><td>MD</td><td>DES-Amber 3.20</td><td>TIP4P</td><td>PDB 2OUE</td><td>500 ns (×3)</td></tr><tr><td>MD</td><td>$\chi_{OL3}\epsilon\zeta_{OL1}$</td><td>TIP3P</td><td>$L_1^{0,pucker}$</td><td>1200 ns (×1), 700 ns (×2)</td></tr><tr><td>MD</td><td>$\chi_{OL3}$</td><td>TIP3P</td><td>$L_1^{0,pucker}$</td><td>500 ns (×3)</td></tr><tr><td>MD</td><td>$\chi_{OL3}\epsilon\zeta_{OL1}$</td><td>OPC</td><td>$L_1^{0,pucker}$</td><td>500 ns (×1), 700 ns (×2)</td></tr><tr><td>MD</td><td>$\chi_{OL3}$</td><td>OPC</td><td>$L_1^{0,pucker}$</td><td>700 ns (×1), 500 ns (×2)</td></tr><tr><td>MD</td><td>$\chi_{OL3}$ + gHBfix21</td><td>OPC</td><td>$L_1^{0,pucker}$</td><td>500 ns (×3)</td></tr><tr><td>MD</td><td>DES-Amber 3.20</td><td>TIP4P</td><td>$L_1^{0,pucker}$</td><td>500 ns (×3)</td></tr><tr><td>REST2</td><td>$\chi_{OL3}\epsilon\zeta_{OL1}$</td><td>TIP3P</td><td>$L_1^{0,pucker}$</td><td>200 [100-200] ns</td></tr><tr><td>REST2</td><td>$\chi_{OL3}$</td><td>OPC</td><td>$L_1^{0,pucker}$</td><td>200 [100-200] ns</td></tr><tr><td>REST2</td><td>$\chi_{OL3}$ + gHBfix21</td><td>OPC</td><td>$L_1^{0,pucker}$</td><td>200 [100-200] ns</td></tr><tr><td>REST2</td><td>DES-Amber 3.20</td><td>TIP4P</td><td>$L_1^{0,pucker}$</td><td>500 [200-500] ns</td></tr><tr><td>MD</td><td>$\chi_{OL3}\epsilon\zeta_{OL1}$</td><td>TIP3P</td><td>$L_1^{0,HB}$</td><td>750 ns (×1), 500 ns (×2)</td></tr><tr><td>REST2</td><td>$\chi_{OL3}\epsilon\zeta_{OL1}$</td><td>TIP3P</td><td>$L_1^{0,HB}$</td><td>1 µs</td></tr></tbody></table>
the characteristics of the minimal HpR's precatalytic state crystallographic structure (PDB: 2OUE).⁵⁰ Those for L₁ were chosen following Mlynsky and co-coworker's description of this basin.²⁸ The requirements for L₂ were chosen following the conformations explored in our own simulations of the hairpin ribozyme, and very similar to those seen in Ref. ²⁹

<figure><img src="image_4.png" /><figcaption>Figure 3: Criteria for the definition of the catalytic site conformations used in this study: view of the distances (in pink for the A-1:O2' H-bonds; in orange for the orientation of the phosphate) and the puckering pseudo-rotation angle of A-1 5-ring sugar cycle (in turquoise) which are collectively utilized to ascertain the conformation of the catalytic site. See Figure 5C-E for snapshots of the L~C~, L~1~ and L~2~ arrangements.</figcaption></figure>

Current picture on the self-cleavage mechanism

Before we discuss the results of our own investigation, we present more structural details
about the minimal hairpin ribozyme, and discuss the results of previous investigations re-
garding its mechanism. The minimal HpR is composed of two independent domains,¹⁹,⁶⁹
each of them containing 2 helices separated by a loop (helices H1 and H2, loop A; helices H3
and H4, loop B respectively) (see Figure 1). The catalytic site is assembled when the two
loops A and B come into contact to form a precatalytic conformation. The cleaved phos-
phoester bond belongs to the catalytic strand, at the intersection of A-1 and G+1, which
are surrounded by G8 (loop A) and A38 (loop B) residues (see Figure 1A). In the crystallographic structure of the HpR's precatalytic state⁵⁰ (Figure 1 B and C), G8 base faces the A-1 residue, possibly forming 2 HBs. A-1:O2' faces G8:N1, and the non-bridging phosphate oxygen pro-Sp (O1P) faces the G8:N2 atom. On the other side, A38:N1 is aligned with the A-1:O3' atom and is close to A-1:O5'. The inline attack angle (IAA), defined as the angle between the O2' nucleophile, phosphorus, and the O5' leaving group, which is required to be as flat as possible to enable the addition-elimination reaction (Figure 1C and Figure 4), is at 155° within the crystal, and the A-1:O2'...G+1:P distance ($d_{O2'–P}$) is 2.93 Å.

It has been shown that A38 and G8 residues are important for the reaction: the HpR reactivity decreases by a significant 14,000-fold reduction when A38 is mutated,⁷⁰ and a notable 350-fold reduction upon the substitution of G8;⁷¹⁻⁷³ G8:N1 seems to have a critical role for the reaction;⁷⁴,⁷⁵ combined with the measured pH dependence of the reaction rate, this led to the idea of a general acid-base mechanism involving these two residues.³⁵,⁷²,⁷³,⁷⁶⁻⁷⁸ The protonation of A38 is supported by a number of results, including an estimation of its pKa⁷⁹ and the comparison between simulations with A38 and A38H⁺,²⁷ supporting prior structural analysis.⁵⁰ The protonation state of G8, with a pKa close to 9.5,⁸⁰ remains debated. In particular, MD simulations with deprotonated G8 resulted in a very unstable active site.²⁷,²⁸ While this does not rule out a role of G8 as a base, this may suggest G8 to be only transiently deprotonated during the reaction.

In the general acid-base mechanism (or dianionic scenario, Figure 4A), the A-1:O2' is deprotonated by G8:N1⁻ (general base) and attacks the phosphorus through an addition—elimination reaction. On the opposite side of the phosphate, A38H⁺ acts as general acid for the protonation of the leaving O5'. An alternative scenario is the inline mono-anionic mechanism⁴⁰ (Figure 4B), in which O2' is deprotonated by one of the non-bridging phosphorus oxygens. Considering which one among pro-Rp (O2P) and pro-Sp (O1P) would deprotonate A-1:O2' likely depends on which of these accept an H-bond from A-1:O2', and thus on the phosphate orientation in the preacatalytic state. This mechanism doesn't directly involve G8
which is believed to remain into its canonical protonation state at biological pH (G8:N1H). A38H⁺ can be either involved as proton shuttle, or acting as general acid.<sup>36-38,40</sup> These chemical pathways have been characterized by semi-empirical Quantum Mechanics/Molecular Mechanics (QM/MM) investigations,<sup>37,38</sup> followed by single-point QM/MM studies.<sup>36,39</sup> The computed energy barrier values are overall in agreement with experiments (20-21 kcal/mol).<sup>19</sup> Notably, it has been suggested that these scenarios may operate concurrently, given their closely matched energy barrier values.<sup>36,39</sup>

A Dianionic scenario or general base-acid mechanism

<figure><img src="image_5.png" /></figure>

B Inline mono-anionic scenario

<figure><img src="image_6.png" /></figure>

Figure 4: Different scenarios of the self-cleavage mechanism of the HpR. (A) The dianionic scenario or general base-acid mechanism. (B) The inline mono-anionic scenario. We chose not to represent A38 to not take a definitive stance on its role, although it is most probably involved in the reprotonation of G+1:O5' terminal end. Similarly, the nature of the non-bridging oxygen that is involved in the substrate-assisted proton transfer (pro-Sp or pro-Rp) is not specified to represent simultaneously both mono-anionic mechanistic pathways.

The pH dependence of the reaction is often put forward as the strongest experimental argument in favor of a dianionic scenario.<sup>35,72,73,76-78</sup> However, one cannot rule out a structural role of the protonation state of these bases, which need not be directly involved in the catalytic mechanism. Since the two mechanisms correspond to different precatalytic active
site structures, it is therefore crucial to assess the explored active site conformations that
would be prone to one mechanism or another.

Results

Exploration starting from the unmethylated PDB structure

We first prepared an initial conformation using the crystallographic structure of the minimal
HpR in the precatalytic state (2OUE),⁵⁰ with A-1:O2' unmethylated, G8 in its canonical
state and A38 protonated on its N1 atom.

We initially chose the Amber14\_ff99bsc0χ<sub>OL3</sub>ϵζ<sub>OL1</sub> force field in TIP3P water, following previous observations that the ϵζ<sub>OL1</sub> corrections were necessary in order to properly reproduce the active site geometry,²⁸ with favorable (flat) in-line attack angle (IAA) between the nucleophilic A-1:O2' oxygen, the phosphorus and the leaving G+1:O5' oxygen. More details about the system preparation are given in the Methods sections.

A first production run, replicated thrice, was propagated for 500 ns. While the system starts from a conformational basin around the crystal structure conformation (that we will refer to as the L<sub>C</sub> basin in what follows), after a few dozens of ns, a major conformational change occurs. It is characterized by a sudden flip of the A-1 puckering pseudo-rotation angle (see Figure 5A, in blue) from C2'-endo (as in the C basin) to C3'-endo (as most frequently observed in RNAs). Based on the characteristics of these new conformations, we define a new conformational basin L<sub>2</sub>.

This reorientation seems at first surprising, because the ϵζ<sub>OL1</sub> corrections have been ex-
plicitly shown<sup>28</sup> to avoid this often-considered artifactual conformation (which is for example
seen with the Amber14\_ff99bsc0χ<sub>OL3</sub> force field), for which the IAA is not as flat as in a so-
called "reactive" conformation<sup>28</sup> (that we refer to as L<sub>1</sub> here, see Methods). Actually, most
simulations of the minimal HpR, with the same sequence, protonation states (A38H<sup>+</sup>/G8)
and in the same environment have also led to the L<sub>2</sub> basin<sup>25-27,29,81</sup> when using other force
<figure><img src="image_7.png" /><figcaption>Figure 5: Brute force MD simulation of HpR. (A) Time evolution of key collective variables reporting on the catalytic site conformation starting from the crystallographic structure, for one of the replicas. (B) Same data for a simulation initiated from $L_1^{0,pucker}$, an artificially-generated structure belonging to the $L_1$ basin (see text). Both A and B runs were parameterized using the Amber14 ff99bsc0$\chi_{OL3}\in\zeta_{OL1}$ RNA force field in TIP3P water. (C) Typical structure of the active site in the $L_C$ basin; (D) in the $L_1$ basin and (E) in the $L_2$ basin.</figcaption></figure>

fields. We will come back later on the differences between our results and that of previous work using the same force field, and will show that this originates from the selection of the initial conformation for the simulation, and to previously limited sampling.²⁸

Typical conformations of the catalytic site in the $L_C$, $L_1$, and $L_2$ basins are represented in Figure 5 C-E. Based on these observations and on the literature,²⁸ we chose five variables that allow to define three criteria to distinguish between the three basins (the H-bond environment of O2', the phosphate orientation, and the A-1 sugar puckering), as described in the Methods section.

In order to assess the rearrangement of the catalytic site, we monitor the time evolution
of 9 variables along the trajectory, among which the IAA angle and the 5 discriminating variables (Figure 5A). The puckering pseudo-rotation angle flip (in blue) coincides with the rearrangement of the H-bond network in the L<sub>C</sub> conformation active site, and to the appearance of the L<sub>2</sub> conformation. This transition also coincides with the H-bond breaking between A-1:O2' and G8:N2 to form new ones with A38:N6 and A10:N6. The phosphate rotates with A-1:O2', maintaining the A-1:O2-pro-Rp HB, and the pro-Sp oxygen ends up facing G8:N1 with which it forms a new H-bond after a few dozens of ns.

Figure 5A highlights the decrease in IAA values when reaching the L₂ basin (average of 121° in the C basin and 108° in the L₂ one). This is clearly visible in the 2D distributions along the IAA and dO₂'P (Figure 6A), suggesting that the puckering conformation indeed has a direct impact on these coordinates and shifts them toward values that are a priori less favorable to reaction, as previously suggested.²⁸

<figure><img src="image_8.png" /><figcaption>Figure 6: 2D distributions along the IAA and the dO2'P distance (see main text). The entire trajectories of one replica per system contribute to these data, with colors corresponding to the different conformational basins visited during these simulations (see Methods for conformational labeling criteria). (A) Simulation run with Amber14 ff99bsc0χOL3εζOL1 force field and initiated from the crystallographic structure and (B) from the L<sub>1</sub><sup>0,pucker</sup> basin.</figcaption></figure>

The transition from L<sub>C</sub> to L<sub>2</sub> is a stochastic event and the limited timescale of our simulations does not guarantee that the system would not further go back to the L<sub>C</sub> basin or evolve towards another conformational basin (such as L<sub>1</sub>) at longer times. In order to better assess
the conformational space exploration starting from the crystallographic structure, we con-
ducted enhanced sampling simulations using Replica Exchange Solute Tempering (REST2),
using 24 replicas for a simulation time of 500 ns (see Methods). The goal of this method
is to accelerate the conformational changes of the ribozyme using a ladder of simulations
that differ in the rescaling of the molecular interactions involving ribozyme atoms, and it
was extensively used to enhance the conformational space exploration of biomolecules,⁸²⁻⁹²
including RNA.⁹³⁻⁹⁸
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: View of the catalytic site behavior during a REST2 simulation initiated from HpR's crystallographic structure (100-500 ns converged frames). (A) Classification of the observed conformations based on the values of 9 collective variables (y axis), organized from most populated (left) to less populated (right) clusters. The width of each cluster is proportional to its weight, indicated on the x axis for the most populated clusters. Conformations with less than 0.01% weight are not shown but are aggregated to evaluate the combined value of these rare arrangements (here, 39.5%). (B) 2D distribution along the IAA and $d_{O2'-P}$, with a color code corresponding to the explored conformational basin (see Methods for conformational labeling criteria).</figcaption></figure>

Similarly to the brute force MD simulation, we observed a rapid change in the puckering of the A-1 sugar from C2'-endo to C3'-endo within a few nanoseconds (because of the enhanced sampling strategy, the transition occurs much faster as compared to the brute force MD).
Figure 7A displays the most populated combination of our set of 9 variables during the converged part of the simulation (100-500 ns) that belongs to the L₂ conformational basin. As observed in our unbiased simulations, the 2D distribution along the IAA angle and the d<sub>O2'–P</sub> distance shows two subpopulations (see Figure 7), with one peak around an IAA of ≈ 120° and d<sub>O2'–P</sub> ≈ 3.3 Å, and another one around ≈ 105° and d<sub>O2'–P</sub> ≈ 3.8 Å.

High values (IAA> 140°) are very rarely explored. The REST2 simulation additionally samples a very small fraction of conformations with very low IAA values, which were not seen in the brute force MD, as well as conformations that do not explore any of the L<sub>C</sub>, L<sub>1</sub> or L<sub>2</sub> basins (3%). Within this simulation framework (in particular, with the chosen force fields), this suggests that L<sub>2</sub> is the most stable conformation of the hairpin ribozyme.

At this stage, we do not find any clear evidence of the exploration of the L₁ basin using the Amber14\_ff99bsc0χ<sub>OL3</sub>εζ<sub>OL1</sub> RNA force field when we start from the PDB structure. To later confirm that these observations are force field independent, we repeated our replicated brute force MD simulations with two other force field combinations.

With the Amber14\_ff99bsc0χ<sub>OL3</sub> force field, the puckering flip and the transition out of the crystallographic basin L<sub>C</sub> are almost immediate (less than 3 ns). During the rest of the 500 ns simulation, the system explores configurations belonging to the L<sub>2</sub> basin (see Figure S2A). Interestingly, the fraction of structures with very low IAA values (60°), which was almost absent with χ<sub>OL3</sub>εζ<sub>OL1</sub> parametrization, is more important, as seen before.<sup>28</sup> When using the DES-Amber 3.20 force field, the active site is much more flexible. The G8 strand moves far from the catalytic strand and the HBs it is involved in (G8:N1...A-1:O2', G8:N2...G+1:pro-Sp, G8:N2...G+1:pro-Rp) are broken (Figure S2B). While most of the structures still belong to the L<sub>2</sub> basin, they thus exhibit different H-bond patterns that locally differ from those explored with Amber14\_ff99bsc0χ<sub>OL3</sub> and Amber14\_ff99bsc0χ<sub>OL3</sub>εζ<sub>OL1</sub>. We will come back later to the peculiar behavior of the DES-Amber force field, that will be better characterized with enhanced sampling. The important conclusion however is that none of the three force fields tested here yield structures in which the catalytic site exhibits
some of the important features of the L₁ basin.

Our simulations thus confirm that when starting from the PDB structure, "reactive" conformations (L₁ basin) prone to the dianionic general acid/base or the pro-Sp monoan-ionic mechanisms, are never observed, regardless of the employed force field. Since previous reports of such structures were based on simulations initially constrained to adopt an aligned conformation,²⁸,³⁰ one could conclude that the transition from the L₂ to the L₁ basin is simply out-of-reach using our simulation framework, because of its large free-energy barrier. In order to study the stability of the L₁ basin, we now turn to simulations that were initially driven into this basin.

Exploration starting from L₁ structures

The L₁ basin is very different from the crystallographic structure. Switching from one to another involves the rotation of the A-1 phosphate (see Figure 5 C and D), which required to break the G8:N2...G+1:pro-Rp H-bond and to form the G8:N2H...G+1:pro-Sp and A38:N6...G+1:pro-Rp HBs. To the best of our knowledge, this conformational change and the formation of the L₁ structure from the crystallographic structure have been achieved using two different strategies. Mlynsky et al. reached the L₁ conformational basin by equilibrating the HpR structure while constraining the A-1 puckering into the C2'-endo conformation.²⁸ Upon our request, they kindly shared their trajectory with us so that we could start from a similar, aligned structure, that we denote L₁<sup>0,pucker</sup>.

On the other hand, Veenis et al. chose to constrain the H-bond network.³⁰ They set up an equilibration protocol in which G8:H1...A-1:O2', G8:H2...G+1:pro-Sp, A38:H1...G+1:O5', A38:H6...G+1:pro-Rp are maintained below 2 Å by a harmonic potential (see Figure 2). One should notice that these distances are present in the crystallographic structure of the TS vanadate mimic⁴¹ but not in the precatalytic state.⁵⁰ We reproduced this procedure and generated 1 µs of biased MD simulation (see Methods). We then released the constraints starting from the last frame (L₁<sup>0,HB</sup>) and simulated the system with no bias for 750 ns in one
replica, and 500 ns in two others. These biased structures differed from the ones provided
by Mlynsky et al. by an average RMSD of 2.78 Å, but also belonged to the L₁ basin.

Starting from $L_1^{0,pucker}$, we propagated one replica for 1.2 µs of force brute MD with the
Amber14 ff99bsc0$\chi_{OL3}\epsilon_{OL1}$ force field (see Figure 5B), and two additional ones for 700 ns.
For several hundreds of nanoseconds, explored conformations fall into the $L_1$ basin, with
A-1:$O2'...G8:N1$ H-bond formed, the A-1 residue puckering pseudo rotation in C2'-endo and
G+1:pro-Sp facing A-1:$O2'$. In the first replica, after about 493 ns, the system overcomes
a sharp transition from the $L_1$ to the $L_2$ basin, which is associated with a change in A-
1 puckering, the rupture of some of the key HBs characteristic of the $L_1$ basin, and the
formation of a new H-bond pattern typical of the $L_2$ basin (Figure 5B). No transition was
observed in the 700 ns of the other two replicas.

The IAA and $d_{O2'_{-P}}$ distributions also change with the transition (Figure 6B) and reach
values very much similar to those observed starting from the crystallographic structure (Fig-
ure 5A and Figure 7B). During the ≈700 ns following the transition, no backward transition
from $L_2$ to $L_1$ nor further rearrangement from the $L_2$ are observed.

Starting from $L_1^{0,HB}$, the system stays in the $L_1$ basin for the duration of the simulation
(Figure S3) in two of the three replicas of the MD trajectory. In one of them, we observe
a sharp increase in distance between A38 and the rest of the active site around 370 ns
(Figure S3B), but this does result in a change in the explored values for $d_{O2'_{-P}}$ and IAA
(Figure S3C). The third replica undergoes a rapid C2'-endo/C3'endo puckering flip.

We chose not to propagate this simulation and rather turned to our enhanced sampling
strategy in order to better characterize the conformational landscape both when starting
from $L_1^{0,HB}$ and $L_1^{0,pucker}$, as now discussed.

In the REST2 simulations starting from $L_1^{0,pucker}$, the transition to $L_2$ unsurprisingly
occurs much faster (20 ns) as compared to the brute force MD. After the steady-state is
reached, the resulting geometries of the active site (Figure S5A) are very similar to those
obtained starting from the crystallographic structure (Fig 7). While a significant population
of L₂ conformations appears around 200 ns, the REST2 simulations starting from L₁⁰,HB predominantly stay in the L₁ basin for much longer (800 ns) (Figure 8). After that time, about one third of conformations belong to the L₂ basin, while the majority of the other conformations do not belong to any of the well-identified basins L<sub>C</sub>, L<sub>1</sub> or L<sub>2</sub>. Further investigation reveals that once a physical replica (which periodically exchanges its rescaling parameter) transitions into the L₂ basin, it never rearranges into the L₁ conformation. A given physical replica thus either remains in the L₁ basin, irreversibly (from the perspective of these simulations) transitions into L₂, or explores other arrangements not well-characterized within the framework of this study. Setting aside the identification of the hidden variable blocking the transition, these observations lead to the conclusion that the L₁ conformation could be classified as metastable, whereas the L₂ conformation is stable.

<figure><img src="image_11.png" /><figcaption>Figure 8: Evolution of the proportion of L<sub>C</sub>, L<sub>1</sub> and L<sub>2</sub> populations along time in the lowest λ = 0 replica trajectory of a REST2 simulation propagated from L<sub>1</sub>⁰,HB, using Amber14 ff99bsc0χ<sub>OL3</sub>εζ<sub>OL1</sub>, suggesting that the initial conformation equilibrates on very long time scales but slowly departs from the L<sub>1</sub> basin</figcaption></figure>

The different behavior of simulations originating from L₁⁰,pucker and L₁⁰,HB, despite the initial structures being close on average, suggests that they differ by a subtle conformational change to which our definition of the L₁ basin is insensitive, and that these are not equivalent ways to generate an initially-aligned conformation. While the identification of the precise differences between these two structures goes beyond the scope of the current work, it would deserve future investigation.
We repeated our brute force and enhanced sampling simulations starting from $L_1^{0,pucker}$ using the variety of RNA force fields (see Table 1 and Method section for our choices of parametrizations). We also varied the water model with both Amber14\_ff99bsc0$\chi_{OL3}$ and Amber14\_ff99bsc0$\chi_{OL3}\epsilon\zeta_{OL1}$ and OPC. We first performed brute force MD to monitor the time scale of the $L_1$ to $L_2$ transition (Figure S4). Each force field combinations are propagated three times, from three random initial velocities. In all 6 combinations, we observed a puckering flip of the A-1 sugar and a departure from the $L_1$ basin within 500 ns for at least one the three MD replicas. Transition times are summarized in Table 2.

<table><caption>Table 2: Pucker reorientation times for brute force MD simulations initiated from $L_1^{0,pucker}$, with different force field parametrizations. These simulations were replicated thrice, leading to the three different values mentioned for each system.</caption><thead><tr><td>Force field</td><td>Water Model</td><td>Transition Time (ns)</td></tr></thead><tbody><tr><td>Amber14\_ff99bsc0$\chi_{OL3}\epsilon\zeta_{OL1}$</td><td>TIP3P</td><td>493, >700, >700</td></tr><tr><td>Amber14\_ff99bsc0$\chi_{OL3}$</td><td>TIP3P</td><td>25, 54, 150</td></tr><tr><td>Amber14\_ff99bsc0$\chi_{OL3}\epsilon\zeta_{OL1}$</td><td>OPC</td><td>20, 696, >700</td></tr><tr><td>Amber14\_ff99bsc0$\chi_{OL3}$</td><td>OPC</td><td>155, 254, 498</td></tr><tr><td>Amber14\_ff99bsc0$\chi_{OL3}\epsilon\zeta_{OL1}$ + gHBfix21</td><td>OPC</td><td>228, 295, >500</td></tr><tr><td>DES-Amber 3.20</td><td>TIP4P</td><td>16, 87, 130</td></tr></tbody></table>

To further enhance the conformational space exploration, we propagated REST2 simulations with three different force field combinations (Figure S5 B-D). In all cases, and as already observed with our original force field choice, the $L_1$ basin is not populated anymore after at most a few tens of nanoseconds, leading to the more stable $L_2$ basin in all cases.

With enough sampling, all four combinations (in addition to what was readily observed) irreversibly depart from this basin. These include the RNA Amber14\_ff99bsc0$\chi_{OL3}\epsilon\zeta_{OL1}$ forcefield, contrasting with prior claims that it allowed to stabilize the $L_1$ basin.²⁸ Upon longer sampling, we show that this is always a metastable state. These multiple tests show that the metastability of the $L_1$ is consistent when varying the water model, the dihedral terms, the atomic charges, the H-bond strength and/or the van der Waals parameters.

We pause to discuss some important differences between the DES-Amber reparametrization and the other force fields. When using Amber14\_ff99bsc0$\chi_{OL3}$, Amber14\_ff99bsc0$\chi_{OL3}\epsilon\zeta_{OL1}$
or Amber14 ff99bsc0χ<sub>OL3</sub>+gHBfix21, we recover the typical distributions of IAA and A-1:O2'-G+1:P distances identified before (Figure S5 A-C). However, the simulation employing DES-Amber exhibits a peculiar behavior following the transition from the L<sub>1</sub> to the L<sub>2</sub> basin (Figure S5D), which was readily observed when starting from the crystallographic structure (Figure S2B). This trajectory is characterized by wide distributions of IAAs and A-1:O2'...G+1:P distances (Figure S5D, right panel). This suggests an increase flexibility of the active site region, also apparent from the large fluctuations of the overall HpR structure, with multiple eRMSD peaks (Figure 9A), significant RMSD values (Figure 9B) and high solvent accessible areas (Figure 9C). In particular, the surrounding G8 residue is pulled out from the catalytic site, releasing its respective HBs with A-1 and G+1 (Figure S5D, left panel), which allows for the phosphate rotation. When such an event occurs, G8 cannot assure its structural role and the number of accessible active site conformations increases. Among these fluctuations, aligned and tight conformations a priori prone to an addition reaction appear (Figure S5D, right panel). However, it is crucial to notice here that these transient arrangements are not accompanied by an H-bond pattern that would correspond to the expectations of a dianionic general acid/base mechanism, especially because G8 is not close enough anymore. The A-1 sugar pucker is consistently in the C3'-endo conformation. These rather surprising observations may suggest that the DES-Amber force field may not be the most adequate choice to study the HpR behavior.

Discussion

Here, we characterize the relevant conformational basins of the minimal hairpin ribozyme in
the ligated state, prior to self-cleavage, using an enhanced sampling trajectories shedding new
light on this mostly-overlooked aspect of ribozyme reactivity. We classified the structures
in three main conformational basins: the crystallographic basin L_c ; a first reactive basin
L_1 that is a priori favorable to either a dianionic scenario (Figure 4A) or a monoanionic
<figure><img src="image_12.png" /><figcaption>Figure 9: Increased flexibility with the DES-Amber force fields. (A) eRMSD<sup>68</sup> distribution with respect to the common initial structure L<sub>1</sub><sup>0,pucker</sup></figcaption></figure>

(B) RMSD distributions with respect to the common initial structure L<sub>1</sub><sup>0,pucker</sup> and (C) total Solvent Accessible Surface Area (panel B) of two converged REST2 simulations with two different parametrizations: Amber14 ff99bsc0χ<sub>OL3</sub>εζ<sub>OL1</sub> (TIP3P) and DES-Amber 3.20 (TIP4P).

scenario (Figure 4B) involving phosphate-assisted proton exchange through pro-Sp<sup>36-39</sup> ;
and a second one, L<sub>2</sub>, that is more favorable to a monoanionic scenario involving phosphate-
assisted proton exchange through pro-Rp (Figure 4B).<sup>36-38,40</sup> Starting from the L<sub>c</sub> basin,
all simulations quickly evolve to and remain stable in the L<sub>2</sub> basin. Other simulations
starting from conformations within L<sub>1</sub>, which are obtained by applying restraints to drive
the active site into this conformation, either based on chemical intuition or based on the
expectations on a transition state structure mimic,<sup>28,30</sup> also evolve toward the L<sub>2</sub> basin,
sometimes on long microsecond timescales. The reverse transition is never observed. The L<sub>1</sub>
basin is thus metastable. The same results are obtained for a variety of force fields, including
Amber14 ff99bsc0χ<sub>OL3</sub>εζ<sub>OL1</sub> which had been claimed to stabilize the L<sub>1</sub> basin as compared
to L<sub>2</sub>,<sup>28</sup> and very recent latest version reparametrization of the DES-Amber family.<sup>59</sup> findings
of this study raise a number of questions that we now discuss.

Are there still force field limitations?

Our results may seem at first surprising. Most simulation studies consider L₁, and not L₂, as
the reactive conformation in the ligated state.²⁸,³⁰ L₁ allows for inline attack angles (IAA)
and A-1:O2' distances ($d_{O2'_{-P}}$) favorable to the O2' addition on the phosphorus, while it is not directly the case in the $L_2$ conformation. In fact, $L_2$, which coincides with the flip of the A-1 ribose puckering pseudo-rotation angle and a major change in the local hydrogen-bonds pattern, was considered as an artifact due to the limitations of most force fields that was supposedly solved by the $\epsilon\zeta_{OL1}$ reparametrization.²⁸ However, we show here that using 5 different RNA force fields, as well as several water models, $L_1$ is always metastable with respect to $L_2$ upon an enhanced exploration of the conformational phase space. No RNA forcefield was found to stabilize $L_1$ over $L_2$ upon proper sampling, contrasting with previous work.

We cannot rule out that current RNA force fields, of which we have tested some of the most recent instances, all fail at reproducing the active site structure of the hairpin ribozyme. In particular, it is quite striking that the ribose of the nucleophilic O2' is most often found in an unusual C2' endo puckering conformation in the precatalytic structure of several self-cleaving ribozymes,⁹⁹⁻¹⁰¹ when our systems here undergo a inevitable flip to C3' endo conformation. There was also some interesting report on model systems with blocked sugar conformations suggesting that the formation of a cyclic 2',3' phosphate was more favorable in the south (C2'-endo) conformation.⁹⁹ This would suggest that the A-1 puckering is actually an important feature of hairpin catalysis that is not correctly reproduced by current force fields.

While the consistency among very significantly different forcefields strengthen these conclusions, further validation is not straightforward. Despite numerous trials, the slow variable for the conformational transitions remain unknown. For example, it is not the puckering angle alone, as we failed to converge extensive, hundreds-of-nanoseconds free-energy calculations along this coordinate. It is highly unlikely that QM/MM calculations would converge on much shorter timescale, would the reaction coordinate be known. Moreover, such calculations would not capture the effects of any forcefield inaccuracy in the MM region, which is actually likely to be the main driver for the conformational transition. Experimental
confirmation is also challenging, as the self-reactive state of the ribozyme cannot be isolated.

We pause to note that the DES-Amber force field<sup>59</sup> we have tested gives surprising and maybe unrealistic flexibility of the active site region, and more generally, results in large fluctuations of the ribozyme structure. It has been shown that the DES-Amber reparametrizations were accurate and successful to reproduce the structural and thermodynamical features of a number of small RNAs,<sup>57,102–104</sup> with still some deficiencies in the description of stacking and non-canonical interactions.<sup>33,57,105</sup> Here, the behavior of the hairpin ribozyme with this set of parameters raises questions and would deserve more in-depth investigation. more generally it would also be interesting to study how other ribozymes behave with this force field, which has never been done.

## What is a reliable starting structure for MD?

As discussed above, we cannot rule out an important deficiency in all tested force fields. Adding to this, we have readily pointed to the major problem that is faced when studying self-cleaving ribozymes: as compared to many other systems, for which a structure of the actual reactive state can be experimentally obtained without chemical modifications of the biomolecule itself, ribozymes need to be chemically altered (typically by removing the O2' or by methylating it) to prevent self degradation of the structure.<sup>20,50</sup> We show that these chemical modifications have important impacts on the active site structure. Upon removal of the methyl group, the $L_C$ basin becomes unstable, in particular because of different steric hindrance and the changes in the H-bond pattern. There is therefore a lack of direct experimental structures which we can compare to our simulation results.

Importantly, the crystallographic structure is very different from the $L_1$ reactive basin. To the best of our knowledge, this reactive conformation has never been obtained spontaneously from the crystal structure of the precatalytic state without additional restraints on the active site geometry, for example on the ribose conformation and by human-guided selection along a restrained trajectory.<sup>28</sup> Another strategy has been to restrain the H-bond network<sup>30</sup> following
chemical intuition or by similarities with the structure of a TS mimic of the active site, but
that is at the same time very different from the H-bond network that can be inferred from
the precatalytic state crystal structure. However, we have shown here that L₁ conformations
prepared in such a way eventually relax into the L₂ basin, sometimes on timescales well
above the microsecond. This thus raises questions about such restraints: while they seem to
follow the expectations of at least the general dianionic acid/general base mechanism, and
assuming that the force fields can be trusted, they do not correspond to the most stable
reactive state conformation, which is obtained much more rapidly by relaxing the crystal
structure in solution upon removal of the O2' methyl group.

We finally note that hairpin might be a relatively simple example of ribozyme as it
does not appear to have a cation-dependent mechanism.¹⁰⁶ Identifying the correct reactive
structure starting from the crystallographic data may be more complex for other ribozymes,
because it can be influenced by the metal ion environment.¹⁰⁷,¹⁰⁸

Insights into the reaction mechanism

Previous comparative studies based on QM and QM/MM calculations suggested that the
monoanionic scenario (with phosphate-assisted proton transfers) and the dianionic one (with
general acid/general base catalysis) had similar free-energy barriers.³⁶ However, it is not
clear from which conformational basin these calculations were performed, as the mechanistic
study was initiated from intermediate structures (phosphoranes for the monoanionic scenario,
and deprotonated G8 residue) equilibrated with force fields which very likely lead to L₂
conformations.

Let us first assume that the free-energy barriers for the two mechanisms are close. As the
L₁ to L₂ conformational transition occurs on a hundreds of nanoseconds to microsecond time
scale (and even more when starting from the L₁<sup>0,HB</sup> structure), a quick transition state theory
estimation of the barrier would lead to ≈ 10 kcal/mol, which is probably overestimated if we
add the effect of a prefactor proportional to a diffusion coordinate along the conformational
transition coordinate. Considering that we never observed the reverse transition, even using
enhanced sampling simulation, we can reasonably assume that the barrier for the L₂ to L₁
transition, at least within the framework of the employed force field, exceeds 10–12 kcal/mol.
As a consequence, if the reaction free-energy barrier (independent of this conformational
change) is indeed similar for both mechanisms, it is thus very unlikely that the system
would first change its conformation from L₂ to L₁ and then follow a dianionic or a pro-Sp
monoanionic scenario. In that case, the most favored mechanism would be a the pro-Rp
monoanionic one.

Indeed, if the L₂ basin is a priori not in a conformation allowing for flat IAA and short
dO2'–P distances that are required for the addition—elimination mechanism, previous simula-
tion work<sup>26,29</sup> showed that the intermediate that would be formed by a first proton transfer
between the O2' and the phosphate pro-Rp non-bridging oxygen (which are H-bonded in
the L₂ conformation) leads to inline arrangements that are very much compatible with a
nucleophilic addition of the deprotonated O2' on the phosphorus. Building on these results,
metadynamics-QM/MM simulations were performed to evaluate the free energy landscape
of the inline mono-anionic scenario.<sup>40</sup> Their mechanism shows that G8 is important to stabi-
lize the intermediate following the first proton-transfer intermediate, as suggested before,<sup>26</sup>
and that A38H<sup>+</sup> could be involved in the proton reshuffling following the O5' departure
and reprotonation. Therefore, the shift towards lower IAAs in the L₂ conformational basin
as compared to the L₁ basin is not necessarily an issue and does not mean that L₂ is an
artefact. In fact, it adopts an ideal geometry for a first proton transfer to the pro-Rp oxy-
gen. However, we notice that in this previous studies,<sup>26,40</sup> the simulations were initiated
from a structure showing Lc features and that a proper QM/MM study starting from L₂
conformation remains to be done.

In this framework, and as discussed before,⁴⁰ the key role played by G8 and A38H⁺ in
terms of electronic stabilization could explain the pH-dependence of the catalytic rate³⁵,⁷²,⁷³,⁷⁶⁻⁷⁸
and the sensitivity of the reaction towards mutations of these residues.²⁷,²⁸ They need not
be directly involved in catalysis. It is interesting to note that the monoanionic scenario, with
phosphate-assisted proton transfers, is actually reminiscent of the uncatalyzed mechanism in
bulk water, that was recently studied by some of us using deep-neural network inter-atomic
potentials in particular,¹⁰⁹,¹¹⁰ and by others before.¹¹¹,¹¹² The role played by the ribozyme
active site residues here would be to provide an electrostatic environment minimizing the
free-energy barriers of the reaction,¹¹³ but this reaction would proceed though a very similar
pathway as compared to the bulk.

Finally, we cannot rule out that the dianionic mechanism would be kinetically more
favorable than the monoanionic mechanism. In that case, one could imagine that starting
from the free-energy minimum of the L₂ basin, a conformational change should first trigger
the formation of the metastable L₁ conformation before the reaction could proceed through a
dianionic mechanism. The apparent reaction free-energy barrier would contain a contribution
both from the chemical step, but also from the conformational change. Such an effect has
been evidenced before, for example to explain the temperature-dependence of an enzymatic
reaction rate constant that was due to a conformational change affecting the active site
rather than an intrinsic temperature-dependence of the chemical step.⁸⁹,¹¹⁴

Conclusions

In this work, we studied the conformations of the minimal hairpin ribozyme in the ligated
precatalytic state for self-splicing using molecular dynamics simulations combined with en-
hanced sampling approaches. The crystal structure⁵⁰ provides an important starting point
for the simulations, but it is chemically altered and might not be representative of reactive
conformations in solution. Upon demethylation of the nucleophilic O2' oxygen, and proper
relaxation and equilibration of the simulations, we find that the system preferentially adopts
a conformation ideal for a monoanionic scenario with phosphate-assisted proton transfer
through the pro-Rp oxygen, while the often sought-after conformation ideal for a general
acid/general base dianionic scenario is not observed.

Instead of starting directly from the crystal structures, restraints were used in some
previous works<sup>28,30</sup> to generate conformations in this second "aligned" basin. We stress that
from the perspective of the active site geometry, the resulting basin is very different from
the crystallographic structure. In any case, we demonstrate that upon further equilibration
of simulations starting from the aligned basin, all spontaneously further evolve into the
other reactive basin on timescales sometimes exceeding the microsecond. This represents an
important novelty over previous work that did not explicitly addressed these long timescales
and may have been stuck in a metastable basin.

We tested five different RNA force fields,<sup>53-55,57,59</sup> in combination with 3 water models. All parameter sets give the same qualitative answer. Our results imply that what is usually considered<sup>28,30</sup> as a canonical reactive conformation with active site geometries and H-bond patterns that are optimal for the addition—elimination reaction, is metastable and is only very marginally populated. Instead, the most stable conformation corresponds to an H-bond pattern that is in fact ideally formed for a phosphate-assisted proton transfer. This conformation is not incompatible with the experimental observations that G8 and A38H<sup>+</sup> are crucial for the reaction,<sup>35,72,73,76-78</sup> because they could play important roles in terms of electronic stabilization that could in particular explain the pH-dependence of the catalytic rate. Even if the relative populations of the different conformations could be affected by force field inaccuracies, our results strongly suggest that two conformations may coexist in the active site, each prone to favor a different mechanism.

Another possibility that we cannot exclude at this point would be that the dianionic
mechanism would be kinetically more favorable than the monoanionic mechanism. In such
a scenario, the initial step of the reaction would entail a conformational transition before
proceeding further, that we hope to be able to quantify in the near future. It would also be
important to characterize the reaction mechanisms starting from these two basins within the
same QM/MM simulation framework, and eventually compare the two pathways. Lastly, we
emphasize once more that despite thorough validation and consensus among all examined force fields, our conclusions may still be subject to artefacts inherent to all these force fields. An essential supplementary examination would involve evaluating the stability of the crystal structure through simulations of the (O2' methylated) crystal, as its behavior might diverge significantly from that observed in simulations in solution.

Our study highlights the challenges associated with identifying a "suitable" initial conformation for reactivity investigations involving ribozymes. These challenges stem from limitations in both available experimental structures and the accuracy of force fields, compounded by the necessity for comprehensive sampling even more critically than for the (usually) less flexible enzyme active sites. Our findings underscore the critical need for meticulous consideration of this aspect, which is frequently overlooked in mixed quantum/classical studies that predominantly concentrate on the chemical reaction itself. Building upon some of our recent work in enzymatic reactions,⁸⁹,¹¹⁴ it becomes evident that accurately identifying the most pertinent reactant state conformations is imperative for understanding the apparent reaction rate constant. This aspect holds equal importance alongside the accurate determination of the thermodynamics and kinetics of the chemical step.

## Funding

The research leading to these results has received funding from the European Research Council under the European Union's Eighth Framework Program (H2020/2014-2020)/ERC Grant Agreement No. 757111 (GS). This work was also supported by the "Initiative d'Excellence" program from the French State (Grant "DYNAMO", ANR-11-LABX-0011-01 to GS and EDD). MJ was funded by a PhD fellowship from the OCAV program of the Paris Science et Lettres University.
## Data availability

Additional details about simulation parameters, simulation convergence, results with additional force fields can be found in the Supporting Information. All input data required for running the simulations discussed here can be found online on Zenodo (zenodo.11033971).

## Acknowledgments

The simulations presented here benefited from a local computing platform administered by G. Letessier, and benefited from the HPC resources of TGCC under the allocation A0130811005 made by GENCI (Grand Equipement National de Calcul Intensif). We thank Petr Jurecka and Marie Zgarbova for sharing their trajectory in the $L_1$ conformational basins and for very useful discussions.

## References

(1) Elbashir, S. M.; Harborth, J.; Lendeckel, W.; Yalcin, A.; Weber, K.; Tuschl, T. Duplexes of 21-nucleotide RNAs mediate RNA interference in cultured mammalian cells. *Nature* **2001**, *411*, 494–498.

(2) Steitz, T. A. A structural understanding of the dynamic ribosome machine. *Nat. Rev. Mol. Cell Bio.* **2008**, *9*, 242–253.

(3) Jinek, M.; Chylinski, K.; Fonfara, I.; Hauer, M.; Doudna, J. A.; Charpentier, E. A programmable dual-RNA-guided DNA endonuclease in adaptive bacterial immunity. *Science* **2012**, *337*, 816–821.

(4) Cech, T. R.; Steitz, J. A. The noncoding RNA revolution - Trashing old rules to forge new ones. *Cell* **2014**, *157*, 77–94.
(5) Setten, R. L.; Rossi, J. J.; ping Han, S. The current state and future directions of RNAi-based therapeutics. *Nat. Rev. Drug Disc.* **2019**, *18*, 421–446.

(6) Cech, T. R.; Zaug, A. J.; Grabowski, P. J. In vitro splicing of the ribosomal RNA precursor of tetrahymena: Involvement of a guanosine nucleotide in the excision of the intervening sequence. *Cell* **1981**, *27*, 487–496.

(7) Kruger, K.; Grabowski, P. J.; Zaug, A. J.; Sands, J.; Gottschling, D. E.; Cech, T. R. Self-splicing RNA: Autoexcision and autocyclization of the ribosomal RNA intervening sequence of tetrahymena. *Cell* **1982**, *31*, 147–157.

(8) Doudna, J. A.; Cech, T. R. The chemical repertoire of natural ribozymes. *Nature* **2002**, *418*, 222–228.

(9) Lilley, D. M. The origins of RNA catalysis in ribozymes. *Trends Biochem. Sci.* **2003**, *28*, 495–501.

(10) Fedor, M. J.; Williamson, J. R. The catalytic diversity of RNAs. *Nat. Rev. Mol. Cell Bio.* **2005**, *6*, 399–412.

(11) Jimenez, R. M.; Polanco, J. A.; Lupták, A. Chemistry and Biology of Self-Cleaving Ribozymes. *Trends Biochem. Sci.* **2015**, *40*, 648–661.

(12) Higgs, P. G.; Lehman, N. The RNA World: molecular cooperation at the origins of life. *Nat. Rev. Genet.* **2014**, *16*, 7–17.

(13) Webb, C.-H. T.; Riccitelli, N. J.; Ruminski, D. J.; Lupták, A. Widespread Occurrence of Self-Cleaving Ribozymes. *Science* **2009**, *326*, 953–953.

(14) Lilley, D. M. J. Mechanisms of RNA catalysis. *Philos. Trans. R. Soc. B-Biol. Sci.* **2011**, *366*, 2910–2917.
(15) Lad, C.; Williams, N. H.; Wolfenden, R. The rate of hydrolysis of phosphomo-
noester dianions and the exceptional catalytic proficiencies of protein and inositol
phosphatases. *Proc. Natl. Acad. Sci.* **2003**, *100*, 5607–5610.

(16) Schroeder, G. K.; Lad, C.; Wyman, P.; Williams, N. H.; Wolfenden, R. The time required for water attack at the phosphorus atom of simple phosphodiesters and of DNA. *Proc. Natl. Acad. Sci.* **2006**, *103*, 4052–4055.

(17) Szostak, J. W. The eightfold path to non-enzymatic RNA replication. *J. Syst. Chem.* **2012**, *3*, 2.

(18) Hanna, R.; Doudna, J. A. Metal ions in ribozyme folding and catalysis. *Curr. Opin. Chem. Biol.* **2000**, *4*, 166–170.

(19) Fedor, M. J. Structure and function of the hairpin ribozyme. *J. Mol. Biol.* **2000**, *297*, 269–291.

(20) Fedor, M. J. Comparative Enzymology and Structural Biology of RNA Self-Cleavage.
*Annu. Rev. Biophys.* **2009**, *38*, 271–299.

(21) Spitale, R. C.; Volpini, R.; Mungillo, M. V.; Krucinska, J.; Cristalli, G.; Wedekind, J. E. Single-Atom Imino Substitutions at A9 and A10 Reveal Distinct Effects on the Fold and Function of the Hairpin Ribozyme Catalytic Core. **2009**, *48*, 7777–7779.

(22) Roberts, J. M.; Beck, J. D.; Pollock, T. B.; Bendixsen, D. P.; Hayden, E. J. RNA sequence to structure analysis from comprehensive pairwise mutagenesis of multiple self-cleaving ribozymes. *eLife* **2023**, *12*, e80360.

(23) Banáš, P.; Jurečka, P.; Walter, N. G.; Šponer, J.; Otyepka, M. Theoretical studies of RNA catalysis: Hybrid QM/MM methods and their comparison with MD and QM. *Methods* **2009**, *49*, 202–216.
(24) Van Der Kamp, M. W.; Mulholland, A. J. Combined quantum mechanics/molecular mechanics (QM/MM) methods in computational enzymology. *Biochemistry* **2013**, *52*, 2708–2728.

(25) Ditzler, M. A.; Sponer, J.; Walter, N. G. Molecular dynamics suggest multifunctionality of an adenine imino group in acid-base catalysis of the hairpin ribozyme. *RNA* **2009**, *15*, 560–575.

(26) Heldenbrand, H.; Janowski, P. A.; Giambaşu, G.; Giese, T. J.; Wedekind, J. E.; York, D. M. Evidence for the Role of Active Site Residues in the Hairpin Ribozyme from Molecular Simulations along the Reaction Path. *J. Am. Chem. Soc.* **2014**, *136*, 7789–7792.

(27) Mlýnský, V.; Banáš, P.; Hollas, D.; Réblová, K.; Walter, N. G.; Šponer, J.; Otyepka, M. Extensive Molecular Dynamics Simulations Showing That Canonical G8 and Proto-nated A38H + Forms Are Most Consistent with Crystal Structures of Hairpin Ri-bozyme. *J. Phys. Chem. B* **2010**, *114*, 6642–6652.

(28) Mlýnský, V.; Kührová, P.; Zgarbová, M.; Jurečka, P.; Walter, N. G.; Otyepka, M.; Šponer, J.; Banáš, P. Reactive Conformation of the Active Site in the Hair-pin Ribozyme Achieved by Molecular Dynamics Simulations with ε/ζ Force Field Reparametrizations. *J. Phys. Chem. B* **2015**, *119*, 4220–4229.

(29) Kumar, N.; Marx, D. Mechanistic role of nucleobases in self-cleavage catalysis of hair-pin ribozyme at ambient *versus* high-pressure conditions. *Phys. Chem. Chem. Phys.* **2018**, *20*, 20886–20898.

(30) Veenis, A. J.; Li, P.; Soudackov, A. V.; Hammes-Schiffer, S.; Bevilacqua, P. C. Investigation of the pKa of the Nucleophilic O2' of the Hairpin Ribozyme. *J. Phys. Chem. B* **2021**, *125*, 11869–11883.
(31) Gaines, C. S.; York, D. M. Ribozyme Catalysis with a Twist: Active State of the Twister Ribozyme in Solution Predicted from Molecular Simulation. *J. Am. Chem. Soc.* **2016**, *138*, 3058–3065.

(32) Šponer, J.; Bussi, G.; Krepl, M.; Banáš, P.; Bottaro, S.; Cunha, R. A.; Gil-Ley, A.; Pinamonti, G.; Poblete, S.; Jurečka, P.; Walter, N. G.; Otyepka, M. RNA Structural Dynamics As Captured by Molecular Simulations: A Comprehensive Overview. *Chem. Rev.* **2018**, *118*, 4177–4338.

(33) Liebl, K.; Zacharias, M. The development of nucleic acids force fields: From an unchallenged past to a competitive future. *Biophys. J.* **2023**, *122*, 2841–2851.

(34) Gaines, C. S.; Piccirilli, J. A.; York, D. M. The L-platform/L-scaffold framework: a blueprint for RNA-cleaving nucleic acid enzyme design. *RNA* **2020**, *26*, 111–125.

(35) Wilson, T. J.; Lilley, D. M. Do the hairpin and VS ribozymes share a common catalytic mechanism based on general acid-base catalysis? A critical assessment of available experimental data. *RNA* **2011**, *17*, 213–221.

(36) Mlýnský, V.; Banáš, P.; Walter, N. G.; Šponer, J.; Otyepka, M. QM/MM Studies of Hairpin Ribozyme Self-Cleavage Suggest the Feasibility of Multiple Competing Reaction Mechanisms. *J. Phys. Chem. B* **2011**, *115*, 13911–13924.

(37) Nam, K.; Gao, J.; York, D. M. Quantum Mechanical/Molecular Mechanical Simulation Study of the Mechanism of Hairpin Ribozyme Catalysis. *J. Am. Chem. Soc.* **2008**, *130*, 4680–4691.

(38) Nam, K.; Gao, J.; York, D. M. Electrostatic interactions in the hairpin ribozyme account for the majority of the rate acceleration without chemical participation by nucleobases. *RNA* **2008**, *14*, 1501–1507.
(39) Mlýnský, V.; Banáš, P.; Šponer, J.; van der Kamp, M. W.; Mulholland, A. J.; Otyepka, M. Comparison of ab Initio, DFT, and Semiempirical QM/MM Approaches for Description of Catalytic Mechanism of Hairpin Ribozyme. *J. Chem. Theory Comput.* **2014**, *10*, 1608–1622.

(40) Kumar, N.; Marx, D. Deciphering the Self-Cleavage Reaction Mechanism of Hairpin Ribozyme. *J. Phys. Chem. B* **2020**, *124*, 4906–4918.

(41) Torelli, A. T.; Krucinska, J.; Wedekind, J. E. A comparison of vanadate to a 2'-5' linkage at the active site of a small ribozyme suggests a role for water in transition-state stabilization. *RNA* **2007**, *13*, 1052–1070.

(42) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. *SoftwareX* **2015**, *1-2*, 19–25.

(43) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. *Comput. Phys. Commun.* **2014**, *185*, 604–613.

(44) The PLUMED consortium, Promoting transparency and reproducibility in enhanced molecular simulations. *Nat. Methods* **2019**, *16*, 670–673.

(45) Humphrey, W.; Dalke, A.; Schulten, K. VMD – Visual Molecular Dynamics. *J. Mol. Graphics* **1996**, *14*, 33–38.

(46) Hénin, J.; Lopes, L. J. S.; Fiorin, G. Human Learning for Molecular Simulations: The Collective Variables Dashboard in VMD. *J. Chem. Theory Comput.* **2022**, *18*, 1945–1956.

(47) Michaud-Agrawal, N.; Denning, E. J.; Woolf, T. B.; Beckstein, O. MDAnalysis: A toolkit for the analysis of molecular dynamics simulations. *J. Comput. Chem.* **2011**, *32*, 2319–2327.
(48) Gowers, R. J.; Linke, M.; Barnoud, J.; Reddy, T. J. E.; Melo, M. N.; Seyler, S. L.; Domanski, J.; Dotson, D. L.; Buchoux, S.; Kenney, I. M.; Beckstein, O. MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations. **2019**, 98–105.

(49) Bottaro, S.; Bussi, G.; Pinamonti, G.; Reißer, S.; Boomsma, W.; Lindorff-Larsen, K. Barnaba: software for analysis of nucleic acid structures and trajectories. *RNA* **2019**, *25*, 219–231.

(50) Salter, J.; Krucinska, J.; Alam, S.; Grum-Tokars, V.; Wedekind, J. E. Water in the Active Site of an All-RNA Hairpin Ribozyme and Effects of Gua8 Base Variants on the Geometry of Phosphoryl Transfer. *Biochemistry* **2006**, *45*, 686–700.

(51) Cheatham, T. E.; Cieplak, P.; Kollman, P. A. A Modified Version of the Cornell *et al*. Force Field with Improved Sugar Pucker Phases and Helical Repeat. *J. Biomol. Struct. Dyn.* **1999**, *16*, 845–862.

(52) Pérez, A.; Marchán, I.; Svozil, D.; Sponer, J.; Cheatham, T. E.; Laughton, C. A.; Orozco, M. Refinement of the AMBER Force Field for Nucleic Acids: Improving the Description of α/γ Conformers. *Biophys. J.* **2007**, *92*, 3817–3829.

(53) Zgarbová, M.; Otyepka, M.; Šponer, J.; Mládek, A.; Banáš, P.; Cheatham, T. E.; Jurečka, P. Refinement of the Cornell et al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles. *J. Chem. Theory Comput.* **2011**, *7*, 2886–2902.

(54) Zgarbová, M.; Luque, F. J.; Šponer, J.; Cheatham, T. E.; Otyepka, M.; Jurečka, P. Toward Improved Description of DNA Backbone: Revisiting Epsilon and Zeta Torsion Force Field Parameters. *J. Chem. Theory Comput.* **2013**, *9*, 2339–2354.

(55) Fröhlking, T.; Mlýnský, V.; Janeček, M.; Kührová, P.; Krepl, M.; Banáš, P.; Šponer, J.;
Bussi, G. Automatic Learning of Hydrogen-Bond Fixes in the AMBER RNA Force Field. *J. Chem. Theory Comput.* **2022**, *18*, 4490–4502.

(56) Joung, I. S.; Cheatham, T. E. Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. *J. Phys. Chem. B* **2008**, *112*, 9020–9041.

(57) Tan, D.; Piana, S.; Dirks, R. M.; Shaw, D. E. RNA force field with accuracy comparable to state-of-the-art protein force fields. *Proc. Natl. Acad. Sci.* **2018**, *115*, 1346–1355.

(58) Piana, S.; Robustelli, P.; Tan, D.; Chen, S.; Shaw, D. E. Development of a Force Field for the Simulation of Single-Chain Proteins and Protein-Protein Complexes. *J. Chem. Theory Comput.* **2020**, *16*, 2494–2507.

(59) Tucker, M. R.; Piana, S.; Tan, D.; LeVine, M. V.; Shaw, D. E. Development of Force Field Parameters for the Simulation of Single- and Double-Stranded DNA Molecules and DNA-Protein Complexes. *J. Phys. Chem. B* **2022**, *126*, 4442–4457.

(60) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. *The Journal of Chemical Physics* **1983**, *79*, 926–935.

(61) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. *J. Chem. Phys.* **2007**, *126*, 014101.

(62) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. *J. Appl. Phys.* **1981**, *52*, 7182–7190.

(63) Nosé, S.; Klein, M. Constant pressure molecular dynamics for molecular systems. *Mol. Phys.* **1983**, *50*, 1055–1076.

(64) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An $N \log(N)$ method for Ewald sums in large systems. *J. Chem. Phys.* **1993**, *98*, 10089–10092.
(65) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A linear constraint solver for molecular simulations. *J. Comput. Chem.* **1997**, *18*, 1463–1472.

(66) Liu, P.; Kim, B.; Friesner, R. A.; Berne, B. J. Replica exchange with solute tempering: A method for sampling biological systems in explicit water. *Proc. Natl. Acad. Sci.* **2005**, *102*, 13749–13754.

(67) Wang, L.; Friesner, R. A.; Berne, B. J. Replica Exchange with Solute Scaling: A More Efficient Version of Replica Exchange with Solute Tempering (REST2). *J. Phys. Chem. B* **2011**, *115*, 9431–9438.

(68) Bottaro, S.; Di Palma, F.; Bussi, G. The role of nucleobase interactions in RNA structure and dynamics. *Nucleic Acids Research* **2014**, *42*, 13306–13314.

(69) Murchie, A. I.; Thomson, J. B.; Walter, F.; Lilley, D. M. Folding of the Hairpin Ribozyme in Its Natural Conformation Achieves Close Physical Proximity of the Loops. *Mol. Cell* **1998**, *1*, 873–881.

(70) Kuzmin, Y. I.; Da Costa, C. P.; Cottrell, J. W.; Fedor, M. J. Role of an Active Site Adenine in Hairpin Ribozyme Catalysis. *J. Mol. Biol.* **2005**, *349*, 989–1010.

(71) Pinard, R. Functional involvement of G8 in the hairpin ribozyme cleavage mechanism. *The EMBO Journal* **2001**, *20*, 6434–6442.

(72) Kuzmin, Y. I.; Da Costa, C. P.; Fedor, M. J. Role of an Active Site Guanine in Hairpin Ribozyme Catalysis Probed by Exogenous Nucleobase Rescue. *J. Mol. Biol.* **2004**, *340*, 233–251.

(73) Lebruska, L. L.; Kuzmine, I. I.; Fedor, M. J. Rescue of an Abasic Hairpin Ribozyme by Cationic Nucleobases. *Chem. Biol.* **2002**, *9*, 465–473.

(74) Kath-Schorr, S.; Wilson, T. J.; Li, N.-S.; Lu, J.; Piccirilli, J. A.; Lilley, D. M. J.
General Acid-Base Catalysis Mediated by Nucleobases in the Hairpin Ribozyme. *J. Am. Chem. Soc.* **2012**, *134*, 16717–16724.

(75) Thomas, J. M.; Perrin, D. M. Active Site Labeling of G8 in the Hairpin Ribozyme: Implications for Structure and Mechanism. *J. Am. Chem. Soc.* **2006**, *128*, 16540–16545.

(76) Bevilacqua, P. C. Mechanistic Considerations for General Acid-Base Catalysis by RNA: Revisiting the Mechanism of the Hairpin Ribozyme. *Biochemistry* **2003**, *42*, 2259–2265.

(77) Nahas, M. K.; Wilson, T. J.; Hohng, S.; Jarvie, K.; Lilley, D. M. J.; Ha, T. Observation of internal cleavage and ligation reactions of a ribozyme. *Nat. Struct. Mol. Biol.* **2004**, *11*, 1107–1113.

(78) Wilson, T. J.; Ouellet, J.; Zhao, Z.-Y.; Harusawa, S.; Araki, L.; Kurihara, T.; Lilley, D. M. Nucleobase catalysis in the hairpin ribozyme. *RNA* **2006**, *12*, 980–987.

(79) Liberman, J. A.; Guo, M.; Jenkins, J. L.; Krucinska, J.; Chen, Y.; Carey, P. R.; Wedekind, J. E. A Transition-State Interaction Shifts Nucleobase Ionization toward Neutrality To Facilitate Small Ribozyme Catalysis. *J. Am. Chem. Soc.* **2012**, *134*, 16933–16936.

(80) Liu, L.; Cottrell, J. W.; Scott, L. G.; Fedor, M. J. Direct measurement of the ionization state of an essential guanine in the hairpin ribozyme. *Nat. Chem. Biol.* **2009**, *5*, 351–357.

(81) Rhodes, M. M.; Reblova, K.; Sponer, J.; Walter, N. G. Trapped water molecules are essential to structural dynamics and function of a ribozyme. *Proc. Natl. Acad. Sci.* **2006**, *103*, 13380–13385.
(82) Abel, R.; Wang, L.; Harder, E. D.; Berne, B. J.; Friesner, R. A. Advancing Drug Discovery through Enhanced Free Energy Calculations. *Acc. Chem. Res.* **2017**, *50*, 1625–1632.

(83) Kamiya, M.; Sugita, Y. Flexible selection of the solute region in replica exchange with solute tempering: Application to protein-folding simulations. *J. Chem. Phys.* **2018**, *149*.

(84) Ramis, R.; Ortega-Castro, J.; Casasnovas, R.; Marino, L.; Vilanova, B.; Adrover, M.; Frau, J. A Coarse-Grained Molecular Dynamics Approach to the Study of the Intrinsically Disordered Protein α-Synuclein. *J. Chem. Inf. Model.* **2019**, *59*, 1458–1471.

(85) Reid, L. M.; Guzzetti, I.; Svensson, T.; Carlsson, A. C.; Su, W.; Leek, T.; von Sydow, L.; Czechtizky, W.; Miljak, M.; Verma, C.; De Maria, L.; Essex, J. W. How well does molecular simulation reproduce environment-specific conformations of the intrinsically disordered peptides PLP, TP2 and ONEG? *Chem. Sci.* **2022**, *13*, 1957–1971.

(86) Liu, X.; Li, X.; Qiao, Q.; Li, F.; Wei, G. ALS-Linked A315T and A315E Mutations Enhance β-Barrel Formation of the TDP-43307-319 Hexamer: A REST2 Simulation Study. *ACS Chem. Neurosci.* **2023**, *14*, 13101320.

(87) Kohn, E. M.; Konovalov, K.; Gomez, C. A.; Hoover, G. N.; Yik, A. K. H.; Huang, X.; Martell, J. D. Terminal Alkyne-Modified DNA Aptamers with Enhanced Protein Binding Affinities. *ACS Chem. Biol.* **2023**, *18*, 1976–1984.

(88) Languin-Cattoën, O.; Sterpone, F.; Stirnemann, G. Binding site plasticity regulation of the FimH catch-bond mechanism. **2023**, *122*, 2744–56.

(89) Maffucci, I.; Laage, D.; Sterpone, F.; Stirnemann, G. Thermal Adaptation of Enzymes: Impacts of Conformational Shifts on Catalytic Activation Energy and Optimum Temperature. *Chem. Eur. J.* **2020**, *26*, 10045–10056.
(90) Languin-Cattoën, O.; Melchionna, S.; Derreumaux, P.; Stirnemann, G.; Sterpone, F. Three Weaknesses for Three Perturbations: Comparing Protein Unfolding under Shear, Force, and Thermal Stresses. *J. Phys. Chem. B* **2018**, *122*, 11922–11930.

(91) Stirnemann, G.; Sterpone, F. Mechanics of protein adaptation to high temperatures. *J. Phys. Chem. Lett.* **2017**, *8*, 5884–5890.

(92) Stirnemann, G.; Sterpone, F. Recovering protein thermal stability using all-atom Hamiltonian replica-exchange simulations in explicit solvent. *J. Chem. Theo. Comput.* **2015**, *11*, 5573–5577.

(93) Röder, K.; Stirnemann, G.; Faccioli, P.; Pasquali, S. Computer-aided comprehensive explorations of RNA structural polymorphism through complementary simulation methods. *QRB Discovery* **2022**, *3*.

(94) Mlýnský, V.; Janeček, M.; Kührová, P.; Fröhlking, T.; Otyepka, M.; Bussi, G.; Banáš, P.; Šponer, J. Toward Convergence in Folding Simulations of RNA Tetraloops: Comparison of Enhanced Sampling Techniques and Effects of Force Field Modifications. *J. Chem. Theory Comput.* **2022**, *18*, 2642–2656.

(95) Röder, K.; Stirnemann, G.; Dock-Bregeon, A.-C.; Wales, D. J.; Pasquali, S. Structural transitions in the RNA 7SK 5 hairpin and their effect on HEXIM binding. *Nucleic Acids Res.* **2019**, *48*, 373–389.

(96) Kührová, P.; Best, R. B.; Bottaro, S.; Bussi, G.; Sponer, J.; Otyepka, M.; Banáš, P. Computer Folding of RNA Tetraloops: Identification of Key Force Field Deficiencies. *J. Chem. Theory Comput.* **2016**, *12*, 4534–4548.

(97) Dabin, A.; Stirnemann, G. Toward a Molecular Mechanism of Complementary RNA Duplexes Denaturation. *J. Phys. Chem. B* **2023**, *127*, 6015–6028.
(98) Dabin, A.; Stirnemann, G. Atomistic simulations of RNA duplex thermal denaturation: Sequence- and forcefield-dependence. *Biophys. Chem. 2024*, **307**, 107167.

(99) Guo, F.; Yue, Z.; Trajkovski, M.; Zhou, X.; Cao, D.; Li, Q.; Wang, B.; Wen, X.; Plavec, J.; Peng, Q.; Xi, Z.; Zhou, C. Effect of Ribose Conformation on RNA Cleavage via Internal Tranesterification. *J. Am. Chem. Soc. 2018*, **140**, 11893–11897.

(100) Klein, D. J.; Ferré-D’Amaré, A. R. Structural Basis of *glmS* Ribozyme Activation by Glucosamine-6-Phosphate. *Science 2006*, **313**, 1752–1756.

(101) Ren, A.; Vušurović, N.; Gebetsberger, J.; Gao, P.; Juen, M.; Kreutz, C.; Micura, R.; Patel, D. J. Pistol ribozyme adopts a pseudoknot fold facilitating site-specific in-line cleavage. *Nat. Chem. Biol. 2016*, **12**, 702–708.

(102) Pinamonti, G.; Paul, F.; Noé, F.; Rodriguez, A.; Bussi, G. The mechanism of RNA base fraying: Molecular dynamics simulations analyzed with core-set Markov state models. *J. Chem. Phys. 2019*, **150**, 154123.

(103) Pietrek, L. M.; Stelzl, L. S.; Hummer, G. Hierarchical Assembly of Single-Stranded RNA. *J. Chem. Theory Comput.* **2024**,

(104) Ashwood, B.; Jones, M. S.; Radakovic, A.; Khanna, S.; Lee, Y.; Sachleben, J. R.; Szostak, J. W.; Ferguson, A. L.; Tokmakoff, A. Direct monitoring of the thermodynamics and kinetics of DNA and RNA dinucleotide dehybridization from gaps and overhangs. *bioRxiv 2023*,

(105) Kührová, P.; Mlýnský, V.; Zgarbová, M.; Krepl, M.; Bussi, G.; Best, R. B.; Otyepka, M.; Šponer, J.; Banáš, P. Improving the Performance of the Amber RNA Force Field by Tuning the Hydrogen-Bonding Interactions. *J. Chem. Theory Comput.* **2019**, *15*, 3288–3305.
(106) Nesbitt, S.; Hegg, L. A.; Fedor, M. J. An unusual pH-independent and metal-ion-independent mechanism for hairpin ribozyme catalysis. *Chem. Bio.* **1997**, *4*, 619-630.

(107) Gaines, C. S.; Giese, T. J.; York, D. M. Cleaning Up Mechanistic Debris Generated by Twister Ribozymes Using Computational RNA Enzymology. *ACS Catal.* **2019**, *9*, 5803-5815.

(108) Ganguly, A.; Weissman, B. P.; Giese, T. J.; Li, N.-S.; Hoshika, S.; Rao, S.; Benner, S. A.; Piccirilli, J. A.; York, D. M. Confluence of theory and experiment reveals the catalytic mechanism of the Varkud satellite ribozyme. *Nat. Chem.* **2020**, *12*, 193-201.

(109) Benayad, Z.; Bova Saint-André, M.; Stirnemann, G. Molecular Mechanisms of Phosphoester Bond Formation in Water Using Tight-Binding Ab Initio Molecular Dynamics. *J. Phys. Chem. B* **2022**, *126*, 8251-8265.

(110) Benayad, Z.; David, R.; Stirnemann, G. Molecular mechanisms of phosphoester bond formation in abiotic conditions with reactive neural network potentials. *ChemRxiv* **2023**.

(111) Florián, J.; Warshel, A. Phosphate Ester Hydrolysis in Aqueous Solution: Associative versus Dissociative Mechanisms. *J. Phys. Chem. B* **1998**, *102*, 719-734.

(112) Duarte, F.; Åqvist, J.; Williams, N. H.; Kamerlin, S. C. L. Resolving apparent conflicts between theoretical and experimental models of phosphate monoester hydrolysis. *J. Am. Chem. Soc.* **2015**, *137*, 1081-1093.

(113) Warshel, A.; Sharma, P. K.; Kato, M.; Xiang, Y.; Liu, H.; Olsson, M. H. M. Electrostatic Basis for Enzyme Catalysis. *Chem. Rev.* **2006**, *106*, 3210-3235.

(114) Maffucci, I.; Laage, D.; Stirnemann, G.; Sterpone, F. Differences in thermal structural
changes and melting between mesophilic and thermophilic dihydrofolate reductase
enzymes. *Phys. Chem. Chem. Phys.* **2020**, *22*, 18361–18373.
<figure><img src="image_13.png" /><figcaption>Figure 10: TOC Graphics</figcaption></figure>