# A multiscale study to characterize SHAPE probe/ RNA molecules pre-reactive complex and its role in initiating the SHAPE reaction

Cécilia Hognon, Ameni Ben Abdeljaoued, Pierre Hardouin, Mélanie Etheve-Quelquejeu, Bruno Sargueil, Damien Laage, Elisa Frezza

► To cite this version:

Cécilia Hognon, Ameni Ben Abdeljaoued, Pierre Hardouin, Mélanie Etheve-Quelquejeu, Bruno Sargueil, et al.. A multiscale study to characterize SHAPE probe/ RNA molecules pre-reactive complex and its role in initiating the SHAPE reaction. Journal of Chemical Information and Modeling, 2025, 10.1021/acs.jcim.5c00812. [hal-05226372](https://hal.science/hal-05226372v1)

HAL Id: hal-05226372

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

Submitted on 27 Aug 2025

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>

Distributed under a Creative Commons Attribution 4.0 International License
# A multiscale study to characterize SHAPE
probe/ RNA molecules pre-reactive complex and
its role in initiating the SHAPE reaction

Cécilia Hognon,†,§ Ameni Ben Abdeljaoued,†,§ Pierre Hardouin,† Mélanie
Etheve-Quelquejeu,‡,|| Bruno Sargueil,† Damien Laage,¶ and Elisa Frezza*,†

†Université Paris Cité, CiTCoM, CNRS, F-75006 Paris, France

‡Université Paris Cité, LCBPT, CNRS, F-75006 Paris, France

¶Laboratory CPCV, Département de Chimie, École Normale Supérieure, PSL University,

Sorbonne University, CNRS, F-75005 Paris, France

§The authors contributed equally to the work

||Institut Universitaire de France

E-mail: elisa.frezza@u-paris.fr

Phone: +33 (0)1 53 73 15 69

## Abstract

Ribonucleic acid (RNA) molecules play a crucial role in nearly every cellular process, with their function closely tied to their three-dimensional (3D) structure. As a result, determining the precise 3D structure of RNAs is essential to understand their biological functions. However, obtaining high-resolution 3D structures remains a significant challenge using traditional biophysical techniques. To address this, chemical probing methods such as SHAPE (Selective 2'-Hydroxyl Acylation analyzed by Primer Extension) have gained widespread popularity. SHAPE reactivities have been introduced in
2D RNA predictors as soft constraints on nucleotide base pairing although the acylation reaction involves the ribose moiety. Little is known about the physical chemistry behind this reaction, leaving several reactivities unexplained. In this context, our aim is to unveil the complex relationship between the local structure of RNA, its dynamics, and the SHAPE chemical reactivity. In this study, using a multiscale approach based on biased molecular dynamics simulations and quantum mechanics/molecular mechanics calculations on the well characterized GAAA RNA tetraloop, we provide new molecular insights on the pre-reactive complex and its binding mode. Our findings emphasize the critical role of the local environment in facilitating the recruitment and proper accommodation of the SHAPE probe. All-atom umbrella sampling (US) molecular dynamics (MD) simulations underscored the significance of the binding angle in forming the pre-reactive complex, which was thoroughly analyzed using geometrical descriptors. In US QM/MM simulations, we used different initial structures to demonstrate the impact of favorable binding on the acylation reaction by characterizing the first step. Moreover, we investigated the deprotonation of the hydroxyl and its role in the acylation reaction, highlighting the role of the local environment and the probe proximity in this process. Our results suggest that the formation of oxyanion prior to the binding is not required to initiate the chemical reaction.
Introduction

Ribonucleic acid (RNA) molecules are involved in a large range of biological processes, in-
cluding gene expression, regulation and biochemical reactions such as protein synthesis or
pre-mRNA (pre-messenger RNA) splicing. ¹,² RNA function relies not only on its sequence
and secondary structure (2D) but also on its unique 3D fold and structural dynamics. ³⁻⁵
However, obtaining high-resolution RNA 3D structures remains challenging due for exam-
ple to the high flexibility and reactivity of RNA, as reflected by the limited number of
RNA 3D structures available in the Protein Data Bank. ⁶,⁷ In this context, chemical prob-
ing techniques are widely used to study RNA molecules. ⁸⁻¹¹ Among these, SHAPE (Selec-
tive 2'-Hydroxyl Acylation analyzed by Primer Extension) is particularly popular due to
its nucleotide-independent nature and compatibility with high-throughput protocols. ⁹,¹²,¹³
SHAPE provides nucleotide-level RNA structure insights using small-molecule electrophilic
probes such as 1-methyl-7-nitroisatoic anhydride⁹ or benzoyl cyanide (BzCN), ¹⁴ which acy-
late the ribose 2'-hydroxyl group to form a 2'-O-adduct. ¹⁵⁻¹⁷ SHAPE data have been inte-
grated as soft constraints to improve 2D structure predictions, a key step for 3D structure
modeling, under the assumption that high reactivity indicates the absence of Watson-Crick
base pairing at the nucleotide. ¹²,¹⁸⁻²⁰ However, SHAPE reactivities are related to nucleotide
flexibility, in particular at the level of ribose, and are dependent on the RNA motifs. ²¹,²²
Although these molecules preferentially react when the nucleotides are single-stranded (ss),
it is impossible to know whether that base is paired or unpaired, especially at lower reactivi-
ties, as confirmed by the fact that several reactivities are not explained when compared with
known experimental structures. ²³ Furthermore, the variation in reactivity for the same nu-
cleotide with different probes remains unclear. ²³ These observations are only a consequence
of a more complex mechanism whose origins have not yet been fully clarified. Moreover,
from the chemical point of view, different mechanisms can be proposed: 1)Moreover, from
the chemical point of view, different mechanisms can be proposed: 1) a direct attack of the
2'-OH onto the electrophilic carbonyl of the SHAPE reagent either after a pre-equilibrium
forming the reactive oxyanion (deprotonated by the buffer as proposed by Kool and co-workers)$^{15}$ or the formation of the 2'-oxyanion following a deprotonation step; 2) a general base catalysis (catalysis by the pyrimidine O2, purine N3, pro-S oxygen, and through-space groups$^{17}$). Despite this, it is important to note that the hydroxyl $pK_a$ is between ~ 12 and 14;$^{24}$ therefore, the local environment and the ribose conformation could play a crucial role in the formation of the oxyanion and change the value of its effective $pK_a$. $^{25}$ In addition, it is important to note that the presence of the oxyanion may significantly impact RNA stability by promoting backbone cleavage. This has been demonstrated in several studies on RNA hydrolysis and ribozymes,$^{26,27}$ which highlight the potential instability of the oxyanion. This shows once more the complexity of SHAPE reactivity and the necessity of gaining a deeper understanding of its chemical origins.

To tackle this problem, in recent years, some research groups, including ours, have focused on correlating SHAPE data with structural fluctuations and parameters obtained through MD simulations of small RNA molecules, utilizing both all-atom and simplified representations.$^{22,23,28-30}$ These studies underscore the critical role of flexibility, particularly ribose puckering, in SHAPE reactivity. Although flexibility is essential, it alone cannot fully explain the observed variations in SHAPE reactivities across different nucleotides and probes.

Another approach to explore the origin of biochemical reactivity involves explicitly taking into account the SHAPE reagent in molecular dynamics simulations, as proposed by Bussi and co-workers.$^{21,31}$ Using all-atom molecular dynamics simulations and enhanced sampling techniques, they investigated BzCN binding to an SRP-RNA motif, highlighting the critical role of local RNA conformations in enabling reagent access to the 2'-OH group that is related to larger SHAPE reactivities and the importance of sugar puckering as a key determinant of SHAPE reactivity. However, this study could only infer pre-reactive conformations based on the distance between the reactive carbon (C*) of BzCN and the O2' oxygen in free energy profiles, without considering other geometrical descriptors to characterize the binding process and distinguish binding modes. Furthermore, to be able to use and combine SHAPE
reactivities obtained with different reagents and unveil the physical chemical origin of its reactivity for RNA 3D structure prediction, it is essential to know the most plausible pre-reactive conformation, in particular its relative orientation with respect to the ribose and the nucleobase of interest since it can influence the free-energy barrier of the reaction. In fact, this conformation is crucial to access the reaction mechanism and sample reaction pathways with its respective free-energy barriers and reaction free energy in order to identify the more plausible path, using for example quantum mechanics/molecular mechanics (QM/MM) approaches.

To fill this gap we considered the GAAA tetraloop, a member of the GNRA tetraloop family which is overrepresented in natural RNA. The structure at the atomic level of such has been characterized many times by NMR, X-Ray crystallography and CryoEM. As its structure relies on non-canonical dynamic interactions, it is difficult to anticipate its SHAPE reactivity profile. Here we investigated the binding free energy profiles for BzCN by taking into account not only the distance as proposed in ref.,²¹ but also a specific angle that describes how the reagent approaches the ribose. By comparing these results with the new experimental data obtained in this work and performing a detailed structural analysis on our simulations, we were able to provide a new molecular-level rationalization of the SHAPE reactivities that goes beyond the capability of the reagent to be close to the nucleobase, allowing us to identify the most favorable conformations of the pre-reactive complex. To validate the importance of these pre-reactive conformations, we performed for the first time on this type of systems QM/MM umbrella sampling simulations using the distance to sample the formation of the adduct and without taking into account the rest of the reaction steps. Our results pave the way to the study of SHAPE reagents reaction mechanisms and pathways with respect to local geometries, a step towards using SHAPE data to predict RNA 3D structures.
# Material and Methods

## System under investigation

To conduct this study, we focused on *Bacillus subtilis* yitJ S box/SAM-I riboswitch RNA, whose structure is available under the PDB ID: 4KQY.³² To compare our results with the experimental ones, we removed the S-adenosyl methionine (SAM) from the crystallography structure. As presented in Figure 1A and B, the aptamer domain of the riboswitch is constituted by four stem-loops called P1 to P4, spatially arranged by 3D interactions. In this work, we mostly focus on the P3a stem-loop as explained below.

## Force field parameterization

The force field for the SHAPE BzCN probe was parameterized following the Restrained Electrostatic Potential (REsP) procedure.³³ The geometry of the SHAPE probe was optimized at the density functional theory (DFT) level, followed by the calculation of the electrostatic potential at Hartree-Fock level with the standard 6-31G* basis set³⁴,³⁵ and the B3LYP functional³⁶ using Gaussian16 software.³⁷ Other parameters were determined using the Generalized Amber Force Field.³⁸ Parameters are available in Table S1 and Figure S1.

## All-atom molecular dynamics simulations

We performed all-atom MD simulations using the GROMACS 2018 package³⁹⁻⁴² with the Amber ff99bsc0 force field with OL3 modification (ff99bsc0χ<sub>OL3</sub>) for RNA molecules.⁴³ The different systems were placed in a cubic box and solvated with TIP4P-EW to a depth of at least 12 Å.⁴⁴ To ensure electroneutrality, K⁺ ions were added to balance the global negative charge and to reach a physiological salt concentration of 0.15 M, K⁺Cl⁻ ion pairs were also added (around 310 ions for the whole structure and around 70 ions for the isolated stem-loop). A cutoff of 10 Å was used and the particle mesh Ewald method⁴⁵,⁴⁶ was employed to treat the long-range electrostatic interactions. The covalent bond lengths involving hydrogen atoms
<figure><img src="image_2.png" /><figcaption>Figure 1: a) 3D structure and b) 2D structure of the *Bacillus subtilis* yitJ S box/SAM-I riboswitch (PDB ID: 4KQY). The grey dotted lines represent the pseudoknots. Loop P1 is in blue, P2 is in sand, P3 is in green, P3a is in dark green and P4 is red.</figcaption></figure>

were restrained using P-LINCS,⁴⁷ allowing a time step of 2 fs. Every 1000 steps, translational movement of the solute was excluded to prevent the accumulation of kinetic energy.⁴⁸

Before the MD production, we carried out the energy minimization and equilibration. During equilibration, at least 10 ns, a Berendsen thermostat ($\tau_T$ = 0.1 ps) and Berendsen pressure coupling ($\tau_P$ = 0.5 ps)⁴⁹ were used. The production part was carried out in an NTP ensemble at a temperature held at 300 K and a pressure held constant at 1 bar using the Bussi velocity-rescaling thermostat ($\tau_T$ = 0.1 ps)⁵⁰ and the Parrinello-Rahman barostat ($\tau_P$ = 0.5 ps).⁵¹ By using positional restraints during minimization and heating, RNA heavy atoms remained fixed and were gradually relaxed along the equilibration. During the last equilibration, where the positional restraints are removed, the three replica were computed using different random initial velocities with *gen-vel* Gromacs option set to yes and *gen-seed* to -1. Moreover, to study the part of the loop P3, P3a, of the *Bacillus subtilis* yitJ S box/SAM-I riboswitch RNA, we applied constraints on the terminal bases, A56 and U69, during all the steps, including MD productions. This ensured the preservation of Watson-Crick hydrogen bonds at the tetraloop ends and prevented any base-pair opening that could
destabilize the structure. For the systems with the whole RNA and P3a stem-loop, MD production were carried out for 1 µs in three replicas and 100 ns in three replica on the system comprising P3a loop and one probe.

## Umbrella sampling in all-atom MD simulations

To conduct umbrella sampling (US) simulations using GROMACS 2018 package³⁹⁻⁴² and PLUMED,⁵² we prepared four different starting structures, depending the nucleotide under investigation. We particularly focused our investigation on the nucleotides present in the tetraloop of the P3a, G61, A62, A63 and A64. In all starting structures, we placed the SHAPE probe at a distance of 14 Å between the reactive carbon of SHAPE probe and the O2' of the ribose. To do this, we used the equilibrated structures extracted from the classical MD simulation on the P3a loop and followed the protocol previously described. We only change, to reduce the computational cost, the type of water box and we solvated the system in a rhombic dodecahedron water box. After the equilibration procedure, the SHAPE probe has been steered to create each window.

The 1D free-energy profiles were obtained sampling the collective variable (CV) defined by the variation of the distance between the reactive oxygen of the ribose of the nucleotide under investigation and the probe's reactive carbon C*. It was calculated in the distance *d* interval comprised between 14 Å to 2 Å (Figure 2a). The instantaneous force values were accumulated in bins of 0.1 Å. For each system, the US simulation was therefore constituted by 121 windows of 5 ns long, representing a total of 605 ns. For comparison, for the nucleotide 62, we also computed 1D replica exchange US simulations using the same protocol with an exchange between the replica each 200 steps. Then, the 2D free energy profiles were also computed sampling the CV previously defined for the 1D profile and the angle between the probe's reactive carbon and C4' and C2' atoms of the nucleotide's ribose, represented in Figure 2b. Depending of the results highlighted by the 1D free energy profiles, the distance interval was comprised between 10 Å to 2 Å or between 8 Å to 2 Å with bins of 0.1 Å.
Regarding the angle collective variable $\alpha$, it was calculated from 0° to 130°, incrementing by 10°. This represented between 854 and 1134 windows of 5 ns long, for a total comprised between 4.27 µs and 5.67 µs. For 1D and 2D US simulations, the potential restraint constant $k_d$ applied along the steered MD simulation and the US production was of 1.08 kcal mol$^{-1}$ Å$^{-2}$ for distance CV and 0.0327 kcal mol$^{-1}$ °$^{-2}$ for angle CV. Finally, the potential of mean force (PMF) profiles were built using the weighted histogram analysis method (WHAM) developed by the Grossfield lab.⁵³ The stability of the PMF profile regardless of the number of data points served as a convergence criteria.

<figure><img src="image_3.png" /><figcaption>Figure 2: a: Definition of the distance $d$ between the reactive carbon (C*) and the O2' of the ribose as collective variable $CV_d$. b: Definition of the valence angle $\alpha$ formed by the reactive carbon C*, the C4' and C2' atoms of the ribose, as collective variable $CV_\alpha$. c: Schematic SHAPE reaction with protonated hydroxyl. d: Schematic SHAPE reaction where the oxyanion explicitly shown.</figcaption></figure>
## NEB-TS calculations

To get insights on the reaction mechanism on single nucleotides, we used the NEB-TS method (for Nudged Elastic Band with transition state (TS) optimization) to perform QM calculations at the level of theory DFT/B3LYP/DEF2-SVP/D4 as implemented in ORCA.⁵⁴ NEB-TS allows to find the minimum energy path of a reaction between reactant and final product, while determining the transition state. We conducted the calculations in gas or using an implicit solvent model. To include solvent effects in the calculations, we used the Universal Solvation Model (SMD) as implicit solvent model.⁵⁵ Here, we investigated the two proposed mechanisms as presented in Figure 2 c and d. To do so, we considered a single nucleotide as in the previous work by Jash and Kool¹⁵ (see Figure 2 d). At first, we performed the geometry optimization of the reactant and the product with the same parameters chosen for the NEB-TS calculations. Second, by providing the optimized structures of the reactant and the reaction product and choosing the number of images *N*, the software automatically generated *N* structures along the reaction path, which it optimized while searching for the structure corresponding to the transition state. Here, we chose 15 images for our calculations (*N*=15).

## Quantum Mechanics/Molecular Mechanics (QM/MM) umbrella sampling calculations

QM/MM calculations were performed using the Sander software of Amber 20 package on the same stem-loop, P3a, and focusing on the nucleotide A62. First, the system was described at all-atom level using the same conditions described in the section above on all-atom MD simulations but utilizing the TIP3P water model. The probe was initially positioned at 4 Å from the nucleotide before switching to the QM/MM description (more details of all-atom system preparation are reported in SI). The QM region was described at the third version of Density Functional Tight Binding (DFTB3) level of theory using the up-to-date 3ob-
3-1 slater-koster parameter files 3ob-freq, 3ob-nhmod and 3ob-ophyd.⁵⁶⁻⁵⁹ The QM region consisted of the investigated nucleotide (A62) without the phosphodiester group at the 5' end, which was cut between C5' and C4' to add hydrogen link atoms. The adjacent nucleotide A63 was included in its entirety, along with the adjacent phosphodiester group of A64 (Figure S2), where a similar cut was made at C5'C4' bond to also incorporate hydrogen link atoms. This resulted in a total of 83 atoms with -2 charge. The rest of the system was described at all-atom level using the force field as detailed in the previous section. Umbrella sampling was conducted from 4 Å to 1 Å using the distance *d* between the reactive carbon C* and the O2' of the ribose under investigation as a collective variable. The constraint was a harmonic potential of 400 kcal mol⁻¹ Å⁻². We performed minimization, heating in NVT then a short equilibration in NPT every 0.5 Å window, and for each one we took the minimization output as a starting point for the next job. The minimization included 5000 cycles, 3000 in steepest descent and 2000 in conjugate gradient. The heating was performed from 0 K to 300 K in 5 ps followed by 5 ps at 300 K using Langevin dynamics with a collision frequency of 1 ps⁻¹. The equilibration in NPT was conducted for 6 ps by adding Monte Carlo barostat (pressure at 1 bar, coupling at 1 ps⁻¹). Final structures were used as a starting point for the neighboring windows in the production. For the umbrella sampling production, we opted for 250 ps long windows with a spacing of 0.05 Å, resulting in overall 61 windows, as a compromise between computational cost and good overlap in CV sampling. All QM/MM calculations were conducted with integration step of 1 fs. An additional constraint was added during minimization and heating to sample starting points at α around 80° (second free-energy basin) and α around 60° (barrier between the two basins) for the windows with *d* at 2.50, 3.00, 3.50, and 4.00 Å but was removed during NPT equilibration and the production. This additional bias was important to force the sampling of different binding poses since in our previous attempts, α values rapidly converged to the first basin as soon as we removed the bias at 4 Å.
To investigate binding at short distances along the collective variable (CV) *d* and the formation of C*-O2' adduct, we employed the protonated form of the 2'-hydroxyl group, rather than enforcing the pre-reactive deprotonated state proposed by Jash and Kool¹⁵ as detailed below from geometrical optimization. As a result, only QM/MM simulations were viable for probing reactivity, starting from a well-equilibrated protonated structure. The choice of the distance *d* as a CV allowed for a more inclusive exploration, enabling naturally sampling proton positions without artificially constraining the system. However, given the known *pK*<sub>a</sub> and the slow protonation equilibrium between the two states, and in an effort to rationalize differences observed in earlier QM/MM umbrella sampling simulations, additional simulations were conducted where we introduced four water molecules near the reactive site using a soft positional restraint (semi-harmonic potential) applied between their oxygen atoms and the ribose C2' atom to maintain their proximity. The extended CV includes the O2'-H2' distance to encourage deprotonation:

$$CV_2 = d_{C*O2'} - d_{O2'H2'} \quad (1)$$

To reduce artefacts in proton dynamics, integration time was set for this specific part to 0.5 fs.

## Structural analysis

To analyze MD simulations, different descriptors were applied and different Python and Visual Molecular Dynamics (VMD)⁶⁰ scripts and libraries were used. The plots were obtained using the software Matlab version R2022a.⁶¹

## RMSD/RMSF and Clustering analysis

The stability of the system was assessed by calculating Root-Mean-Square-Deviation (RMSD) time series for either the riboswitch or the P3a loop, using their initial structures as refer-
ences. RMSF (Root-Mean-Square-Fluctuations) were also determined for all atoms and representing the average fluctuations for each residue, to measure the average deviation of a molecule over time from a reference position. We computed these parameters using the GROMACS 2018 package.³⁹⁻⁴² Then, using the TTClust software,⁶² we conducted a hierarchical cluster analysis to identify clusters and extracted the most representative structure from the MD simulation or windows of the US simulation.

## Puckering analysis

To describe the ribose flexibility and its conformation, we computed the sugar puckering. It is important to note that the sugar ring is well described using pseudorotation parameters. While there are four possible pseudorotation parameters for a five-membered ring (Frezza et al.²⁸), two in particular are useful to characterize the sugar conformation: the so-called phase (Pha) and amplitude (Amp). While the amplitude describes the degree of ring puckering, the phase describes which atoms are most out of the mean ring plane. We calculated these parameters using the equations given below:

$$ \text{Amp} = \sqrt{(a^2 + b^2)}; \text{Pha} = \cos^{-1}\left(\frac{a}{\text{Amp}}\right) \quad (2) $$

where $a = 0.4 \sum_{i=1}^{5} \nu_i \cos[0.8\pi(i - 1)]$ and $b = -0.4 \sum_{i=1}^{5} \nu_i \sin[0.8\pi(i - 1)]$ with $\nu_i$ the ring dihedral angle $i$. This approach has the advantage of processing the ring dihedral angle $\nu_1$ (C1'-C2'-C3'-C4') to $\nu_5$ (O4'-C1'-C2'-C3') in an equivalent manner. Conventionally, sugar ring puckers are divided into 10 families described by the atom which is most displaced from the mean ring plane (C1', C2', C3', C4' or O4') and the direction of such displacement (endo for displacements on the side of the C5 atom and exo for displacements on the other side). In different sugar conformations, the distance between C2 atoms belonging to neighboring nucleobases ($d_{C2-C2}$) and the orientation of P relative to the sugar/bases are different from each other. Based on previous works, the fluctuations of the distance $d_{C2-C2}$ are supposed to
be related to the chemical reactivity of RNA with SHAPE reagents. In this context, it could be useful to group the sugar puckers into two large families. The sugar puckers C1'-exo, C2'-endo, C3' -exo, C4'-endo belong to the B-like family, while C1'-endo, C2'-exo, C3'-endo, C4'-exo belong to the A-like family. For the single angular quantities, because of their periodicity, standard statistics cannot be used to compute the average and the fluctuations of angular values. For example, the arithmetic mean of 0°and 360°is 180°, which is misleading because for most purposes 360°is the same thing as 0°. Hence, we used the directional statistics instead, as explained in.²⁸ To do so, the MD trajectories were analyzed using a Python script based on the packages mdtraj⁶³ and NumPy.⁶⁴

π-stacking analysis

To investigate the stacking between BzCN and the RNA, we considered the ribose and the nucleobase of the nucleotide under investigation in the simulation and its nearest neighbours. In the literature, several methods have been proposed to characterize stacking for example by employing specific descriptors⁶⁵ or distance-based feature vector.⁶⁶ In this work, we defined a simplified definition based on our previous work on protein-RNA complexes.⁶⁷ First, we computed the center of mass (COM) for the probe ($\mathbf{r}_{COM,BzCN}$) by taking into account the phenyl ring and the center of mass for the ribose belonging to the nucleotide $i(\mathbf{r}_{COM,r,i})$ and the nucleobase ($\mathbf{r}_{COM,b,i}$) by taking into account only the heavy atoms and the distance between the probe and either the ribose ($d_{stack,p-r} = \mathbf{r}_{COM,BzCN} - \mathbf{r}_{COM,r}$) or the nucleobase ($d_{stack,p-b} = \mathbf{r}_{COM,BzCN} - \mathbf{r}_{COM,b}$). For the nucleobase, we computed the angle between the normal vectors to the probe and the nucleobase $i$ ($\alpha_{stack}$) defined as follows:

$$ \alpha_{stack,p-b,i} = \arccos \frac{\mathbf{r}_{BzCN} \cdot \mathbf{r}_{b,i}}{|\mathbf{r}_{BzCN}||\mathbf{r}_{b,i}|} \quad (3) $$

where $\mathbf{r}_{BzCN} = (\mathbf{r}_{1,BzCN} - \mathbf{r}_{COM,BzCN}) \times (\mathbf{r}_{2,BzCN} - \mathbf{r}_{COM,BzCN})$ and $\mathbf{r}_{b,i} = (\mathbf{r}_{1,b,i} - \mathbf{r}_{COM,b,i}) \times (\mathbf{r}_{2,b,i} - \mathbf{r}_{COM,b,i})$ are the normal vectors to a ring of the BzCN and the nucleobase $i$, respectively.
tively. The vector $\mathbf{r}_{1,BzCN}$ represents the coordinates $(x, y, z)$ of the atom C5 and $\mathbf{r}_{2,BzCN}$ represents the coordinates $(x, y, z)$ of the atom C6 (see Figure S1). The vector $\mathbf{r}_{1,b,i}$ represents the coordinates $(x, y, z)$ of the atom N3 for G and A and the atom C2 for U and C and $\mathbf{r}_{2,b,i}$ represents the coordinates $(x, y, z)$ of the atom C6 for G and A and C4 for U and C.

For the ribose, since it is not planar, we computed the inertial tensor $\mathbf{I}$ by considering the atoms of the furane ring and we diagonalized the tensor to obtain the eigenvector $\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3$. Therefore, we calculated the angle between the normal vector of the probe ($\mathbf{r}_{BzCN}$) and the eigenvector $\mathbf{e}_r$ perpendicular to the plane of the ribose, as follows:

$$ \alpha_{stack,p-r,i} = \arccos \frac{(\mathbf{r}_{BzCN} \cdot \mathbf{e}_{r,i})}{|\mathbf{r}_{BzCN}|} \quad (4) $$

In addition to the stacking angles, we also projected the distances $d_{stack,p-b}$ (for the nucleobase) and $d_{stack,p-r}$ (for the ribose) onto their respective local coordinate systems to obtain the displacements in the $x, y$ and $z$ directions, denoted as $d_x, d_y$ and $d_z$, and we computed the displacement on the plane $xy$ defined as $d_{x,y} = \sqrt{(dx^2 + dy^2)}$. In conclusion, we considered that a stacking interaction is present if the distance $d_{stack,p-b}$ or $d_{stack,p-r}$ is below 5 Å and the absolute value of $\alpha_{stack,p-b,i}$ or $\alpha_{stack,p-r,i}$ is lower than 45°. Therefore, for each snapshot we can associate a value 0 or 1 based on the absence or the presence of the stacking respectively. All these analyses were conducted using a python script based on the packages `mdtraj`⁶³ and `NumPy`.⁶⁴

## Circular Variance

Inspired by previous works on proteins, we used the circular variance ($CV$), a simple measure based on the vectorial distribution of the neighbors around a point, to describe the exposure of particles to solvent.⁶⁸,⁶⁹ The $CV$ of an atom $i$ (denoted as $CV_i$) is defined as follows:

$$ CV = 1 - \frac{1}{n_i} \left| \sum_{j \neq i, |\mathbf{r}_{ij}| \le r_c} \frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|} \right| \quad (5) $$
where $r_c$ is a cutoff distance (here we chose values between 8 Å up to 12 Å), $n_i$ is the number of neighboring atom within the cutoff distance and $|\mathbf{r}_{ij}|$ is the distance vector between the atom $i$ and the neighboring atom $j$.

## eRMSD and dynamical secondary structure analysis

To quantify the difference in terms of the internal arrangements of bases of two given structures, we computed the eRMSD as defined in the work of Buttano et al⁷⁰ using PLUMED.⁵² The eRMSD can be considered as a vectorial version of contact maps. To calculate the eRMSD, a local reference system in the center of the six-membered ring of each nucleobase in a molecule has to be defined. To do so, the C2,C4,C6 atoms for pyrimidines and the C2, C6, C4 atoms for purines were selected as described previously in.⁷⁰ To cluster the simulations based on eRMSD values and to generate "dynamic secondary structure" figures, we also used the software Barnaba.⁷¹ These figures depict the extended secondary structure along with the population of each interaction within a collection of three-dimensional structures, illustrating the variability in annotations across the set of structures.

## Attack angles

In order to characterize the approach, we also computed the attack angles, $\beta_1$ and $\beta_2$. The former angle, $\beta_1$, also known as Bürgi–Dunitz angle, is defined as the Nu-C\*-O bonding angle (Figure 3a), with Nu the nucleophile group (the hydroxyl of the ribose). Angle $\beta_2$, similar to the Flippin–Lodge angle, is defined by the normal vector of the plane $B$ formed by the atoms C8, C\* and C5, centered on the reactive carbon C\* so that the normal vector points in the same direction of the Re face following the Cahn-Ingold-Prelog priority rules (Figure 3b and c), and the difference vector $d_{C^*-O_2'}$ defined as $d_{C^*-O_2'}$ with $\mathbf{r}_{C^*}$ representing the coordinates $(x, y, z)$ of the reactive carbon C\* belonging to BzCN and $\mathbf{r}_{O_2'}$ representing the coordinates $(x, y, z)$ of the oxygen O2' of the hydroxyl group of the ribose. Therefore, the
attack angle $\beta_2$ is defined:

$$ \beta_2 = \arccos \left( \frac{\mathbf{n}_B \cdot \mathbf{d}_{C^*-O2'}}{|\mathbf{n}_B||\mathbf{d}_{C^*-O2'}|} \right) \qquad (6) $$

where $\mathbf{n}_B$ is the normal vector of the plane $B$ defined as:

$$ \mathbf{n}_B = \mathbf{d}_{C^*-C5} \times \mathbf{d}_{C^*-C8} \qquad (7) $$

where $\mathbf{d}_{C^*-C5} = \mathbf{r}_{C*}$ with $\mathbf{r}_{C5}$ representing the coordinates $(x, y, z)$ of the reactive carbon C5 belonging to BzCN and $\mathbf{r}_{C8}$ representing the coordinates $(x, y, z)$ of the atom C8 belonging to BzCN.

<figure><img src="image_4.png" /><figcaption>Figure 3: Used nucleophilic attack angles $\beta_1$ and $\beta_2$: a) $\beta_1$ defined as the Nu-C*-O bonding angle; b) Definition of Re (face up when clockwise, face down when counterclockwise) Si (face down when clockwise, face up when counterclockwise) faces following Cahn-Ingold-Prelog priority rules; c) $\beta_2$ defined as the angle between the vector Nu-C* and the normal vector n of the plane C5C*C8 pointing to the Re face.</figcaption></figure>

## Characterization of QM/MM simulations

To further investigate differences arising from the distinct initial poses ($\alpha_i$) in QM/MM US simulations, particularly in relation to the nucleobases, we computed the distribution of $\kappa_i$. This parameter reflects the orientations of the nucleobases of the target nucleotide A62 (C1’P63), the adjacent nucleotide A63, and the probe BzCN (C*O2’) relative to the normal vector of a defined plane $R$. The C1’ ribose carbone attached directly to the nucleobase was used here to minimize noise from nucleobase flexibility and rotation. This analysis allowed
us to determine whether the investigated elements reside on the same side of the plane $\mathcal{R}$. Specifically, $\kappa_i$ is defined as the angle with respect to the normal vector of a plane formed by the phosphate atoms of the nucleotide G61, A63, and A64, which points outwards the RNA. If the nucleobases and the probe are on the same side, this will be reflected in the density distribution, where the data points will appear on the same side of the plot. The plane is defined using phosphate atoms of the backbone, with at least one phosphate not involved in the loop (the atom P of the nucleotide G61) to minimize mobility and fluctuations, ensuring a relatively stable reference for analysis. In addition, we calculated the angle $\beta_2$ (Figures 3a and b) in the different simulations to characterize the nucleophilic attack preferences in these different binding scenarios. The normal vector $\mathbf{n}_R$ of the plane $\mathcal{R}$ is defined:

$$ \mathbf{n}_R = \mathbf{d}_{P63-P64} \times \mathbf{d}_{P63-P64} \qquad (8) $$

where $\mathbf{d}_{P63-P64} = \mathbf{r}_{P63} - \mathbf{r}_{P64}$ with $\mathbf{r}_{P63}$ representing the coordinates $(x, y, z)$ of the atom P belonging to the nucleobase A63 and $\mathbf{r}_{P64}$ representing the coordinates $(x, y, z)$ of the atom P belonging to the nucleobase A64, and $\mathbf{d}_{P63-P61} = \mathbf{r}_{P63} - \mathbf{r}_{P61}$ with $\mathbf{r}_{P61}$ representing the coordinates $(x, y, z)$ of the atom P belonging to the nucleobase G61. The angles $\kappa_i$ to for the nucleobase i (A62 or A63), along with BzCN are defined (see Figure 4):

$$ \kappa_{\text{nucleobases},i} = \arccos \left( \frac{\mathbf{n}_R \cdot \mathbf{d}_{C6-C1',i}}{|\mathbf{n}_R||\mathbf{d}_{C6-C1'}|} \right) ; \quad \kappa_{\text{BzCN},i} = \arccos \left( \frac{\mathbf{n}_R \cdot \mathbf{d}_{C^*-O2'}}{|\mathbf{n}_R||\mathbf{d}_{C^*-O2'}|} \right) \qquad (9) $$

where $\mathbf{d}_{C6-C1',i} = \mathbf{r}_{C6,i} - \mathbf{r}_{C1',i}$ with $\mathbf{r}_{C6,i}$ representing the coordinates $(x, y, z)$ of the atom C6 belonging to the nucleotide i and $\mathbf{r}_{C1',i}$ representing the coordinates $(x, y, z)$ of the atom C1' belonging to the nucleotide i and $\mathbf{d}_{C2*-O2',i} = \mathbf{r}_{C^*,i} - \mathbf{r}_{O2',i}$ with $\mathbf{r}_{O2',i}$ representing the coordinates $(x, y, z)$ of the atom O2' belonging to the nucleotide A62.

To characterize the local geometry of RNA in response to the presence of the probe at short range of the collective variable $d$ in our QM/MM US, we computed the angles $\tau$ and $\xi$. The valence angle $\xi$ was measured with respect to the nucleobases, specifically
a)

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

b)

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

Figure 4: Definition of κ angles. a) κ₁ and κ₂ as the angles between the normal vector **n** and the respective nucleobases of A62 (V1) and A63 (V2). b) κ₃ angles between the normal vector **n** and the distance vectors **d**<sub>C*-*O*₂</sub> showcasing differences resulting from which side the probe approaches.

between the atom C1' of the ribose directly attached to the nucleobase of A62 and the
atom C1' attached to the nucleobase of A63, both centered on the phosphate atom P of
the nucleotide A63. This angle provides insight into how the ribose is exposed from the
nucleobase side. The second valence angle, τ, is defined by the phosphate atoms P
belonging to the nucleotide G61, A63 and A64 to monitor the compensatory behavior in the
backbone. In addition to the angular measurements, we also computed the C6C6 distance
between the nucleobases of A62 and A63 to observe the widening or narrowing between the
bases distance, offering a quantitative indicator of the openness of the RNA structure. Since
C1' is attached to the backbone, while C6 is positioned farther away, the C6C6 coupled
with ξ provide valuable information about how the nucleobases move apart in response
to the probe's presence. Furthermore, phosphate-phosphate distances were computed to
investigate the RNA backbone, while the angle τ reflects the relative orientation of the RNA
backbone in relation to nucleobase movement. This allows us to observe the compensation
of the openness on the nucleobase side by corresponding adjustments on the backbone side.
Together, these geometrical descriptors enable a comprehensive understanding of how RNA's
local geometry responds to the probe present at short range, shedding light on the potential
conformational changes upon probe binding.

SHAPE experiments

SHAPE chemical probing was performed with BzCN at three different concentrations (20, 50 and 100mM) to consider the cooperative and concentration effects emphasized in the literature (see the SI for all the experimental details and Figure S3 for the BzCN reactivity on the 2D structure).<sup>23,31,72,73</sup> Experiments are conducted in triplicates to ensure the robustness of our findings. Raw data were treated using QShape<sup>74</sup> embedded in the IPANEMAP Suite workflow which outputs a mean reactivity value and a standard error for each nucleotide at a given porbe concentration<sup>75</sup> Finally IPANEMAP suite provides the reactivity displayed on the secondary structure model for each BzCN concentration (Figure S3).

**Results and Discussion**

In this study, we investigated the interactions between a specific SHAPE probe, BzCN, and
a well-characterized RNA molecule, the *Bacillus subtilis* yitJ S box/SAM-I riboswitch using
all-atom simulations and QM/MM calculations in relation to SHAPE experimental data as
a crucial step for its use in predicting the 3D structure of RNA.

For experimental reference in this study, we performed SHAPE structure probing at three
different concentrations of BzCN (see SI for details). In all-atom MD simulations, we did not
consider the oxyanion state, as the optimized structure of this state with the probe obtained
using ORCA showed a dislocation of the phosphate group — an expected outcome given
that the oxyanion form is known to favor RNA cleavage (see Figure S11). This limitation
was addressed using QM/MM simulations and by the choice of the first CV *d* that does not
restrain the system protonation state, thereby allowing an unbiased exploration of proton
dynamics. To further drive the system toward deprotonation and gain more understanding
of the differences in the QM/MM US outcomes, a second CV was introduced (eq. 1)

*Bacillus subtilis yitJ S box/SAM-I riboswitch characterization using all-atom MD simulations and SHAPE chemical probing*

Before investigating the interactions between the SHAPE probe and the RNA molecule,
three independent MD replicas were carried out in the absence of the probe with the aim
to characterize the RNA molecule, in terms of its stability, flexibility and the solvent acces-
sibility of the ribose, the key element in the acylation reaction. As a preliminary step, we
analyzed the unbiased MD trajectory by computing solvent accessibility of 2'-OH (ASA), a
previously proposed approach,²¹ as well as the circular variance (CV), an alternative metric
proposed to study the solvent exposition in proteins, using various distance cut-offs. Our re-
sults reported in the SI regarding ASA analysis are in agreement with the previously reported
analysis of different RNA molecules²¹,²²,²⁸ (Figure S6A and S6C). No significant correlation
was also observed for the circular variance, highlighting the difficulty to properly define the
local accessibility and the geometry of the ribose. Interestingly, a moderate correlation was
found for both fluctuations of ASA ($\rho_{\sigma_{asa}}$ = 0.4254) and CV ($\rho_{\sigma_{cv}}$ = 0.3380) (Figure S6B and
S6D). Contrarily to what is usually accepted, this suggests that the variability in accessibility
could be an important feature to consider when describing SHAPE reactivity, highlighting
again the importance of RNA dynamics and the need for molecular dynamics simulations.

To select an appropriate model system for further investigations of SHAPE binding and
reactivity, we performed a structural and stability analysis. The chosen model system should
fold independently and be free from 3D contacts with other parts of the molecule, allowing
us to focus on a simplified case while still extracting data from a larger biologically relevant
RNA molecule. To do so, we conducted, using TTClust,⁶² a hierarchical clustering analysis
based on the RMSD of the RNA, and the RMSD time series (see Figure S7 and S9 for
simulations of the whole RNA and different stem-loops, respectively) along with the RSMF
profile to investigate the flexibility. Based on this analysis (Figure S7 and Figure S8), two
regions - the P3 (nucleotide from 55 to 75) and P4 stem- loops (nucleotide from 90 to 102) —
responsible for a moderate RMSD value (between 5.7 Å and 7.2 Å) exhibit high flexibility and
undergo significant conformational changes in the 3 replicas. The high flexibility observed
in these regions corresponds to the large movement impacting the whole RNA structure,
which induces a high global change in the RMSD. For example, in the stem-loop P3 the
changes at the level of RMSD are related to the change of the orientation of the nucleotide
A55. On the contrary, the stem-loop P3a presents a low flexibility and a high stability
along all the replicas with an RMSD value around 1.0 Å and its standard deviation equal
to 0.3 Å. Its stability is also confirmed by three independent simulations conducted on this
isolated stem-loop (Figure S10). Our results clearly indicate that P3a stem loop is stable
and does not form any tertiary contacts with other parts of the molecule. Worth noting
P3a stem-loop is neither involved in direct contact with SAM, nor in the SAM-dependant
folding of the aptamer.³² Consequently SHAPE reactivity for nucleotides within P3a does
not vary upon SAM addition⁷⁶,⁷⁷ further confirming the overall stability of this stem-loop.
All together, these findings indicate that the P3a stem-loop is an ideal model system for
investigating SHAPE binding and reaction, as well as establishing a computational protocol
to study these processes.

Then, we conducted SHAPE experiments at three different BzCN concentrations and in
the absence of SAM to prevent any interference with the chemical reaction. At the level of the
P3a stem-loop, we observed notable differences between the three concentrations — 20, 50,
and 100 mM (see Figure S3 and Figure S3-4). Specifically, concentration-dependent effects
were obtained for A63 and A64, which exhibited distinct behaviors and greater variability
in their values. These findings underscore the importance of carefully selecting the probe
concentration in experiments to ensure meaningful comparisons with computational studies.
Consequently, we focused on results obtained at 20 mM for BzCN to minimize the existence of cooperative effects³¹ (see Figure 1C and D for the results). Based on that, G61 is identified as non-reactive, while A62 and A64 exhibit high reactivity, and A63 shows moderate reactivity (see Figure 5).

<figure><img src="image_7.png" /><figcaption>Figure 5: a). BzCN reactivity for the P3a loop represented in the 2D structure. b) Reactivity map of P3a loop for BzCN probe at 20mM. The horizontal lines indicate the level to reach moderate (yellow) or high (red) reactivity.</figcaption></figure>

Investigating BzCN binding using the model system P3a

To gain molecular insights into SHAPE reactivity, we applied a multiscale approach under
the hypothesis that the reaction requires the probe to successfully bind to the RNA first,
as also proposed by Bussi and co-workers.²¹ To do so, we investigated the binding of BzCN
using umbrella sampling, initially through all-atom MD simulations, followed by QM/MM
calculations for distances below 4 Å to properly account for bond formation. We conducted
all-atom US MD simulations using at first only the distance *d* between the reactive carbon
C* and the O2' of the ribose of the investigated nucleotide as a collective variable, as defined
above. To evaluate favorable or unfavorable binding states, we computed the potentials of
mean force (PMF) using WHAM⁵³ and then analyzed the profiles, in particular the free-
energy barriers of the binding states ($\Delta G_{bind}^{\ddagger}$) at a distance around 4 Å. Bound states are
defined here as structures for which the distance *d* < 4.0 Å, as proposed by Bussi and co-
workers in previous works.²¹

For the tetraloop under investigation, based on the high $\Delta G_{bind}^{\ddagger}$ binding free-energy
<figure><img src="image_8.png" /><figcaption>Figure 6: Potential of Mean Force (PMF) obtained by 1D umbrella sampling MD simulations using the distance between the O2' of the ribose and the reactive carbon C* as collective variable for G61 (top left, non reactive experimentally), A62 (top right, highly reactive), A63 (bottom left, moderately reactive) and A64 (bottom right, highly reactive). The vertical gray line represents the bound distance of 4 Å. For sake of clarity, the uncertainties are not shown since they are very small.</figcaption></figure>

barrier, the SHAPE probe, BzCN, cannot favorably bind the experimentally non-reactive nucleotide G61 probably due to its engagement in stacking interactions with the neighboring nucleobase instead at larger distances (Figure 6). A different scenario is obtained for the other three nucleotides (A62, A63 and A64) where a favorable state is observed based on the value of $\Delta G_{bind}^{\ddagger}$ (see Figure 6) ($G_{bind}^{\ddagger} < 1\,\text{kcal mol}^{-1}$). However, it is not possible to distinguish between the medium and the high reactive nucleotides based on these 1D PMF profiles, an important step to key features to explain SHAPE reaction. A similar free-energy profile is also obtained for A62 using 1D Replica Exchange US MD simulations (Figure S12), suggesting that the observed results are not due to insufficient sampling. In fact, by taking into account only the distance in the prediction of free energy profiles, we neglect the fact that for the same distance different relative positions of the probe with respect to nucleotide under investigation are possible. In other words, for the same binding distance the probe can
interact either with the ribose of interest or the neighboring nucleotides (see Figure 7). Due
to competition between hydrolysis (half-life 0.25 s at 37 °C) and the acylation reaction, this
could have a critical impact on the possibility of accommodating the probe in a favorable
position to undergo the acylation reaction in time since favorable binding should affect
subsequent nucleophilic attack. This aspect has been investigated via QM/MM calculations
for short distances, as presented below.

<figure><img src="image_9.png" /><figcaption>Figure 7: Exemplification to explain the different relative positions of the probe with respect to the nucleotide under investigation.</figcaption></figure>

To overcome these limitations and properly describe the orientation with respect to the
ribose, we defined a valence angle α based on the reactive carbon C*, the C4' and C2'
atoms of the ribose to challenge the probe recruitment with respect to the ribose without
biasing in advance the sampled nucleophilic attack angles as shown in Figure 2B. Before
computing the 2D US MD, we first assessed the impact of the probe on RNA structures,
specifically evaluating whether its presence induces significant structural changes. To do
so, we calculated the eRMSD<sup>70</sup> along the windows of the 1D US MD simulations. We
considered as reference the starting structure where the probe is not bound to the RNA.
Therefore, we obtained for each given pair of distance d and angle α in each window a value
of eRMSD. Figure 8 summarizes our results. When the non-reactive nucleotide G61 is under
investigation, we observe that the RNA molecule explores distinct structural ensembles,
as evidenced by the presence of multiple clusters (see "dynamic secondary structures"). Moreover, the eRMSD exceeds 0.7 in the same region, which is the proposed threshold to distinguish significantly similar structures, indicating substantial conformational variability for given values of *d* and *α*. (see Figure 8). In these simulations, RNA structure changes significantly when the probe approaches the ribose at distances smaller than 5 Å, explaining why probe binding is not a favorable event. A different scenario is obtained for the other three nucleotides. First, the binding of the probe does not dramatically affect the RNA structures, as shown by the eRMSD maps (Figure 8). Moreover, for highly reactive nucleotides (A62 and A64), BzCN can sample a larger variability of pairs of valence angles *α* and distance *d* with respect to the intermediate nucleotide A63.

Based on these results, we carried out 2D US simulations to gain a deeper understanding of the free-energy landscape with respect to the distance *d* and *α* that were used as collective variables. In particular, introducing the binding angle *α* enabled more effective sampling of the various possible probe orientations when binding to a given nucleotide. This is important because a single probe-nucleotide distance can correspond to multiple orientations, each potentially leading to distinct pre-reactive complexes with varying reactivity. Simulations were performed across 5-ns windows, covering the range from 854 to 1134 windows, with a total simulation time ranging from 4.27 µs to 5.67 µs per base and the results are summarized in Figure 9. For the non-reactive nucleotide G61 (Figure 9b and c), two separated basins are observed representing two different relative positions: outside the groove (lower energy conformations) and inside the groove (high energy conformations). However, as highlighted by the 1D PMF, the free-energy barriers obtained from the 2D PMF, $\Delta G_{bind}^{\ddagger}$, is large, therefore not compatible with a favorable bound state ($\Delta G_{bind,AB}^{\ddagger} \sim 3.2$ kcal mol$^{-1}$). In contrast, for reactive nucleotides, free-energy barriers (around 1 kcal mol$^{-1}$) separate two minima and favorable bound states are obtained (Figure 9a, d and e). Interestingly, for the highly reactive nucleotides, BzCN seems to bind easily to the RNA, as indicated by a flatter energy landscape. Additionally, the first binding minimum (indicated with the letter B in
<figure><img src="image_10.png" /><figcaption>Figure 8: Dynamic extended secondary structures representing base-pairing and base-stacking interactions following the Leontis-Westhof notation,⁷¹,⁷⁸ for the main clusters obtained for the 1D US MD simulations on the P3a tetraloop: a)-e) for G61, f) for A62, g) for A63 and h)-i) for A64. The percentage of the cluster appearance is indicated below each 2D representation. In the 2D structure, the nucleotide under investigation is colored in red. j), k), l) and m): eRMSD plotted based on the values of distance between the reactive carbon C* and the O2' of the ribose and the angle α obtained for G61, A62, A63 and A64 respectively. The value with an eRMSD larger than 0.7 are shown in gray.</figcaption></figure>
Figure 9a) occurs at lower values of the valence angle $\alpha$. This observation can allow us to speculate that a nucleotide to be reactive should not only be able to approach the ribose but also bind in a specific way. However, it is important to point out that the low free energy of the bound state does not presume whether the reaction barrier will be low or not.

## Investigating the impact of full structure on BzCN binding free energy surface

To confirm our findings and assess the reliability of the isolated stem-loop P3a as a suitable model system, we carried out 1D and 2D US MD simulations on the same nucleotides using the previously described protocol while considering the entire structure. The 1D PMFs differ from those obtained for the isolated stem loop P3a, since for large distances there are contributions of the rest of the structure for the alternative binding modes (Figure S13). In other words, for the same distance $d$ different binding angles $\alpha$ can be explored, as explained above. Although we can still distinguish between nonreactive and reactive nucleotides based on favorable and unfavorable binding determined by the free-energy barrier values $\Delta G_{bind}^{\ddagger}$, the analysis can become more complex because favorable bound structures can be related to distant RNA regions interacting with the probe other than the investigated nucleotide.

Subsequently, we compared the 2D PMF profiles obtained for the loop in the isolated structure and in the entire molecule of the stem loop P3a (Figure S14). First, this comparison highlights that the free energy surfaces offer a clearer understanding and facilitate more meaningful comparisons. For the non-reactive nucleotide G61, the 2D PMF profile obtained by considering the entire structure differs from that of the isolated stem-loop. In this case, the free energy surface reveals that, due to the presence of the remaining structure, only a portion of the free energy landscape (corresponding to large values of $\alpha$ also described as the high free energy region in Figure 9b) falls below $3.25 \text{ kcal mol}^{-1}$ highlighting the importance of comparing the free-energy profiles obtained from the whole structure and the isolated model to understand its limitation and the critical role of neighboring nucleotide
<figure><img src="image_11.png" /><figcaption>Figure 9: Potential of Mean Force (PMF) surface obtained by 2D umbrella sampling MD simulations using the distance *d* between the O2' of the ribose and the reactive carbon C* and the valence angle α as collective variables. a): 2D PMF obtained for A62 with a representative structures for the three minima; b) and c): 2D PMF obtained for G61; d): 2D PMF obtained for A63. e): 2D PMF obtained for A64.</figcaption></figure>
considered in the analysis. Also for the whole structure, a free-energy barrier higher than
1.0 kcal mol⁻¹ is observed further indicating that binding is unfavorable,²¹ as drawn from
the simulations for the isolated stem-loop. In contrast, for the reactive nucleotides, the 2D
PMF profiles for the high and medium reactive nucleotides obtained for the full structure
are similar to those obtained with the isolated the stem-loop (Figure 9), showing again a
favorable binding at *d* = 4 Å. These results confirm that the rest of the structure does not
significantly affect probe binding, enabling us to study the loop by considering only the
isolated P3a stem-loop at a considerably reduced computational cost. Furthermore, our
findings suggest that comparing the free energy surfaces as a function of the distance *d* and
the valence angle α should provide valuable insights into the impact of 3D contacts on the
binding process, making these profiles particularly well-suited for this analysis and, more
generally, to describe the binding of SHAPE probes. Based on all that, in the following, we
will refer only to the analysis on the biased simulations conducted on the isolated the P3a
stem-loop.

RNA local structural characterization in the presence of the probe

To analyze the binding in more detail in the biased simulations of the P3a stem-loop, we first
computed the puckering for each snapshot and then examined the puckering conformations
at the bound state (see Figure S15). Although the most common puckering is the C3'-endo,
the ribose can adopt various conformations; notably, highly reactive nucleotides exhibit
increased variability, highlighting the inherent flexibility of the ribose. These findings align
with previous studies,²¹,²⁸ which demonstrate that SHAPE reactivity is closely related to
ribose flexibility, an essential element for the chemical reaction.

Furthermore, for each nucleotide under investigation and its neighboring residues, we
evaluated whether BzCN was stacked. This was determined based on two criteria: (i) the
distance between the centers of mass (less than 5 Å) and (ii) the stacking angle α<sub>stack</sub> defined
as the angle between the perpendicular vectors to the planes (less than 45°), as outlined in
the Materials and Methods section. To further characterize the stacking, particularly the
relative position of the probe with respect to the ribose and nucleobase, we also computed the
in-plane displacement *d*<sub>xy</sub> if the probe is stacked. Additionally, we calculated the reweighted
3D PMF using the distance *d*, the bond angle α and *d*<sub>xy</sub> as CVs for the reactive nucleotides in
the tetraloop. For the highly reactive nucleotide A62 (see Figure 10), stacking was strongly
observed not only with its ribose and nucleobase but also with the adjacent ribose of A63.
At this stage, both basins of α appear favorable for probe stacking onto A62. However, only
basin B (low values of α) allows additional stacking with the adjacent ribose (A63). This
picture becomes clearer when A63 is analyzed independently: the probe stacks exclusively
onto the ribose of A63, but within a narrow α range (20° to 60°). Within this range, and
at higher α values, the probe also stacks onto the ribose of A62. This observation aligns
with experimental reactivity data for A62 and A63— when the probe is positioned nearby, it
has access to both nucleotides, introducing potential competition. The narrower α window
for A63 suggests a more constrained stacking configuration, which might explain its lower
observed experimental reactivity compared to A62 (Figure S16 and S17). In contrast, for
A64, favorable α values correspond to stacking interactions either with the ribose or the
nucleobase. However, no stacking with the adjacent nucleotide A63 is observed, indicating
a lack of competition, unlike the case of A62/A63 (Figure S18 and S19).

Although our all-atom US simulations allow us to investigate the binding event, the
chemical reaction can only proceed if specific values of the attack angles β₁ (Bürgi–Dunitz
angle) and β₂ (analogous to the Flippin–Lodge angle) are adopted. To do so, we computed
the 3D PMF reweighted on the distance *d*, the bond angle α and the attack angle β₁ or β₂
and we investigated whether there is a relationship between the attack angles sampled at the
binding distance of 4 Å and the bond angle α (see Figure 11). At the binding site, the attack
angle β₁ exhibits a broad distribution, indicating that nearly all values of the attack angle
can be favorably sampled. More in detail for A62, no differences were observed regardless of
the bond angle α. In contrast, for A63 and A64, an energetic preference was observed with
<figure><img src="image_12.png" /><figcaption>Figure 10: Stacking values obtained for each pair of distance *d* and *α* along the 2D US MD simulations on A62. The stacking with respect to the ribose (left) and the nucleobase (right) belonging to A62 (top) and to 63 (bottom) are shown. The reweighted 3D PMF as a function of the distance, the bond angle *α* and the displacement *d*<sub>xy</sub> if the stacking is present is shown at binding distance. A vertical dash line separates the results on the ribose (left) and the nucleobase (right).</figcaption></figure>

low values of the bond angle α. To access complementary pieces of information, we analyzed
the 3D PMF for the attack angle β₂ at the binding distance of 4 Å. As for the attack angle
β₁, the attack angle β₂ can be favorably sampled over a wide range of values for the three
nucleotides under investigation. Four energetic basins are observed and commonly separated
by low free-energy barriers. The basins in terms of the attack angle β₂ are larger for the
minima with low value of the bond angle α (around 20° for A62 and for A64, and 30° for A63).
These findings show that when the binding is favorable, there is no clear preference between
optimum and non-optimum values. Investigating the dynamics of these attack angles and
their possible convergence toward the values that allow the chemical reaction to occur would
require conducting the simulation below 4 Å, which is not possible when using only an MM
description.
<figure><img src="image_13.png" /></figure>

Figure 11: Reweighted 3D PMF as a function of the distance, the bond angle $\alpha$ and the attack angle $\beta_i$ ($i = 1, 2$) projected on the 2D surface for $\alpha$ and $\beta_i$ for $d$ between 3.9 and 4 Å for the nucleotide A62 (A and D for the attack angle $\beta_1$ and $\beta_2$ respectively), for the nucleotide A63 (B and E for the attack angle $\beta_1$ and $\beta_2$ respectively) and for the nucleotide A64 (C and F for the attack angle $\beta_1$ and $\beta_2$ respectively).

Investigating binding scenarios outcome on C\*O2' bond formation using QM/MM US calculations

Before investigating in details the shorter range of the distance $d$ using QM/MM, we performed NEB-TS calculations on a single nucleotide based on the work by Jash and Kool¹⁵ (see SI for more details) by taking into account the two possible protonation states of the hydroxyl group (see Figure 2c and d). These calculations highlight that when the hydroxyl group is protonated a favorable reaction can take place where at the transition state four distances ($d_1$ between the reactive carbon and the O2' of the hydroxyl group, $d_2$ between the O2' of the hydroxyl group and the hydrogen of the hydroxyl group, $d_3$ between the reactive
carbon and the carbon of the cyanide group and $d_4$ between the hydrogen of the hydroxyl group and the carbon of the cyanide group) are involved (see Figure S11). In contrast, for the oxyanion case, the optimized structures of the reactants show phosphate group cleavage occurring without the involvement of the SHAPE probe. Additionally, the transition state structure reveals significant disruption of the overall geometry. Only the product remains unaffected, most likely because the oxyanion function is no longer present (see Figure S11). Although these calculations do not capture the full dynamics of the system or account for multiple transition states, they do confirm the relevance of the chosen collective variable. Moreover, they highlight the importance of including the protonated state of the hydroxyl group in our QM/MM calculations. This allows the system to adjust its protonation state during the simulation, addressing the instability of the oxyanion form that arises when the hydroxyl group is deprotonated by the buffer—especially in the absence of prior knowledge about the hydrogen's location.

To further investigate the shorter range of distances $d$ and evaluate the formation of the bond between the reactive C\* of BzCN and the O2' atom of the ribose in our case system, we switched to a QM/MM description to perform 1D QM/MM US simulations on nucleotide A62. To study the impact of binding on the formation of the C\*O2' bond while keeping computational costs manageable, we selected three binding states as starting points. Each state is characterized by a distinct initial value of the bond angle $\alpha$: $\alpha_i30$, $\alpha_i60$ and $\alpha_i80$. When we evaluated the evolution of the sampled values of the bond angle $\alpha$ across the three QM/MM US simulations, we observed that, without applying any bias on $\alpha$, the values of the binding angle immediately decreased in the cases of $\alpha_i60$ and $\alpha_i80$ simulations, approaching the lower values corresponding to the minimum B as determined by the 2D PMF (Figure 9). To be able to evaluate the different binding states, we extended the additional bias on the bond angle $\alpha$ between 4 Å and 2.5 Å during the minimizing and heating steps. As soon as we remove the additional bias on $\alpha$, the angle values decrease again, approaching those corresponding to the minimum B especially when getting close to
2.5 Å, with larger variability around 4 Å. This evolution closely resembles the one observed in the QM/MM simulation initiated from αᵢ30 (Figure S20). These results suggest that a more complex relationship is present between the binding states and the occurrence of the chemical reaction, explaining why in MM US we were unable to fully capture the relationship between the bond angle α and the attack angles β₁ and β₂.

As in the all-atom US MD simulations, we computed the PMF profiles from the three simulations which exhibit strong similarities (Figure 12). In every case, the free-energy begins to increase around 2.5 Å where we propose that the reaction starts. Notably, around 1.4–1.5 Å, corresponding to the value for the C\*O2' bond at equilibrium, the αᵢ60 simulation displays a distinct local energy minimum around 1.4 Å (Figure 12). In contrast, the free energy profiles for αᵢ30 and αᵢ80 continue to increase, suggesting the presence of different sampled transition states.

To better characterize the chemical changes, we computed three key distances with respect to the distance *d* (the 1D CV). In BzCN, the evolution of the distances between C\* and O1 (*d*<sub>C\*O1</sub>) and between C\* and C8 (*d*<sub>C\*O1</sub>) is similar across the three QM/MM US simulations: the distance *d*<sub>C\*O1</sub> increases suggesting the reduction of the carbonyl (C=O) to a single C-O bond. In contrast, the evolution of the distance between the O2' atom of the ribose and H2' in the αᵢ60 simulation differs from that observed under the other conditions, suggesting that this variation may be responsible for the differences observed in the PMF profile. In contrast to the αᵢ30° and αᵢ80° simulations, the atom O2' can be deprotonated when the probe approaches leading to an oxyanion, and it was constantly deprotonated by the O1 atom of BzCN at CV values below 1.5 Å, corresponding to a local minimum associated with hemiacetal formation. Between 1.5 and 1.6 Å, partial deprotonation of the O2' atom was observed, facilitated by the O3 atom of the phosphodiester group and N3 of the nucleobase of the nucleotide A62 (Figure S22). In fact, consistent oxyanion formation is observed only in the αᵢ60 trajectory. This process is facilitated by proton hopping, beginning with N3, O3', and OP2, and ultimately stabilizing on O1 of the probe.
Analyzing the structures with values of *d* at the equilibrium of bond formed by the oxygen of the hydroxyl (O2') and the reactive carbon (around 1.4 Å and 1.5 Å) we can notice that for the simulations α<sub>i</sub>30 and α<sub>i</sub>80 the 2'-hydroxyl group remains protonated. In these cases, the electron density appears deficient near the O2' atom (3 covalent bonds), while O1 adopts a partial anionic character due to reduction of the carbonyl group of BzCN, located adjacent to the nitrile moiety. This suggests that resonance stabilization (mesomeric effect) is not fully established in these cases. These observations may be attributed to the nature of the collective variable (CV), which does not explicitly capture the entire reaction pathway. However, the observed differences likely reflect the varying polarization effects introduced by the probe. To better characterize these effects and assess both similarities and differences between simulations, we first performed a detailed geometrical analysis.

<figure><img src="image_14.png" /><figcaption>Figure 12: Potential of Mean Force (PMF) obtained by 1D QM/MM umbrella sampling using the distance *d* between the O2' of the ribose and the reactive carbon C* as collective variable for each different starting poses α<sub>i</sub>30 (red), α<sub>i</sub>60 (blue) and α<sub>i</sub>80 (orange) a): 1D PMF obtained for A62 with a representative structures; b) pose at 4 Å for α<sub>i</sub>60; c): pose at 4 Å for α<sub>i</sub>30; d): pose at 4 Å for α<sub>i</sub>80.</figcaption></figure>

To characterize the relative orientation of the probe BzCN with respect to the nucleobases
of the target nucleotide A62 and the adjacent nucleotide A63, we calculated the distribution
of $\kappa_i$ (see the definition in Material and Methods) as a function of the CV. For the $\alpha_i$ 80 simulation, the probe remains on the opposite side of the nucleobases both when approaching the ribose and after the C\*O2' bond formation (Figure S21b). The opposite scenario was observed in the $\alpha_i$ 60 simulation where the probe remained consistently positioned between the two nucleobases of the target nucleotide (A62) and the adjacent A63 (Figure S21c). Interestingly, in the $\alpha_i$ 30 simulations, the probe was initially positioned on the opposite side before shifting its relative position to lie between the two nucleobases as the distance $d$ decreases (Figure S21a). Additionally, as $d$ decreased and C\*O2' bond formed, a small shift of $\kappa$ of the nucleobase of A62 was observed suggesting changes in dynamics and orientation with respect to the backbone due to a restricted area resulting from the presence of the probe. It is also at these conditions, that a distinct protonated state was observed (Figure S21a) where N3 (of the nucleobase) captured the proton. In the case of $\alpha_i$ 80, a protonated state was also observed at an unphysically short C\*O2' bond length (1.3 Å), suggesting an electronically unstable configuration. Therefore, it was not considered further.

To investigate the mechanistic differences that follow from these orientations, we analyzed the nucleophilic attack geometry. Since the reactive carbon C\*, is a prochiral center, we chose to characterize the attack based on the Re/Si face distinction using the $\beta_2$ angle, also used to analyze the all-atom US MD simulations, as described in the Materials and Methods section. The attack angle $\beta_2$ varies without a specific preference when the distance $d$ between the reactive carbon C\* and the O2' atom is between 3 and 4 Å, as already highlighted by the all-atom simulations. On the contrary, when $d < 3$ Å, the values of the attack angle $\beta_2$ converge within a specific range allowing to distinguish between an attack on the Si or Re face. In the $\alpha_i$ 30 and $\alpha_i$ 80 simulations, the nucleophilic attack occurs on the Re face, whereas in the $\alpha_i$ 60 simulation, the attack is directed at the Si face (Figure 13).

To get a better understanding and delve deeper into RNA local geometric changes associated with the nucleophilic attack in each of the three $\alpha_i$ US QM/MM simulations, we analyzed the distribution of the angles $\xi$ and $\tau$, along with the nucleobase-nucleobase and
<figure><img src="image_15.png" /><figcaption>Figure 13: The evolution of the flipping Lodge angle $β_2$ across the three US $α_i$30 (red), $α_i$60 (blue) and $α_i$80 indicating whether the nucleophilic attack occurs on the Re or Si face.</figcaption></figure>

phosphate-phosphate distances as defined in the section Material and Methods. For the $α_i$30 simulation, at larger distances $d$, the nucleobases are close one to each other as also shown with an inter-base distance ($C6_{A62}C6_{A63}$) of approximately 7.9 Å. This configuration restricts access to the ribose, as indicated by $ξ$ values near 25° (Figure 14a). However, as the CV decreases to 2–3 Å, the distribution becomes bimodal with $ξ$ shifting to approximately -50° and the interbase distance increasing to around 10 Å. This suggests that the local RNA structure begins to open before eventually stabilizing into a monomodal distribution at $ξ$ around -50° (Figure 14a), which corresponds to a fully open structure with an inter-base distance of 17.1 Å (Table S1). This openness indicates that the RNA accommodates the probe at short distances, thus altering the local accessibility, a feature difficult to observe without the explicit presence of the probe despite the well-documented importance of the fluctuations in the distance between consecutive nucleobases.²²,²⁸ Moreover, this nucleobase/probe opening is accompanied by a compensatory adjustment in the backbone. The angle $τ$, which characterizes backbone rigidity, exhibits an opposite trend compared to $κ$ it starts at a relatively open value (around -25°) at approximately 4 Å before shifting toward 0°. In the $α_i$60
simulation, with respect to the αᵢ₃₀ simulation, the RNA begins to exhibit local conformational changes even when the probe is distant (d > 3 Å) as evidenced by a τ distribution already near 0° and a broad ξ distribution spanning from closed to open states (Figure S23a). In contrast, the αᵢ₈₀ simulation shows no significant shifts, except for a slight transition from an initially open conformation toward a tighter base-base interaction, particularly during the formation of the C\*O2' bond probably due to the presence of a restricted area after the bond formation (Figure S23b).

Finally, to have the complete picture, we can also combine the pieces of information obtained by the distribution of κᵢ on the relative orientation of the probe BzCN with respect to the target nucleotide A62 and the adjacent nucleotide A63, and the local geometrical adjustments of RNA via the angles ξ and τ. Analysis of the κᵢ reveals that, in the αᵢ₆₀ simulation—where the probe initially resides on the same side as the nucleobases—rapid local structural adjustments occur. In contrast, the αᵢ₃₀ simulation exhibits dynamic changes, as evidenced by shifts in the distributions of the τ and ξ angles. These findings highlight the importance of explicitly modeling the probe when investigating SHAPE reactivity, since local geometric alterations may arise from probe interactions rather than from pre-existing RNA conformational preferences.

## Effect of protonation and deprotonation of O2' hydroxyl in the first step of the acylation reaction

As highlighted above, the different QM/MM simulations highlight different protonation / deprotonation states based on the starting structures. In fact only a consistent oxyanion formation is observed only in the αᵢ₆₀ when it is not observed for the simulations αᵢ₃₀ and in the αᵢ₈₀. This process is facilitated by proton hopping, beginning with N3, O3', and OP2, and ultimately stabilizing on O1 of the probe. Moreover, analyzing the geometry of the nucleophilic attack allows us to investigate not only the mechanistic differences arising from these orientations but also their relationship with the protonation state. Analyzing the
<figure><img src="image_16.png" /><figcaption>Figure 14: Density distribution of ξ (orange) and τ (blue) as a function of the collective coordinate *d* obtained for the α<sub>i</sub>30 simulation. The angle ξ is defined by the atoms P63, C1 of nucleobase A62, and C1 of nucleobase A63, while τ is measured as the angle between the vectors P61-P63 and P63-P64. (a) Distribution for *d* > 2 Å, representing configurations before adduct formation; (b) Distribution for *d* < 2 Å, capturing structural adjustments near and after bond formation.</figcaption></figure>

Bürgi–Dunitz angle (ideally ~ 107) using the β₁ definition, all three simulations sampled
geometries within the ideal range. However, only the Si-face trajectory, αᵢ60, resulted in
deprotonation of the 2'-hydroxyl group and stabilization of the hemiacetal intermediate.
Further analysis of interactions with the major proton acceptor, N3—specifically at distances
greater than 2.5 Å—showed that the frequency of hydrogen bonding between N3 and the
2'-hydrogen was significantly higher in the αᵢ60 case, occurring more than twice as often as
in αᵢ30 and αᵢ80 (11%, 3.4%, and 4.9%, respectively).

Another key aspect for the protonation state is the local environment that can result
from the relative position of the BzCN with respect to the nucleobase. In αᵢ₆₀, the BzCN
probe intercalates between the two nucleobases at large distance (Figure S21), forming a
tightly stacked environment. In contrast, in αᵢ₈₀ BzCN is mostly not stacked with the
neighboring nucleobases, resulting in a weak polarizability effect on the nucleobase side where
the hydroxyl group is typically oriented, while in α<sub>i</sub>30 the probe intercalation occurs only
mid-reaction, which may explain the presence of some deprotonated structures, although
these were outnumbered by protonated ribose forms. These findings suggest that early
intercalation of the probe enhances local polarization (shifting the effective pKa) which
facilitates proton transfer and stabilizes the oxyanion. The fact that O1 (from the probe)
acts as the final proton acceptor highlights the probe's direct role in shaping the reactive
environment. Furthermore, the oxyanion state appears unstable without favorable structural
and electronic conditions and often reverts to the neutral hydroxyl form. When the probe is
distant, the environment becomes rigid, limiting access to major proton acceptors. Only O3
remains accessible, and proton transfer via this route typically coincides with phosphodiester
bond cleavage. For instance, N3 of the nucleobase is unavailable as a proton acceptor due
to persistent stacking with A63, as shown by dynamic secondary structure analysis from
all-atom 1D umbrella sampling (Figure 8f). This restriction is relieved when the probe
approaches, as reflected in shifts in the τ angle and C6–C6 distance (see Table S2).

To further investigate the equilibrium between the protonated and deprotonated states
as suggested by Kool and Jash¹⁵ and rule out the importance or not of the oxyanion in the
first step of the SHAPE reaction, first in the case of αᵢ60, where different deprotonation
states are already observed prior to the first reaction step using only the collective variable
d, we assessed the stability of the protonated N3 and the 2'-oxyanion state observed at d
around 3.7 Å. By extending the simulation window around this state (d = 3.7 Å) for 400
ps, we found that the O2' oxyanion was stabilized by 3 (47.7%) to 4 (34.2%) hydrogen
bonds with nearby water molecules (Figure S24). These interactions exhibited a prominent
peak at r = 1.55 Å in the radial distribution function g(r). This first-shell peak at 1.55
Å—shorter than the typical hydrogen bond length in the TIP3P water model and closer to
proton transfer distances—suggests that the system was at the threshold of proton transfer,
a process likely limited by the classical treatment of water molecules (Figure S25).

To enhance the sampling of the oxyanion species, aware of the slow equilibrium between
the protonated/deprotonated states of the 2'-hydroxyl group, we also considered a new collective variable $CV_2$ (see eq. (1)) in the presence and absence of four water molecules (potential proton acceptors) in the QM region with semi-harmonic potential to restrain the oxygen mobility for the latter. The presence of water molecules allowed us to investigate water-mediated proton transfer, especially since in $\alpha_i30$ and $\alpha_i80$, hydrogen bonds with the major proton acceptor (N3) occurred at approximately half the frequency observed in $\alpha_i60$. For $\alpha_i30^\circ$, we analyzed the range of collective variable values $d$ between 2.0 and 1.4 Å, along within the new extended collective variable $CV_2$ (Eq. (1)) (see Figure S26a). Interestingly, starting slightly beyond the equilibrium distance $d$ (around 2 Å) corresponding to a value of $CV_2$ around 0.5 Å, a deprotonation of 2'-hydroxyl group occurs. However, nearby water molecules quickly acted as proton donors, restoring the protonated state of the 2'-hydroxyl (Figure S26b). At the equilibrium distance $d$ (between 1.5 and 1.4 Å), the buffer participates in capturing the proton and with the O1 atom of BzCN protonated via the buffer also, an hemiacetal was formed as well as in $\alpha_i60$ case suggesting that the formation of oxyanion prior to binding and acylation initiation is not required for the reaction to start and that the probe presence at short range brings enough polarization to deprotonate the 2'-hydroxyl group (Figure S26c). Depending on the local environment, this deprotonation is facilitated by nearby acceptors or by the buffer. These results allow us to propose that the differences in the PMF profiles of $\alpha_i60$ and $\alpha_i30$ (Figure 12) could be due to the deprotonation state. These findings, together with those discussed above, suggest that deprotonation of the 2'-hydroxyl group occurs only at short C*–O2' distances (below 2 Å and before equilibrium is reached). This indicates that the presence of the probe is essential for inducing local polarization and promoting deprotonation.
Conclusions

In this work, we obtained an unprecedented molecular description of the interactions between
the SHAPE probe BzCN and the RNA by conducting a multiscale study that combines stan-
dard all-atom molecular dynamics simulations, all-atom umbrella sampling simulations and
QM/MM US simulations, along with new SHAPE experiments. We investigated the *Bacillus
subtilis* yitJ S box/SAM-I riboswitch RNA, in particular the P3a stem-loop caracterized by
the well-known GAAA motif. Because it contains nucleotides with low, medium and high
reactivity, P3a was an ideal case system. Standard all-atom MD simulations on the whole
structure and the isolated stem-loop highlighted its stability and the lack of interactions with
the rest of the structure, a result that is in agreement with SHAPE experimental data.

Based on previous work conducted by Bussi and co-workers where they suggested that
binding has an impact on overall reaction barrier,²¹ we set-up a multiscale approach to
link the binding to the first step of the SHAPE acylation reaction. To investigate the
binding process and the approach of the probe prior to the formation of the chemical bond
between the reactive carbon C* of the probe and the O2' atom of the ribose, we highlighted
the necessity of combining both all-atom and QM/MM descriptions. We conducted 1D
then 2D US MD simulations using, for the first case, the distance *d* between the reactive
carbon C* of the probe and the O2' atom of the ribose then by combining *d* and the bond
angle α, that accounts for the relative position of the probe with respect to the ribose,
as collective variables. For the non-reactive nucleotide using both strategies no favorable
binding is observed. On the contrary for the reactive nucleotides, our study highlighted that
the introduction of the bond angle α as CV allows to identify specific favorable bound states
at the binding distance *d* = 4 Å. The analysis of the free-energy surface also allowed us to
better characterize key interactions between the probe and the RNA by computing different
geometrical descriptors. The probe is recruited mainly by stacking interactions with either
the nucleobase or the ribose or both of them. Moreover, at the binding distance, a large
variability of ribose puckering is observed suggesting the importance of ribose flexibility in
the acylation reaction.²¹⁻²³ In order to study larger systems, a possible strategy would be to compute replica exchange 1D US sampling simulations using either the temperature or the hamiltonian and reweigh the simulations based on the two CVs to anticipate the energetic barriers and conduct if needed 2D US simulations on selected nucleotides.

To study the first step of the reaction, the analysis of attack angles in all-atom and QM/MM US simulations revealed that, at larger distances *d* (*d* > 3 Å) there is no clear preference for the attack angles β₁ and β₂, that can be largely sampled at the bound minima. However, QM/MM US simulations conducted on three different starting α angles showed that a preference emerges when *d* < 3 Å. Additionally, these simulations demonstrated that, despite starting from initial α angles different from the minimum, the α angle tends to converge toward smaller α values. This suggests an indirect relationship between the bound α angle, the attack angles and the likelihood of the reaction occurring. The free-energy profiles obtained from the three simulations show slight differences, particularly in terms of the transition state flowing also the differences observed with β₂. To better assess these variations, we introduced specific descriptors to capture local changes on the RNA level along with the relative position of the probe with respect to the target nucleotide and its neighboring one, implicitly accounting for the dynamics of the ribose accessibility. Moreover, we investigated the protonation dynamics of the ribose 2'-hydroxyl group and its role in the reaction, focusing on the influence of both the probe and the buffer. Our results demonstrate that the presence of the probe is critical for inducing local polarization and enabling deprotonation. Importantly, this suggests that oxyanion formation does not need to precede probe binding—which could otherwise lead to RNA cleavage, as indicated by our QM calculations. Instead, the probe itself provides the necessary polarization that could initiate deprotonation highlighting the possibility of different deprotonation/protonation pathways at short distances between the reactive carbon and the O2' of the ribose.

Based on all these findings, we hypothesize that binding scenarios result in subsequent different reaction paths. A low binding free-energy barrier reflect a favorable proximity
between the probe and the nucleotide where RNA can also undergo local structural changes.
Moreover, these conformations allow larger sampling of attack angles facilitating the reaction
at shorter distance. All these events prior to the SHAPE acylation reaction are in competition
with the probe hydrolysis. Therefore, it is essential for the chemical reagent to adopt the
most favorable bound state that minimizes the need for structural readjustments prior to
the reaction. Therefore, our work paves the way to studying the complete SHAPE reaction
mechanism —an essential step toward uncovering the complex chemical origins of SHAPE
data. In particular our results allow to better understand the unexpected SHAPE reactivity
profile of the overrepresented GAAA tetraloop. Our results are instrumental for the future
utilization of SHAPE data for 3D RNA structure prediction.

Acknowledgement

This work was supported by ANR-22-CE29-0009 and by ANR-21-CE45-0034603 (INSANNE).
Work at CiTCoM is also funded by CNRS and Université Paris Cité. ABA is a recipient of a
MESRI (Ministere français de l'Enseignement Supérieur, de la Recherche et de l'Innovation)
fellowship awarded by the MTCI doctoral school of Université Paris Cité. PH was funded
by ANR-20-CE12-0026-02 (PRERIBO60S). CH was funded by ANR-22-CE29-0009 (MA-
GENTA). The Grand Équipement National De Calcul Intensif (GENCI) is acknowledged
for funding this research with grant A0120711431, A0140714119 and A0160714119. We also
thank Université Paris Cité being part of the Couperin Consortium for the read-and-green
agreement.

Data and Software Availability

The simulations are available at Zenodo 10.5281/zenodo.15176027, as well as part of the
scripts used in the analysis https://github.com/elisafrezza/SHAPE_MD_BzCN
## Contributions

Writing original manuscript: C.H., A.B.A., P.H., E.F. Revised manuscript: C.H., A.B.A., P.H., M.E.Q., B.S., D.L., E.F. All-atom MD simulations: C.H. NEB-TS calculations: C.H. QM/MM simulations: A.B.A. Data analysis: C.H., A.B.A., P.H., E.F. SHAPE experiments: P.H. Supervision: E.F., D.L., B.S., M.E.Q., Conceptualization: E.F. Funding Acquisition: B.S. and E.F.

## Supporting Information Available

The Supporting Information is available free of charge on the ACS Publications website at: Computational details on Force field parametrization and QM/MM calculations; Details on SHAPE experiments: RNA preparation, SHAPE chemical probing. Additional results on SHAPE BzCN P3a tetraloop reactivity. Additional results from standard all-atom MD simulations: Comparison correlation between SHAPE reactivity for BzCN and ASA and CV; Stability of *Bacillus subtilis* yitJ S box/SAM-I riboswitch (RMSD and RMSF plots on the whole structure and the P3a tetra-loop). 1D Umbrella Sampling results on the entire structure. 2D Umbrella Sampling results on the entire structure. Structural analysis on 2D Umbrella Sampling MD simulations on the isolated P3a stem-loop: puckering at the binding distance, stacking and displacement. Structural analysis on 1D QM/MM Umbrella Sampling simulations on the isolated P3a stem-loop: 2D distribution of the parameter $\kappa_i$; time-series of key distances; Distribution of the $\xi$ and $\tau$ parameters; Table for the inter-phosphate distances and the distance between the nucleobase A62 and A63.

## References

(1) Sharp, P. A. The centrality of RNA. *Cell* **2009**, *136*, 577–580.
(2) Morris, K. V.; Mattick, J. S. The rise of regulatory RNA. *Nat. Rev. Genet.* **2014**, *15*, 423–437.

(3) Ganser, L. R.; Kelly, M. L.; Herschlag, D.; Al-Hashimi, H. M. The roles of structural dynamics in the cellular functions of RNAs. *Nat. Rev. Mol. Cell Biol.* **2019**, *20*, 474–489.

(4) Incarnato, D.; Oliviero, S. The RNA epistructurome:: Uncovering RNA function by studying structure and post-transcriptional modifications. *Trends in Biotechnology* **2017**, *35*, 318–333.

(5) Lilley, D. M. J. How RNA acts as a nuclease: some mechanistic comparisons in the nucleolytic ribozymes. *Biochem. Soc. Trans.* **2017**, *45*, 683–691.

(6) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. The protein data bank. *Nucleic acids research* **2000**, *28*, 235–242.

(7) Cao, X.; Zhang, Y.; Ding, Y.; Wan, Y. Identification of RNA structures and their roles in RNA functions. *Nature Reviews Molecular Cell Biology* **2024**, 1–18.

(8) Weeks, K. M. Advances in RNA structure analysis by chemical probing. *Curr. Opin. Struct. Biol.* **2010**, *20*, 295–304.

(9) Hajdin, C. E.; Bellaousov, S.; Huggins, W.; Leonard, C. W.; Mathews, D. H.; Weeks, K. M. Accurate SHAPE-directed RNA secondary structure modeling, including pseudoknots. *Proc. Natl. Acad. Sci. U. S. A.* **2013**, *110*, 5498–5503.

(10) Mitchell, D., 3rd; Assmann, S. M.; Bevilacqua, P. C. Probing RNA structure in vivo. *Curr. Opin. Struct. Biol.* **2019**, *59*, 151–158.

(11) Douds, C. A.; Babitzke, P.; Bevilacqua, P. C. A new reagent for in vivo structure
probing of RNA G and U residues that improves RNA structure prediction alone and
combined with DMS. **RNA 2024**, *30*, 901–919.

(12) Spasic, A.; Assmann, S. M.; Bevilacqua, P. C.; Mathews, D. H. Modeling RNA secondary structure folding ensembles using SHAPE mapping data. *Nucleic acids research* **2018**, *46*, 314–323.

(13) Hardouin, P.; Lyonnet du Moutier, F.-X.; Sargueil, B. RNA Design: Methods and Protocols; Springer, 2024; pp 177–191.

(14) Mortimer, S. A.; Weeks, K. M. Time-resolved RNA SHAPE chemistry: quantitative RNA structure analysis in one-second snapshots and at single-nucleotide resolution. *Nature protocols* **2009**, *4*, 1413–1421.

(15) Jash, B.; Kool, E. T. Conjugation of RNA via 2'-OH acylation: Mechanisms determining nucleotide reactivity. *Chemical Communications* **2022**, *58*, 3693–3696.

(16) Mortimer, S. A.; Weeks, K. M. A fast-acting reagent for accurate analysis of RNA secondary and tertiary structure by SHAPE chemistry. *Journal of the American Chemical Society* **2007**, *129*, 4144–4145.

(17) McGinnis, J. L.; Dunkle, J. A.; Cate, J. H.; Weeks, K. M. The mechanisms of RNA SHAPE chemistry. *Journal of the American Chemical Society* **2012**, *134*, 6617–6624.

(18) Saaidi, A.; Allouche, D.; Regnier, M.; Sargueil, B.; Ponty, Y. IPANEMAP: integrative probing analysis of nucleic acids empowered by multiple accessibility profiles. *Nucleic acids research* **2020**, *48*, 8276–8289.

(19) Calonaci, N.; Jones, A.; Cuturello, F.; Sattler, M.; Bussi, G. Machine learning a model for RNA structure prediction. *NAR genomics and bioinformatics* **2020**, *2*, lqaa090.

(20) Washietl, S.; Hofacker, I. L.; Stadler, P. F.; Kellis, M. RNA folding with soft constraints:
reconciliation of probing data and thermodynamic secondary structure prediction. *Nucleic acids research* **2012**, *40*, 4261–4272.

(21) Mlýnský, V.; Bussi, G. Molecular dynamics simulations reveal an interplay between SHAPE reagent binding and RNA flexibility. *The journal of physical chemistry letters* **2018**, *9*, 313–318.

(22) Pinamonti, G.; Bottaro, S.; Micheletti, C.; Bussi, G. Elastic network models for RNA: a comparative assessment with molecular dynamics and SHAPE experiments. *Nucleic acids research* **2015**, *43*, 7260–7269.

(23) De Bisschop, G.; Allouche, D.; Frezza, E.; Masquida, B.; Ponty, Y.; Will, S.; Sargueil, B. Progress toward SHAPE constrained computational prediction of tertiary interactions in RNA structure. *Non-coding RNA* **2021**, *7*, 71.

(24) Thaplyal, P.; Bevilacqua, P. C. *Methods in enzymology*; Elsevier, 2014; Vol. 549; pp 189–219.

(25) Sen, S.; Pal, U.; Maiti, N. C. p K a Determination of d-Ribose by Raman Spectroscopy. *The Journal of Physical Chemistry B* **2014**, *118*, 909–914.

(26) Lippert, B. Ligand-pKa Shifts through Metals: Potential Relevance to Ribozyme Chemistry. *Chemistry & Biodiversity* **2008**, *5*, 1455–1474.

(27) Breaker, R. R.; Emilsson, G. M.; Lazarev, D.; Nakamura, S.; Puskarz, I. J.; Roth, A.; Sudarsan, N. A common speed limit for RNA-cleaving ribozymes and deoxyribozymes. *Rna* **2003**, *9*, 949–957.

(28) Frezza, E.; Courban, A.; Allouche, D.; Sargueil, B.; Pasquali, S. The interplay between molecular flexibility and RNA chemical probing reactivities analyzed at the nucleotide level via an extensive molecular dynamics study. *Methods* **2019**, *162*, 108–127.
(29) Hurst, T.; Xu, X.; Zhao, P.; Chen, S.-J. Quantitative understanding of SHAPE mechanism from RNA structure and dynamics analysis. *The Journal of Physical Chemistry B* **2018**, *122*, 4771–4783.

(30) Kirmizialtin, S.; Hennelly, S. P.; Schug, A.; Onuchic, J. N.; Sanbonmatsu, K. Y. *Methods in enzymology*; Elsevier, 2015; Vol. 553; pp 215–234.

(31) Calonaci, N.; Bernetti, M.; Jones, A.; Sattler, M.; Bussi, G. Molecular dynamics simulations with grand-canonical reweighting suggest cooperativity effects in RNA structure probing experiments. *J. Chem. Theory Comput.* **2023**, *19*, 3672–3685.

(32) Lu, C.; Ding, F.; Chowdhury, A.; Pradhan, V.; Tomsic, J.; Holmes, W. M.; Henkin, T. M.; Ke, A. SAM recognition and conformational switching mechanism in the Bacillus subtilis yitJ S box/SAM-I riboswitch. *Journal of molecular biology* **2010**, *404*, 803–818.

(33) Wang, J.; Cieplak, P.; Kollman, P. A. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? *J. Comput. Chem.* **2000**, *21*, 1049–1074.

(34) Petersson, G. A.; Bennett, A.; Tensfeldt, T. G.; Al-Laham, M. A.; Shirley, W. A.; Mantzaris, J. A complete basis set model chemistry. I. The total energies of closed-shell atoms and hydrides of the first-row elements. *J. Chem. Phys.* **1988**, *89*, 2193–2218.

(35) Petersson, G. A.; Al-Laham, M. A. A complete basis set model chemistry. II. Open-shell systems and the total energies of the first-row atoms. *J. Chem. Phys.* **1991**, *94*, 6081–6090.

(36) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. *J. Chem. Phys.* **1993**, *98*, 5648–5652.

(37) Frisch, M. J. et al. Gaussian 16 Revision C.01. 2016; Gaussian Inc. Wallingford CT.
(38) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. *J. Comput. Chem.* **2004**, *25*, 1157–1174.

(39) Lindahl, E.; Hess, B.; Van Der Spoel, D. GROMACS 3.0: a package for molecular simulation and trajectory analysis. *Molecular modeling annual 2001*, **7**, 306–317.

(40) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. *Bioinformatics* **2013**, *29*, 845–854, Place: England Publisher: Science for Life Laboratory, Stockholm and Uppsala, 171 21 Stockholm, Sweden.

(41) Hess, B.; Kutzner, C.; Van Der Spoel, D.; Lindahl, E. GROMACS 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation. *Journal of chemical theory and computation* **2008**, *4*, 435–447, Publisher: ACS Publications.

(42) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: fast, flexible, and free. *Journal of computational chemistry* **2005**, *26*, 1701–1718, Publisher: Wiley Online Library.

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

(44) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew. *The Journal of Chemical Physics* **2004**, *120*, 9665–9678.

(45) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N log (N) method for Ewald sums in large systems. *The Journal of chemical physics* **1993**, *98*, 10089–10092.
(46) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. *The Journal of chemical physics* **1995**, *103*, 8577–8593.

(47) Hess, B.; Bekker, H.; Berendsen, H. J.; Fraaije, J. G. LINCS: a linear constraint solver for molecular simulations. *Journal of computational chemistry* **1997**, *18*, 1463–1472.

(48) Harvey, S. C.; Tan, R. K.-Z.; Cheatham III, T. E. The flying ice cube: velocity rescaling in molecular dynamics leads to violation of energy equipartition. *Journal of computational chemistry* **1998**, *19*, 726–740.

(49) Berendsen, H. J.; Postma, J. P. M. v.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. *The Journal of chemical physics* **1984**, *81*, 3684–3690.

(50) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. *The Journal of chemical physics* **2007**, *126*, 14101, Publisher: AIP Publishing.

(51) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. *Journal of Applied physics* **1981**, *52*, 7182–7190.

(52) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A portable plugin for free-energy calculations with molecular dynamics. *Comput. Phys. Commun.* **2009**, *180*, 1961–1972.

(53) Grossfield, A. WHAM: the weighted histogram analysis method. Version 2.0.11, http://membrane.urmc.rochester.edu/wordpress/?page_id=126.

(54) Neese, F. Software update: the ORCA program system, version 5.0. *WIRES Comput. Molec. Sci.* **2022**, *12*, e1606.
(55) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions. *The Journal of Physical Chemistry B* **2009**, *113*, 6378–6396.

(56) Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method (SCC-DFTB). *Journal of Chemical Theory and Computation* **2011**, *7*, 931–948.

(57) Gaus, M.; Goez, A.; Elstner, M. Parametrization and Benchmark of DFTB3 for Organic Molecules. *Journal of Chemical Theory and Computation* **2012**, *9*, 338–354.

(58) Gaus, M.; Lu, X.; Elstner, M.; Cui, Q. Parameterization of DFTB3/3OB for Sulfur and Phosphorus for Chemical and Biological Applications. *Journal of Chemical Theory and Computation* **2014**, *10*, 1518–1537.

(59) Gaus, M.; Jin, H.; Demapan, D.; Christensen, A. S.; Goyal, P.; Elstner, M.; Cui, Q. DFTB3 Parametrization for Copper: The Importance of Orbital Angular Momentum Dependence of Hubbard Parameters. *Journal of Chemical Theory and Computation* **2015**, *11*, 4205–4219.

(60) Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. *J. Mol. Graph. 1996*, *14*, 33–8, 27–8.

(61) The MathWorks Inc., MATLAB.

(62) Tubiana, T.; Carvaillo, J.-C.; Boulard, Y.; Bressanelli, S. TTClust: a versatile molecular simulation trajectory clustering program with graphical summaries. *Journal of chemical information and modeling* **2018**, *58*, 2178–2182.

(63) McGibbon, R. T.; Beauchamp, K. A.; Harrigan, M. P.; Klein, C.; Swails, J. M.; Hernández, C. X.; Schwantes, C. R.; Wang, L.-P.; Lane, T. J.; Pande, V. S. MD-
Traj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories.
*Biophysical Journal* **2015**, *109*, 1528 – 1532.

(64) Harris, C. R. et al. Array programming with NumPy. *Nature* **2020**, *585*, 357–362.

(65) Taghavi, A.; Riveros, I.; Wales, D. J.; Yildirim, I. Evaluating geometric definitions of stacking for RNA dinucleoside monophosphates using molecular mechanics calculations.
*Journal of chemical theory and computation* **2022**, *18*, 3637–3653.

(66) Sedova, A.; Banavali, N. K. Geometric patterns for neighboring bases near the stacked state in nucleic acid strands. *Biochemistry* **2017**, *56*, 1426–1443.

(67) Sabei, A.; Hognon, C.; Martin, J.; Frezza, E. Dynamics of Protein–RNA Interfaces Using All-Atom Molecular Dynamics Simulations. *The Journal of Physical Chemistry B* **2024**, *128*, 4865–4886.

(68) Mezei, M. A new method for mapping macromolecular topography. *Journal of Molecular Graphics and Modelling* **2003**, *21*, 463–472.

(69) Ceres, N.; Pasi, M.; Lavery, R. A protein solvation model based on residue burial.
*Journal of Chemical Theory and Computation* **2012**, *8*, 2141–2144.

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

(71) 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.

(72) Arnold, E. B.; Cohn, D. D.; Bose, E. M.; Wolfe, G.; Jones, A. N. bioRxiv.

(73) Weeks, K. M.; Mauger, D. M. Exploring RNA structural codes with SHAPE chemistry.
*Accounts of chemical research* **2011**, *44*, 1280–1291.
(74) Karabiber, F.; McGinnis, J. L.; Favorov, O. V.; Weeks, K. M. QuShape: rapid, accurate, and best-practices quantification of nucleic acid probing information, resolved by capillary electrophoresis. *Rna* **2013**, *19*, 63–73.

(75) Hardouin, P.; Pan, N.; Lyonnet du Moutier, F.-X.; Chamond, N.; Ponty, Y.; Will, S.; Sargueil, B. IPANEMAP Suite: a pipeline for probing-informed RNA structure modeling. *NAR Genomics and Bioinformatics* **2025**, *7*, lqaf028.

(76) Hennelly, S. P.; Novikova, I. V.; Sanbonmatsu, K. Y. The expression platform and the aptamer: cooperativity between Mg2+ and ligand in the SAM-I riboswitch. *Nucleic Acids Research* **2012**, *41*, 1922–1935.

(77) Fernandez-de Cossio-Diaz, J.; Hardouin, P.; Lyonnet du Moutier, F.-X.; Di Gioacchino, A.; Marchand, B.; Ponty, Y.; Sargueil, B.; Monasson, R.; Cocco, S. Designing Molecular RNA Switches with Restricted Boltzmann Machines. *bioRxiv* **2024**,

(78) Leontis, N. B.; Westhof, E. Geometric nomenclature and classification of RNA base pairs. *Rna* **2001**, *7*, 499–512.
TOC Graphic

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