# Targeting RNA with small molecules using state-of-the-art methods provides highly predictive affinities of riboswitch inhibitors

Narjes Ansari, Chengwen Liu, Florent Hédin, Jérôme Hénin, Jay W Ponder, Pengyu Ren, Jean-Philip Piquemal, Louis Lagardère, Krystel El Hage

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

## To cite this version:

Narjes Ansari, Chengwen Liu, Florent Hédin, Jérôme Hénin, Jay W Ponder, et al. Targeting RNA with small molecules using state-of-the-art methods provides highly predictive affinities of riboswitch inhibitors. Communications Biology, 2025, 8 (1), pp.1405. 10.1038/s42003-025-08809-y. hal-05498925

HAL Id: hal-05498925

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

Submitted on 7 Feb 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Targeting RNA with small molecules using state-of-the-art methods provides highly predictive affinities of riboswitch inhibitors

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

Check for updates

Narjes Ansari Ⓔ<fnref n="1" />✉, Chengwen Liu<fnref n="1" />,², Florent Hédin Ⓔ<fnref n="1" />, Jérôme Hénin Ⓔ³, Jay W. Ponder Ⓔ<fnref n="1" />,⁴, Pengyu Ren<fnref n="1" />,², Jean-Philip Piquemal Ⓔ<fnref n="1" />,⁵, Louis Lagardère Ⓔ<fnref n="1" />,⁵✉ & Krystel El Hage Ⓔ<fnref n="1" />✉

Targeting RNA with small molecules represents a promising yet relatively unexplored avenue for the design of new drugs. Nevertheless, challenges arise from the lack of computational models and techniques able to accurately model RNA systems, and predict their binding affinities to small molecules. Here, we tackle these difficulties by developing a tailored state-of-the-art approach for absolute binding free energy calculations of RNA-binding small molecules. For this, we combine the advanced AMOEBA polarizable force field to the newly developed lambda-Adaptive Biasing Force scheme associated to refined restraints allowing for efficient sampling. To capture the free energy barrier associated to challenging RNA conformational changes, we combine machine learning-based collective variables with enhanced sampling simulations. Applying this computational protocol to a complex Riboswitch-like RNA target demonstrates quantitative predictions. These results pave the way for the routine application of free energy simulations in RNA-targeted drug discovery, thus providing a significant reduction in their failure rate.

RNA-targeting therapeutics offer an attractive alternative to reach traditionally undruggable proteins and expand the druggable target space in the human genome. Indeed, RNA can adopt three-dimensional (3D) structures, conferring varied functional roles in human biology as well as disease-related dysfunction. They are unique molecules capable of interacting with all three major forms of biological macromolecules: proteins, DNAs, and RNAs. Due to their diverse mechanisms of action, including gene silencing, splice modulation, and protein interaction, RNA therapeutics offer a versatile platform for drug development<fnref n="1" />⁻⁴.

Therapeutic strategies to target RNA include nucleotide-based agents such as (i) antisense oligonucleotides (ASOs), small interfering RNA (siRNA), and microRNAs (miRNAs) that can directly target messenger RNAs (mRNAs) and noncoding RNAs (ncRNAs) through Watson-Crick base pairing, (ii) clustered regularly inter-spaced short palindromic repeats (CRISPR) gene editing that can directly modify target RNA sequences to treat specific disorders, (iii) and small molecules that recognize RNA structures². While ASOs and CRISPR editing have been invaluable to the

field of chemical biology, their applications in therapeutics lacks efficient delivery strategies and presents significant adverse reactions⁵,⁶. These faced difficulties are due to their physico-chemical properties (size, charge,...), off-target effects, interactions with the immune system, rapid clearance and nuclease degradation. Alternatively, small molecules offer a promising alternative with potential for oral bio-availability and blood-brain barrier penetration. In particular, the extensive medicinal chemistry knowledge allows to systematically optimize pharmacokinetics and potency. Moreover, functional RNAs or RNA motifs are highly structured and form binding pockets or clefts that are accessible by small molecules.

Overall, targeting RNA with small molecules is emerging as a promising drug discovery approach. Such molecules interact with various target RNA substrates, including RNA structural motifs (such as hairpins, bulges, and internal loops), RNA enzymes (ribozymes), and specific RNA sequences⁷,⁸. Several binders have been identified and exhibit various modes of action (MOAs), from simple binding to direct cleavage or recruitment of endogenous nucleases. These molecules modulate diverse biological

<footnote n="1">¹Qubit Pharmaceuticals, 75014 Paris, France. ²Department of Biomedical Engineering, The University of Texas at Austin, Austin, TX, USA. ³Université Paris Cité, CNRS, Laboratoire de Biochimie Théorique, Paris, France. ⁴Department of Chemistry, Washington University in Saint Louis, Saint Louis, MO, USA. ⁵Sorbonne Université, Laboratoire de Chimie Théorique, UMR 7616 CNRS, Paris, France. ✉e-mail: narjesa@qubit-pharmaceuticals.com; louis.lagardere@sorbonne-universite.fr; krystel.elhage@qubit-pharmaceuticals.com</footnote>
processes, such as inhibiting bacterial and viral translation (e.g., ribocil and a
riboswitch in *Escherichia coli*, 2-aminobenzimidazole derivatives and the
hepatitis C internal ribosome entry site (IRES)¹⁰, respectively), and directing
splicing by acting as molecular glues with cellular proteins (e.g., branaplam¹¹
and the FDA-approved risdiplam¹²).

Challenges arise when tackling these systems throughout the drug discovery process and more specifically at the early stages, where the computer-aided drug design tools and simulation approaches are essential to characterize the target. Understanding its structural dynamics as well as its mechanism of action and interplay with other partners will guide the implementation of a pertinent inhibition or modulation mechanism. However, current methods face significant challenges in addressing the key obstacles associated with structure determination, accurate computational models, and reliable affinity prediction methods, which are the key elements required to guide the initial molecular design and rational drug optimization phases.

RNA's chemical composition is distinct from proteins due to a highly
electronegative surface potential associated to a limited buried surface area.
It is highly charged, dynamic, and can adopt different conformations. In
addition, divalent metal ions play an important role in its structural stability.
Since RNA is also surrounded by ions and polarizable water molecules, it
makes its targeting very challenging without having a reliable computational
model that encompasses an accurate treatment of electrostatics¹³,¹⁴. This can
be addressed by using second-generation force fields that account for many-
body polarization effects¹⁵,¹⁶. Indeed, their inclusion in simulation play a key
role in predicting biophysical properties governing the interactions between
a target, its environment, and a potential drug candidate¹⁷. This is the case of
AMOEBA, an advanced multipolar polarizable force field that employs
atomic induced dipoles to include polarization effects and atomic multi-
poles, up to quadrupoles, to represent the anisotropy of permanent
electrostatics¹⁸,¹⁹. These non-additive effects are essential for accurately
ranking affinities of small molecules¹³,²⁰,²¹. Although polarizable force fields
were traditionally considered as computationally slow, limiting their early
adoption, recent advances in computing hardware, especially GPUs ( Gra-
phics Processing Unit) have greatly improved the situation²². Further cou-
plings with enhanced sampling techniques strongly increased the molecular
dynamics (MD) timescales accessible to these simulations¹⁴,²³⁻²⁸.

Computing free energies of binding in silico has become a routinely
used tool in drug discovery projects²⁹,³⁰. Among the existing techniques,
alchemical ones, where the transformation of interest is made through a
non-physical path, has proven to be both efficient and reliable, especially for
relative free energies of binding (RBFE) between compounds³¹. Computing
Absolute Binding Free Energies (ABFE) is more challenging but also more
directly useful for drug discovery as it does not rely on the knowledge of a
reference compound and is not limited to ligands sharing a similar binding
mode. Computing binding affinities of small molecules bound to RNA or
DNA with such methods remain a daunting challenge when compared to
proteins³²⁻³⁶. Given the particular nature of RNA, it is considerably more
sensitive to the protocol adopted than in the case of small molecule-protein
complexes, especially when the ligand is large, charged, with highly flexible,
elongated structure. In addition, both the RNA target and the ligand remain
exposed to the solvent and the counterions, in contrast to ligand-protein
complexes in which the ligand is in a partly buried protein cavity and
interacts with only a limited number of "discrete" water molecules without
intrusions from a counterion, which is common in RNA systems.

Furthermore, the shape of the ligand and its binding mode, i.e.,
intercalation or groove binding, will impact the overall protocol used to get
ABFEs, especially in terms of the choice and the number of restraints used to
keep the ligand in place during simulations and facilitate convergence³⁷,³⁸.
Other aspects, such as the discretization of the alchemical path (i.e., the so-
called "lambda windows" used in standard methods) and the sampling time
are critical factors that should be taken into account when developing a
tailored protocol within a well-defined workflow for ABFE calculations for
RNA (and DNA) small molecule binders²⁷,³²,³⁹. Moreover, the free energy
barrier associated with the significant conformational change that may exist

between the Apo and Holo structures of the target is notoriously challenging
to estimate¹⁰.

A recently developed technique, the lambda-Adaptive Biasing Force (lambda-ABF) approach, grounded on λ-dynamics combined with a newly developed multiple-walker adaptive biasing force bias, bypasses the discretization of the alchemical path and enables efficient sampling of orthogonal degrees of freedom through the variable nature of the alchemical variable. In combination with a recently introduced restraint scheme³⁸, distance-to-bound-configuration (DBC), it has shown improved performance compared to more standard methods²⁷.

In this work, the inhibition of the hepatitis C internal ribosome entry site (HCV-IRES) IIa subdomain by nineteen 2-aminobenzimidazole derivatives is taken as a case study to address the aforementioned challenges and to illustrate how advanced and state-of-the-art computational techniques can be employed to compute the affinity of a class of small molecules for their RNA target with high accuracy.

The HCV-IRES is essential for the initiation of viral protein synthesis and adopts well-defined folds that are potential targets for antiviral translation inhibitors. The 3D structure of the IRES domain IIa in complex with a benzimidazole translation inhibitor at 2.2 Å resolution (PDB ID 3TZR) was determined and provides insights into inhibiting its function¹⁰,⁴¹. This system is chosen for its complexity that makes it one of the most challenging systems for atomistic simulations and absolute binding affinity calculations. First, the ligand-binding site has an intricate architecture, organized by a metal spine at the back of the cavity. Additionally, it contains three magnesium ions as intrinsic structural components (Fig. 1). Moreover, the inhibitors are 2-aminobenzimidazole derivatives containing 2 to 3 positive charges, as well as several aromatic and non-aromatic cycles with lengthy arms (see Fig. 2)⁴². Finally, the comparison of the structures of bound and unbound RNA shows that the RNA undergoes a dramatic ligand-induced conformational adaptation from an L-shaped to an extended form, creating a deep pocket that resembles the substrate binding site in riboswitches. Also, the aforementioned three magnesium ions undergo adaptive reorganization upon binding of the benzimidazole ligand⁴¹.

In order to comprehend the binding mode of these inhibitors and reproduce their binding affinities, we have resorted to a computational model that accounts for the underlying effects of polarization and electron anisotropy, AMOEBA<sup>18,43,44</sup>, recently parameterized for DNA and RNA using quantum chemistry (QC) methods<sup>19</sup> and tested in several oligonucleotide simulations<sup>32,33,45</sup>. We have also used the massively parallel GPU-accelerated molecular dynamics Tinker-HP software package<sup>22,23,46,47</sup> coupled to the Colvars library<sup>27,48</sup>. We have put in place a novel state-of-the-art approach for ABFE calculations tailored for riboswitch-like RNA using the lambda-ABF scheme<sup>27</sup> combined with positional, orientational, and conformational restraints<sup>38</sup>. We also propose a strategy to estimate the free energy difference associated to the conformational change between the Apo and the Holo structure of the riboswitch through advanced simulation techniques leveraging machine-learned collective variables (CVs)<sup>49</sup>.

Results

In this section, we first assess the binding mode of the studied small
molecules to their RNA target, followed by computing their binding affi-
nities. Finally, we explore the conformational changes between the Apo and
Holo states of the RNA.

Binding mode analysis

The RNA architecture of the ligand-binding pocket consists of the phos-
phate of U56, which is rotated into the RNA helix major groove, and two
magnesium ions anchored between the phosphate oxygen atoms and the
bases of A57 and U59 (Fig. 1b). The extreme contortion of the RNA
backbone that directs the U56 phosphate into the helix interior is further
pinned in place by hydrogen bonds between the C55 phosphate and 2'-
hydroxyl groups of the flanking residues A54 and U56. The elaborate net-
work of backbone and metal interactions forces the base of A54 to project
away from the RNA whereas residues C55 and U56 tightly pack in the RNA
**Fig. 1** | Structural features of the riboswitch-like RNA-Benzimidazole complex. a Overall view of the riboswitch-like RNA complex with Benzimidazole (ligand 4) after equilibration. The ligand is depicted in fluorescent blue. Three Mg²⁺ ions, shown in green, are coordinated by surrounding water molecules. b Schematic representation of interactions within the ligand-binding site. Hydrogen bonds are depicted by dashed lines. Non-Watson-Crick base pairs are shown with solid lines and symbols according to Leontis et al.⁸⁷. Stacked lines (≡) denote base stacking and intercalation of the ligand, which is colored pink. Residues that interact with the benzimidazole are in red. c Detailed view of the ligand-binding site. The bases of G52 and A53, which form the intercalation site for the benzimidazole scaffold, are shown in pink. The purine base of G110 is colored yellow. Hydrogen bond interactions between G110 and the ligand are depicted with dashed lines and labeled as HB₁ and HB₂. Possible hydrogen bonds between the NH⁺ group of the ligand and A109 and G52 are also shown by dashed lines. d Surface representation highlighting the ligand-binding pocket. The interaction regions of A53 and G52 are colored pink, G110 is shown in yellow, and A109 is depicted in orange.

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

helix major groove as well as against the backbone of C58 and U59 (Fig. 1d). The third magnesium ion closes part of the solvent exposed part of the ligand site by bridging the RNA strands with interactions at the base of A53 and the C108 phosphate. The top and bottom part of the binding pocket are stabilized by base triples which form through cross-bracing interactions along the RNA helix with A53 docking at the Hoogsteen edge of A109 and A57 interacting with the sugar edge of C111. Additional hydrogen bonds stabilize the distorted internal loop in the benzimidazole complex and include interactions of the U56 phosphate with the C58 base and of the C55 2'-hydroxyl with the A53 base (Fig. 1b).

Understanding this intricate architecture is key to analyzing how the ligands interact with the RNA structure. It is well-known that the contribution of hydrogen bonding and stacking interactions are predominant in RNA recognition⁵⁰. Following the detailed relaxation process described in the simulation setup, we analyzed the binding modes obtained from the 40 ns MD simulations. Figure 3 illustrates the binding modes of six representative ligands, out of nineteen, within the target RNA. These ligands, ligand 1 to 6, cover different ranges of affinities and include a non-binder, ligand 1 (see Fig. 2). Three primary interactions are identified as crucial for stabilizing the ligand within the binding pocket: (i) a π-π stacking interaction occurs between the ligand and the nucleic bases G52 and A53, (ii) two hydrogen bonds (labeled HB₁ and HB₂) formed between the central core of the ligand (benzimidazole) and G110 (see Fig. 1c), and (iii) an additional ionic interaction is established between the NH⁺ moiety of the dimethylazanium group at the end of each arm of the ligand and the non-bridging phosphate oxygen atoms (PO) of A109 (or G110) and the PO of G52 (or A53) on the opposite strand (depicted as N<sub>a</sub> and N<sub>b</sub>, respectively in Fig. 3).

With these key interactions established, their evolution as a function of time was computed throughout the trajectory and for each ligand; the normalized distributions are shown in Fig. 4. For the π-π stacking interactions, we computed the distance between the center of mass (labeled as COM) of the purine rings of G52 and A53, and the COM of the benzimidazole scaffold of the ligand (referred to as LIG<sub>COM</sub> in Fig. 4). For the distances between G52<sub>COM</sub> and LIG<sub>COM</sub>, the results reveal that for most ligands, the distribution is narrow and centers around 4.1 to 4.4 Å, consistent with the typical π-π interaction range of 3.3 to 4.4 Å. In contrast, for ligands 1 and 2, the distribution is broader, with ligand 1 centering around 4.8 Å while showing a shoulder at 4.0 Å. As for ligand 2, the distribution centers around 4.3 Å with a shoulder at 5 Å. This is concomitant with the fact that ligand 1 is a non-binder, while ligand 2 is the weakest binder in the list. Similarly, the distribution of distances between A53<sub>COM</sub> and LIG<sub>COM</sub> shows a comparable trend for all ligands. Here, the distribution also centers around a narrow range, specifically between 3.6 and 3.8 Å for ligands 3 to 6. This indicates that the overall trend in interaction distances is consistent with the π-π stacking interactions observed with G52<sub>COM</sub>, albeit with a slightly shorter average distance for the interactions involving A53<sub>COM</sub>.

Next, the hydrogen bond interactions, HB₁ and HB₂, were analyzed. For these interactions, ligand 1 exhibits a broader distribution, suggesting a weaker and more transient interaction. In comparison, ligand 6 exhibits a narrower distribution, centered around 2 Å with the highest intensities, reflecting strong and consistent interactions. Ligands 3, 4, and 5 show a peak similar to ligand 6 for HB₂, although with slightly lower intensity for HB₁. Ligand 2 shows a similar peak as the other ligands; however, it has a tail at low probability, suggesting occasional weaker interactions.
<figure><img src="image_6.png" /><figcaption>**Fig. 2** | 2D chemical structures of the nineteen studied translation inhibitors. 2D structure of the scaffold with common atoms labeled 1 to 7. The potential growing regions are highlighted as R₁ and R₂. For each ligand, the dissociation constant (K_D) and the corresponding size of RNA are reported from Seth et al.⁴².</figcaption></figure>

Finally, regarding the ionic interaction, the N<sub>a</sub> atom of ligands 2, 3, and 4 forms a hydrogen bond with the phosphate (PO) atoms of A109, as indicated by a sharp peak around 2 Å. Ligand 6 is unique in forming an HB with the PO atom of G110. For the second arm of the ligand (N<sub>b</sub> atom), ligands 4, 5, and 6 alternate between the PO atom of G52 and that of A53, as evidenced by multiple peaks in the distribution, one of which is centered around 2 Å. Ligands 1, 2, and 3 are also the only ligands that do not form a HB with the PO atom of G52, nor with A53, except for a rare HB with the PO atom of G52 at 2 Å for ligand 3.

## Binding affinity calculations

To accurately predict the binding affinities of extended charged ligands (bearing 2 to 3 charges) bound to a long 29-mer riboswitch-like RNA with three structural Mg<sup>2+</sup> ions as intrinsic components in the vicinity of the binding pocket, we have developed a specific protocol utilizing the lambda-ABF method and DBC restraints.

As free energy is a state function, the computed ΔG<sub>o</sub> are fundamentally independent of the intermediate states sampled during an alchemical simulation. This allows for tailored modifications in these states to ensure better sampling and convergence as it is routinely done in restraining the ligand to keep it in the binding site while scaling down its non-bonded interactions with the environment. As mentioned earlier, different types of restraints exist that can either be positional or positional and orientational. However, these two methods suffer from distinct limitations<sup>27,38</sup>, especially when considering elongated ligands containing two flexible arms. What is

required in such systems is positional, orientational as well as conformational restraints, that can be defined through the newly developed DBC collective variable, which is the root mean square deviation (RMSD) of the ligand in the moving reference frame of the binding site<sup>38</sup>. Therefore, to facilitate sampling and maintain the ligand in the binding site, we employed DBC restraints between the benzimidazole ring of all ligands and heavy atoms from the RNA binding site. The selected atoms are highlighted in Supplementary Fig. S1.

Besides, it is well known that starting alchemical simulations from Holo structures may lead to incomplete relaxation towards the Apo structure<sup>51</sup>. This incomplete relaxation can result in a significant free energy difference, potentially causing errors of several kcal/mol in the final prediction. Typically, to address this, all predicted affinities are shifted by a common value to account for this discrepancy<sup>31</sup>. This adjustment assumes that all systems experience the same degree of relaxation from the Holo to the Apo state. To ensure the validity of this assumption, facilitate the sampling and maintain structural stability of the RNA-ligand complex, we implemented slight positional restraints with a force constant of 1 kcal/mol/A<sup>2</sup>) on the phosphate (P) atoms of RNA and Mg<sup>2+</sup> ions. These adjustments were designed to avoid impacting the fully bound state, ensuring the integrity of the complex throughout the simulations. This approach allows us to independently estimate and account for the relaxation effect, thus improving the accuracy of our absolute binding affinity calculations.

In the case of our system, it is known that benzimidazole inhibitors can induce a conformational change characterized by a widened interhelical
<figure><img src="image_7.png" /><figcaption>**Fig. 3** | 3D representation of the binding modes of six representative ligands and their interactions with the environment. a–f show the binding modes of Ligands 1 to 6, respectively. Ligands are depicted in licorice style and colored in cyan, while the phosphate backbone is represented by a yellow meshed ribbon. Green arrows</figcaption></figure>

indicate hydrogen bond donors, and purple circles indicate π-π stacking interactions between the G52, A53 bases and the ligand. The positively charged nitrogen (N) atoms are shown in red with a superscript (+). The chirality labels are displayed as S or R for chiral centers.

angle in the IRES subdomain Ila<sup>10,41,52</sup>. This conformation facilitates the undocking of subdomain IIB from the ribosome, leading to the inhibition of IRES-driven translation in HCV-infected cells. As illustrated in Fig. 5, the L-shaped RNA structure in the Apo state is known to adopt a straighter configuration<sup>10,41,52</sup>. Apart from the restructuring of several nucleic bases, this conformational change also includes a rearrangement of Mg<sup>2+</sup> ion positions (Fig. 5a, b). To mitigate any potential bias from such known conformational changes, we applied the restraints described earlier. This approach helps ensure that such known rearrangements are properly accounted for, thereby improving the reliability of our ABFE calculations.

To obtain the alchemical free energy surface for each system at T = 300 K, we ran three replicas, each utilizing four walkers. Calculations were conducted in both the complex (bound state) and solvent (ligand in the bulk) phases. For the complex phase, simulations ran for ~50 ns for both van der Waals (vdW) and electrostatic (ELE) legs. In contrast, a 20 ns run was enough for convergence in the solvent phase for both VDW and ELE components. Raw data from each replica of the six representative ligands with varying ranges of affinities (Ligand 1 to 6) are provided in the Supplementary Information, Section II and Supplementary Tables S2 and S3. Convergence plots are also provided as Supplementary Figs. S2 to S25. The free energy cost associated with the release of the DBC restraint is readily computed in the gas phase using TI<sup>28,53</sup>. This is achieved by gradually releasing the restraint to a compatible harmonic distance restraint, the cost of which can be computed analytically (see Supplementary Table S1).

As described in the Methods section, we have employed two distinct setups for ABFE calculations, one consisting of a neutralized system (Neut.) for all ligands and a second including a physiological ion concentration (Phys.) for the six representative ligands (Ligand 1 to 6). The computed ΔG values of all ligands in neutralized conditions are reported in Table 1 and compared with the experimental data. Supplementary Table S4, compares the computed ΔG values in neutralized condition (ΔG<sub>Comp-Neut</sub>) and in the presence of physiological ion concentration (ΔG<sub>Comp-Phys</sub>) for the six representative ligands. The raw absolute binding free energies (ΔG<sub>Comp-raw</sub>) consistently exhibit stronger values (favorable) compared to the experimental results. In order to account for the Apo-Holo

conformational change, the raw binding free energy (ΔG<sub>Comp-raw</sub>) is converted to the final binding free energy ΔG<sub>Comp</sub> by adding a constant shift, which is calculated as follows ΔG<sub>shift</sub> = (Σi=1 to N Gi<sup>Exp</sup> - Σi=1 to N Gi<sup>Comp</sup>)<sup>51</sup>, so that the mean of computed free energies is equal to the mean of experimental free energies. This shift yields a value of 10.31 ± 1.1 kcal/mol for the neutral setup taking into account all ligands. Taking into account the six representative ligands, a value of 10.35 ± 0.91 kcal/mol is obtained for the neutral setup and 11.33 ± 1.19 kcal/mol in the presence of physiological ion concentration. Despite this slight difference, the results are consistent within the error margins.

The correlation plots are illustrated in Fig. 6 and Supplementary Fig. S26. Upon including the systematic shift between the calculated and experimental binding free energies to accommodate the Apo-Holo protein reorganization free energy, we observe a correlation between computed and experimental results across the entire dataset of 0.6 (Pearson r), with a root mean square error (RMSE) of 1.03 kcal/mol and a mean absolute error (MAE) of 0.85 kcal/mol for the neutralized setup (Fig. 6). It's worth noting that ligands 1, 9, 11, 12, 13, 14, 17, 18, 19 were excluded from the analysis as they are predicted to be non-binders in agreement with the experimental measurement which showed no binding at 100 or 200 µM. As for ligand 11, also a non-binder, it leaves the pocket after 70 ns of plain MD, and hence was not considered for ABFE calculations. These results show a solid reproduction of the binding affinities and a good ranking, with ligand 6 and ligand 4 identified as the most potent and ligands 1, 9, 11, 12, 13, 14, 17, 18, 19 as non-binders. Taking into account only the six representative ligands, we observe a similar correlation between the two set-ups: a correlation of 0.7, with an RMSE of 0.81 kcal/mol and an MAE of 0.70 kcal/mol for the neutralized setup (Supplementary Fig. S26a); a correlation of 0.7 with an RMSE of 1.05 kcal/mol and an MAE of 0.84 kcal/mol for the set-up with a physiological concentration of salt (Supplementary Fig. S26b). Overall, both setups show similar satisfactory results.

Even though the experimental K<sub>D</sub> measurements were done with 2 different constructs, this shouldn't have a large effect on the affinities, as seen for the example of ligand 3, where K<sub>D</sub>=10 µ M in 40-mer compared to 17 µ M in the 29-mer<sup>42</sup>, which translate into 0.3 kcal/mol of difference, or on
<figure><img src="image_8.png" /><figcaption>**Fig. 4** | Normalized distributions of key interaction quantities for six representative ligands. The distributions for distances between a G52<sub>COM</sub> and LIG<sub>COM</sub>, b A53<sub>COM</sub> and LIG<sub>COM</sub>, c HB<sub>1</sub> distances, d HB<sub>2</sub> distances, distances between</figcaption></figure>

e A109<sub>PO</sub> and LIG<sub>Na</sub>, f G110<sub>PO</sub> and LIG<sub>Na</sub>, g LIG<sub>Nb</sub> and G52<sub>PO</sub>, h LIG<sub>Nb</sub> and A53<sub>PO</sub>. The colors represent different ligands as indicated in the legend.
<figure><img src="image_9.png" /><figcaption>**Fig. 5 | Conformational differences between Apo and Holo forms of RNA.**</figcaption></figure>

a L-shaped structure of the Apo state (PDB ID: 2nok) and b straight structure of the Holo state (PDB ID: 3tzr). c Alignment of the Apo and Holo structures. The

positions of Mg²⁺ ions in each structure are illustrated, denoted as *i*<sub>a</sub> and *i*<sub>h</sub> for the Apo and Holo states respectively, where *i* = 1 – 3. The key bases within the binding pocket are highlighted, with G110 in dark yellow and G52 and A53 in pink.

ranking the selected compounds in our study, it nonetheless might slightly effect the observed correlations. This is also without considering the error on the experimentally measured affinities.

**Performance.** To illustrate the computational cost required for one ABFE calculation, we have conducted benchmarking tests of lambda-ABF simulations on three different GPU architectures: NVIDIA H100, A100 and V100. This was done on the complex phase which is the most computationally expensive one. The performance on different GPU types is reported in Supplementary Table S5. As expected, it is lower on A100 and V100 GPUs due to differences in computational power and memory bandwidth. The performance on H100 GPUs is 1.5 times faster than A100 and 2.4 times faster than V100 GPUs. Details are provided in Section IV of Supplementary Information. All in all, the total amount of simulation time for the ABFEs ran in the present study encompasses more than 28 µs.

**Apo-Holo conformational change**

As mentioned in the above subsection, the $\Delta G_{\text{shift}}$ resulting from the energy difference between the computed $\Delta G$ and experimental values is estimated to be between 10.31 and $11.35 \pm 1.0$ kcal/mol, primarily reflects the Apo-Holo conformational change. To substantiate these findings, employing more direct computational methods is crucial. Given the substantial free energy barrier associated with the Apo-Holo transition, capturing this difference accurately requires advanced sampling techniques. As illustrated in Fig. 5, the Apo-Holo conformational change involves complex collective motions that are challenging to capture with standard methods. Simple CVs, such as distances or dihedrals, are insufficient for this system. Instead, more advanced machine learning-based CVs can be employed. To address this, we employed the OPES-Explore⁵⁴,⁵⁵ enhanced sampling method in conjunction with Deep-Linear Discriminant Analysis (Deep-LDA)⁴⁹, a machine learning-based technique specifically designed to identify effective CVs for complex transitions.

To identify the Deep-LDA CV, we conducted an unbiased simulation of 40 ns from both Apo and Holo structures, incorporating 47 descriptors that included distances and angles between various P atoms and three Mg²⁺ ions. The setup and training of the neural network is detailed in Supplementary Information Section V. Subsequently, we performed exploratory OPES simulations, each lasting 10 to 15 ns, with applied biases ranging from 2 to 15 kcal/mol on the Deep-LDA CV. Notably, a transition from the Apo to the Holo state was observed around 10 ns when the applied bias exceeded 13 kcal/mol (see Supplementary Figs. S27 and S28), which aligns well with the

energy change estimated from our lambda-ABF calculations. However, the conformational changes of specific bases is not captured by the P and Mg atoms in our Deep-LDA CV as well as the still approximate nature of the OPES strategy.

Despite the inherent challenges in achieving precise free energy calculations given the system's complexity, our approach provides a valuable qualitative approximation of the free energy barrier associated with the Apo-Holo conformational change.

**Discussion**

The complexity and dynamic nature of RNA, combined with its critical role in various biological processes, makes it a very difficult target in drug discovery. However, with these challenges, it becomes essential to harness the latest computational techniques in order to predict the binding affinities between drug candidates and the target of interest with high precision.

In this study, we present a state-of-the-art approach for calculating ABFEs of RNA-binding small molecules for application in drug discovery, showcasing its effectiveness in a challenging RNA system: the HCV-IRES riboswitch-like RNA. This system presents a highly charged and flexible RNA structure, an intricate binding site with three magnesium ions as intrinsic components, ligands with multiple positive charges and several aromatic and non-aromatic cycles with lengthy arms, as well as a substantial ligand-induced conformational change between the Apo and Holo states.

All in all, by combining the computational performance of the GPU-accelerated Tinker-HP package to enhanced sampling and machine learning, while leveraging the accuracy of the AMOEBA force field, we are able to model challenging drug discovery targets such as RNA and to compute absolute binding affinities to guide drug design with the lambda-ABF methodology. This synergy not only overcomes the limitations of available free energy methods but also establishes a new standard in computational modeling of RNA-ligand interactions. The present work represents a new achievement in computational methodologies for drug design, where the complexities of RNA and other challenging targets can be addressed with unparalleled accuracy and speed, thus guiding the development of RNA-targeted therapeutics with precision and confidence.

**Methods**

**System preparation**

The initial input for the simulation comprises the 3D X-ray structure of riboswitch-like RNA in the presence of a benzimidazole derivative (ligand 4 shown in Figs. 1 and 2d), obtained from the Protein Data Bank (PDB ID: 3TZR, resolution 2.21 Å). Using this structure as a starting point, we
<table><caption>**Table 1 | Experimental and computed ΔG values for each ligand**</caption><thead><tr><th>Ligand</th><th>K<sub>D</sub>(μM)</th><th>ΔG<sub>Exp</sub><br>(kcal/mol)</th><th>ΔG<sub>Comp</sub>(ΔG<sub>Comp-raw</sub>)<br>(kcal/mol)</th></tr></thead><tbody><tr><td>1</td><td>&gt; 200</td><td>&gt; -5.04</td><td>-0.12 (-10.43) ± 0.92</td></tr><tr><td>2</td><td>40</td><td>-6.01</td><td>-6.65 (-16.96) ± 0.68</td></tr><tr><td>3</td><td>8</td><td>-6.96</td><td>-7.30 (-16.88) ± 1.54</td></tr><tr><td>4</td><td>0.86</td><td>-8.27</td><td>-7.17 (-17.48) ± 0.65</td></tr><tr><td>5</td><td>3.50</td><td>-7.44</td><td>-7.25 (-17.56) ± 0.81</td></tr><tr><td>6</td><td>0.72</td><td>-8.39</td><td>-9.64 (-19.95) ± 1.54</td></tr><tr><td>7</td><td>1.70</td><td>-7.88</td><td>-7.26 (-17.57) ± 1.78</td></tr><tr><td>8</td><td>100</td><td>-5.45</td><td>-4.49 (-14.79) ± 0.56</td></tr><tr><td>9</td><td>&gt; 200</td><td>&gt; -5.04</td><td>-3.56 (-13.87) ± 0.53</td></tr><tr><td>10</td><td>125</td><td>-5.32</td><td>-7.30 (-17.61) ± 0.53</td></tr><tr><td>11</td><td>&gt; 200</td><td>&gt; -5.04</td><td>—</td></tr><tr><td>12</td><td>&gt; 200</td><td>&gt; -5.04</td><td>-0.71 (-11.01) ± 1.28</td></tr><tr><td>13</td><td>&gt; 200</td><td>&gt; -5.04</td><td>1.07 (-9.24) ± 1.38</td></tr><tr><td>14</td><td>&gt; 200</td><td>&gt; -5.04</td><td>2.92 (-7.39) ± 1.66</td></tr><tr><td>15</td><td>50</td><td>-5.86</td><td>-6.73 (-17.04) ± 1.03</td></tr><tr><td>16</td><td>10</td><td>-6.82</td><td>-5.32 (-15.63) ± 1.50</td></tr><tr><td>17</td><td>&gt; 100</td><td>&gt; -5.04</td><td>-2.81 (-13.12) ± 1.11</td></tr><tr><td>18</td><td>&gt; 100</td><td>&gt; -5.04</td><td>1.40 (-8.90) ± 0.92</td></tr><tr><td>19</td><td>&gt; 100</td><td>&gt; -5.04</td><td>3.04 (-7.27) ± 0.83</td></tr></tbody></table>

Computed ΔG<sub>Comp</sub> values are reported under neutralized conditions. The computed ΔG<sub>Comp</sub> and errors represent the mean and standard error of the mean from three replicas for each ligand (n = 3). Raw computed ΔG values, ΔG<sub>Comp-raw</sub>, are shown in parentheses. For Ligands 4, 5, 6 and 7 the experimental ΔG<sub>Exp</sub> is measured on a 40-mer RNA, while for the rest, the 29-mer is used. All computed ΔG were measured using a 29-mer construct. ΔG values are reported in kcal/mol. K<sub>D</sub> values are reported in μM.

generated the input for all other ligands illustrated in Fig. 2. Concerning ligands 5 and 6 which are diastereomers they were designed to be conformationally restricted benzimidazoles to increase the binding, and the desired conformation was separated by HPLC<sup>42</sup>. Hence, for our simulations the desired 3D conformation was generated based on ligand 4 (SSO in 3TZR). For consistency, we used a 29-mer RNA, as depicted in Fig. 1b, with all compounds. In the work of Seth et al.<sup>42</sup>, ligands 4, 5, 6 and 7 were measured with a 40-mer construct, while ligands 1, 2, and 3 were measured with the 29-mer construct. However, this does not significantly affect the K<sub>D</sub> value or the ranking, given that the results for ligand 3 reported in Tables 2 and 3 of Seth et al.<sup>42</sup> are consistent. Ligand 3 was tested with both the 29-mer and 40-mer RNA and it was found to have a slightly stronger binding to the 29-mer (K<sub>D</sub> = 10μM in the 40-mer compared to 17 μM in the 29-mer). Given the experimental measurement error, we expect the results to be within the same range.

### AMOEBA force field parameters

The systems studied here are highly charged and dynamic, while being surrounded by ions and polarizable water molecules. To accurately model these systems, we employed the AMOEBA force field, which accounts for polarization and anisotropy through a multipolar representation of electrostatics. AMOEBA has been recently parameterized for DNA and RNA based on QC methods and validated in several oligonucleotide simulations<sup>19</sup>. Additionally, we used the Tinker-HP software on GPUs<sup>46</sup> for efficient computational performance.

### Generating AMOEBA force field parameters for small molecules.

Ligand parameterization is performed using QC calculations and the Poltype package<sup>56</sup> which facilitate the automated generation of AMOEBA parameters. The process begins with an SDF file defining the atom

coordinate, bond order, ionic state, etc. Prior to geometry optimization, one conformation was selected from 500 conformers generated by the RDKit tool. Structure characteristics including intra-molecular hydrogen bonds, radius of gyration, solvent-accessible surface area were considered during the selection. Key torsion angles were kept frozen to the values of the selected conformer during subsequent MP2/6-31G(d) geometry optimization to maintain the preferred conformation. Two single-point (SP) calculations were performed on the optimized structure using the Psi4 software<sup>57</sup>: MP2/6-311G(d,p) for low-level and MP2/aug-cc-pvttz for high-level. Electron density from the low-level SP was used with the GDMA program<sup>58</sup> (with switch value 0) to generate initial atomic multipoles (charges, dipoles, and quadrupoles). Electron density from high-level SP was employed to calculate the electrostatic potential (ESP) on defined grids surrounding the molecule. The POTENTIAL program within Tinker 8 software<sup>59</sup> further refined the dipoles and quadrupoles by fitting them to reproduce the ESP. Valence (bond, angle, stretch-bend, and opbend), vdW, and torsion terms were matched using SMARTS patterns and assigned from published parameter databases<sup>60</sup>. Notably, vdW parameters for aromatic rings were assigned using those from a previous study<sup>61</sup> due to their improved interactions with DNA and RNA. For any missing torsion parameters, a conventional spin-and-fit procedure was used. This involved restrained optimization using the xtb program<sup>62</sup> with GFN-2 model<sup>63</sup> and SP energy calculations at ωb97x-d<sup>64</sup>/6-311+G(d) level in Psi4 software<sup>57</sup>. To minimize computational cost and potential steric clashes during rotating the bond of the whole ligand, a fragment molecule containing the fitted torsions was extracted from the ligand and used for torsion parameter fitting.

### MD simulations

All systems were neutralized by adding the appropriate number of counterions (K<sup>+</sup>) in a water box of 80 Å × 80 Å × 80 Å such that the distance between the FLAP and its periodic images is at least 16 Å. The systems contain ~52,000 atoms and were subject to NVT and NPT simulations for equilibration before production.

The systems underwent initial minimization using Tinker-HP “minimize” program, followed by a gradual heating process in 250 ps intervals. The temperature was increased from 50 K to 100 K, 200 K, and 300 K under the NVT ensemble, with RNA-ligand atoms restrained using a force constant of 10 kcal/mol/A<sup>2</sup>). During this step, MD simulations were conducted using the Verlet integrator with a time step of 1.0 fs. Subsequently, the systems were equilibrated in the NPT ensemble at 1 atm using the same restraints with force constants of 10 and 5 kcal/mol/A<sup>2</sup>), each for 0.5 ns, using the RESPA integrator<sup>65</sup> and a time step of 2 fs. To relax the water molecules around the ions at 300 K, a 10 ns restrained MD simulation was run for each system with a force constant of 1 kcal/mol/A<sup>2</sup>), followed by simulations where the restraint was applied only to phosphate (P) atoms and Mg<sup>2+</sup> ions with a force constant of 1 kcal/mol/A<sup>2</sup>), run for 5 ns. The ligands were stabilized by employing distance-to-bound-configuration (DBC) restraints<sup>38</sup> through the Colvars library<sup>48</sup>, constraining their RMSD from the bound configuration. For DBC definition, all the ligands' heavy atoms and the RNA heavy atoms within a 4 Å radius of the ligand were included. Temperature and pressure were controlled by the Bussi<sup>66</sup> thermostat and Berendsen<sup>67</sup> barostat, respectively. Van der Waals interactions utilized a 9 Å cut-off, while electrostatic interactions were treated by particle mesh ewald (PME)<sup>68</sup> with a real-space cutoff of 7 Å. Induced dipoles were calculated with a Preconditioned Conjugate Gradient solver with a convergence tolerance of 1 × 10<sup>-5</sup> Debye<sup>69</sup>.

The polyanionic nature of RNA attracts water molecules and counterions, influencing its stability. Additionally, Mg<sup>2+</sup> ions can shed their hydration shells and establish strong interactions with specific sites in the RNA, promoting local structural rearrangements. To mitigate artifacts during the production phase, positional restraints were applied to the phosphate (P) atoms and Mg<sup>2+</sup> ions with a force constant of 1 kcal/mol/A<sup>2</sup>). To prevent fraying, distance restraints were added to enforce the hydrogen bonds between the terminal base pairs. The partially restrained simulations
<caption>**Fig. 6 | Experimental vs computed ΔG values.**</caption>

Computed ΔG and errors are the mean and standard error of the mean from three replicas for each ligand in the neutralized complex (n = 3). The dark shaded region spans ±1 kcal/mol; the lighter region spans ±2 kcal/mol. The color bar represents the absolute value difference between experimental and computed values. Pearson r, RMSE, and MAE are reported in the right bottom corner of each plot. The non-binders are not shown in the plots.

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

were run for 40 ns for production; however, for all non-binders, simulations were run for 100 ns.

## DBC coordinate for binding pose determination and restraint definition

Accurate alchemical absolute free energy calculations require the use of restraints to ensure proper convergence. Typically, these restraints involve harmonic restraints between centers of mass of the ligand and the binding site (positional restraints)⁷⁰ or more sophisticated methods like 'Boresch' restraints⁷¹, which act on distances, angles, and dihedrals (positional and orientational restraints). In this study, we employed a more refined approach using DBC collective variables. DBC measures the RMSD of ligand coordinates relative to the moving frame of the receptor's binding site³⁸. This metric captures positional, orientational and conformational deviation of the ligand in a single collective variable. The restraint is flat-bottomed with a starting value that matches the distribution characteristic of the binding mode. This distribution can be efficiently monitored using the Colvars Dashboard⁷² within the VMD software³⁵. Depending on the binding scenario, this distribution will be narrow for poses that are tightly bound, while it will be broader for systems with a diverse set of loosely bound configurations. The primary goal of these restraints is to restrict sampling during the decoupling process to those configurations that are relevant in the fully coupled state and thus facilitate convergence.

## lambda-ABF for ABFE simulations

As stated above, we have resorted here to alchemical free energy simulations in which the ligand is progressively “alchemically” decoupled from its environment, first in complex with the riboswitch host, and separately from the bulk solvent. The standard free energy of binding is obtained through a thermodynamic cycle. Such techniques require the sampling of the alchemical Hamiltonians. State-of-the-art methods typically use a pre-defined λ schedule with simulations at fixed λ values and a post-processing

stage in which free energy differences are computed using estimators such as Thermodynamic Integration (TI)⁷⁴ or Free Energy Perturbation⁷⁵. However, these fixed-λ simulations can hinder the relaxation of orthogonal degrees of freedom that may benefit from a variable λ exploration which can be made possible through additional enhanced sampling techniques such as Hamiltonian Replica Exchange⁷⁶ or Expanded Ensemble methods⁷⁷. Here, we used the newly developed lambda-ABF approach²⁷ which enables efficient relaxation of orthogonal degrees of freedom due to the variable nature of λ, as well as ergodic sampling thanks to the ABF bias. Its sampling efficiency has been shown to be superior compared to more traditional techniques. We utilized its implementation within Tinker-HP, in combination with the Colvars library¹⁸, which facilitates user-friendly convergence estimation without requiring post-processing, and allows integration with other CV-based methods.

For calculating the free energy surface of our systems at T = 300 K, we used four walkers with the lambda-ABF method and decoupled vdW and ELE interactions sequentially in separate legs. Each walker was run for 50 ns per leg. We ran three independent replicas for each setup. As for the free energy cost associated with the release of the DBC restraint, it is computed in the gas phase through TI by progressively releasing it to a compatible harmonic distance restraint which can then be computed analytically²⁷,³⁸,⁵³.

During the alchemical decoupling, the phosphate atoms of the RNA backbone and the Mg²⁺ ions were restrained with a slight force constant of 1 kcal/mol/A²). It allows to, (i) prevent any destructuration, (ii) measure the ligand binding contribution independently from any Apo-Holo conformation changes by fixing the endpoint of the thermodynamic cycle.

Concerning the ion concentration used, it is well known that it can significantly impact the results of binding affinities⁷⁸, especially for charged systems such as RNA-ligand complexes⁷⁹. Reproducing all experimental conditions is crucial to get meaningful results in free energy simulations. However, methodological issues can arise when computing electrostatics using standard methods like PME in periodic boundary conditions.
Although finite size effects diminish with larger simulation boxes, they may still be non-negligible in practice⁸⁰. A notable issue occurs with PME in non-neutralized boxes which happen by default when scaling down electrostatics of ligands during alchemical simulations. In the reference experimental setup for mass spectrometry (MS) based binding assay that was used to measure the K<sub>D</sub>, the incorporation of a physiological ion concentration was not evoked.⁴² Given the positively charged nature of all the ligands reported here (2 to 3 positive charges) we cannot scale down the charge of anions to maintain a neutral charge of the simulation box in this setup. We have thus employed two distinct protocols for ABFE calculations:

* Neutralized setup (Neut.): K⁺ counter-ions were added to neutralize the system, which led to non neutral boxes in the full decoupling state (total charge of 2 to 3, depending on the ligand). This setup allows us to assess the impact of non neutral boxes and the absence of physiological ions on binding affinities.

* Physiological Ion Concentration (Phys.): 0.15 Molar of KCl were added to the system 2 to 3 Cl⁻ ions had their charge scaled down during ligand scaling (neutral boxes maintained). This setup aims at replicating physiological conditions more closely and evaluating the effect of physiological ionic strength on the computed binding affinities.

## OPES explore

OPES Explore⁵⁵ is a CV-based enhanced sampling technique which is an evolution of metadynamics⁸¹. It aims at broadening the system's sampling to a target probability distribution known as the Well-tempered distribution, defined as

$$p^{\text{tg}}(\mathbf{s}) \propto [P(s)]^\frac{1}{\gamma}$$

where P(s) represents the unbiased marginal distribution of chosen CVs, and γ is the bias factor controlling the broadening.

To achieve this, Gaussian kernels are employed to reconstruct p<sup>tg</sup>(s), which in turn determines the bias potential V(s) through the following recursive strategy at step n:

$$V(s) = (\gamma - 1)k_B T \log \left( \frac{p^{\text{tg}}(s)}{Z} + \epsilon \right) \\ p^{\text{tg}}(s) = \frac{1}{n} \sum_{k=1}^{n} G_k(s, s_k)$$

where k<sub>B</sub> is the Boltzmann constant, T is the temperature set by the thermostat, Z is a normalization factor, and G<sub>k</sub>(s, s<sub>k</sub>) represents the Gaussian kernel deposited at step k.

The initial Gaussian kernel width, often denoted as SIGMA, typically corresponds to the standard deviation of s in the initial basin. The bias factor γ is typically set using γ = ΔE/(k<sub>B</sub>T), where ΔE is the barrier parameter. The regularization term ε ensures stability and is related to ΔE through ε = e⁻<sup>-ΔE/(k<sub>B</sub>T(1-1/γ))</sup>.

The barrier parameter ΔE sets a limit on the maximum bias energy that OPES⁵⁴ can introduce into the system. It should be large enough to facilitate transitions away from the initial basin while preventing the system from accessing irrelevant high-energy states that are difficult to reverse.

For further details, refer to the original literature introducing the OPES method⁵⁴,⁵⁵.

## Enhanced CV design using Deep-LDA

Deep-LDA⁴⁹ is a way to define CV that correspond to a transition between two conformational states. It is rooted in LDA⁸², a classical technique in classification tasks. Here, it is tailored to differentiate between Apo and Holo states observed in unbiased system simulations.

Standard LDA optimizes Fisher's ratio w<sup>T</sup>S<sub>w</sub>w, where w defines the discriminant direction s(R) = w<sup>T</sup>d(R). Here, S<sub>w</sub> and S<sub>w</sub> represent between- and within-class scatter matrices respectively, calculated from descriptor distributions in Apo and Holo states.

Deep-LDA extends LDA by employing a neural network (NN) to process N<sub>d</sub> descriptors d, optimizing LDA within the last hidden layer h. The resulting Deep-LDA CV s = w<sup>T</sup>h offers a non-linear, expressive projection that has proven effective in various applications⁸³⁻⁸⁵. To enhance usability in simulations, we apply a cubic transformation, s<sub>w</sub> = s + s<sup>383,86</sup>, ensuring smoother distributions.

## Statistics and reproducibility

All the free energy simulations performed in this study have been produced three times by varying the random seed that generates the initial velocities. The final results are the average of these three repeats and the standard deviation is calculated among these replicates. Each leg of the complex phase (riboswitch in complex with the ligands of interest) has been run for 50 ns and each leg of the solvent phase (where the ligand lies alone in bulk water) has been run for 20 ns.

## Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

## Data availability

The input files used in this study are available on GitHub: https://github.com/krystel-elhage/targeting_rna. Numerical source data for graphs are provided as Supplementary Data 1; it contains the source data behind the graphs in Figures 4 and 6, and Supplementary Figs. S26, S27 and S28.

## Code availability

The code used during the study is available on GitHub: https://github.com/TinkerTools/tinker-hp, https://github.com/TinkerTools/poltype2, https://github.com/Colvars/colvars.

Received: 8 November 2024; Accepted: 28 August 2025;

Published online: 01 October 2025

## References

1. Harries, L. W. RNA biology provides new therapeutic targets for human disease. *Front. Genet.* **10**, 205 (2019).

2. Zhu, Y., Zhu, L., Wang, X. & Jin, H. RNA-based therapeutics: an overview and prospectus. *Cell Death Dis.* **13**, 644 (2022).

3. Yu, A.-M., Choi, Y. H. & Tu, M.-J. RNA drugs and RNA targets for small molecules: principles, progress, and challenges. *Pharmacol. Rev.* **72**, 862–898 (2020).

4. Lightfoot, H. L. & Smith, G. F. Targeting RNA with small molecules-a safety perspective. *Br. J. Pharmacol.* **182**, 4201–4220 (2023).

5. Bennett, C. F. Therapeutic antisense oligonucleotides are coming of age. *Annu. Rev. Med.* **70**, 307–321 (2019).

6. Luther, D., Lee, Y., Nagaraj, H., Scaletti, F. & Rotello, V. Delivery approaches for crispr/cas9 therapeutics in vivo: advances and challenges. *Expert Opin. Drug Deliv.* **15**, 905–913 (2018).

7. Wu, P. Inhibition of RNA-binding proteins with small molecules. *Nat. Rev. Chem.* **4**, 441–458 (2020).

8. Childs-Disney, J. L. et al. Targeting RNA structures with small molecules. *Nat. Rev. Drug Discov.* **21**, 736–762 (2022).

9. Howe, J. A. et al. Selective small-molecule inhibition of an RNA structural element. *Nature* **526**, 672–677 (2015).

10. Dibrov, S. M. et al. Hepatitis C virus translation inhibitors targeting the internal ribosomal entry site: miniperspective. *J. Med. Chem.* **57**, 1694–1707 (2014).

11. Palacino, J. et al. Smn2 splice modulators enhance u1-pre-mRNA association and rescue sma mice. *Nat. Chem. Biol.* **11**, 511–517 (2015).

12. Ratni, H. et al. Discovery of risdiplam, a selective survival of motor neuron-2 (smn2) gene splicing modifier for the treatment of spinal muscular atrophy (SMA). *J Med Chem.* **61**, 6501–6517 (2018).
13. Gresh, N. et al. Addressing the issues of non-isotropy and non-additivity in the development of quantum chemistry-grounded polarizable molecular mechanics. in *Quantum Modeling of Complex Molecular Systems*, 1-49 (Springer, 2015).

14. Jing, Z. et al. Polarizable force fields for biomolecular simulations: recent advances and applications. *Annu. Rev. Biophys.* **48**, 371–394 (2019).

15. Shi, Y., Ren, P., Schnieders, M. & Piquemal, J.-P. Polarizable force fields for biomolecular modeling. *Rev. Comput. Chem.* **28**, 51–86 (2015).

16. Melcr, J. & Piquemal, J.-P. Accurate biomolecular simulations account for electronic polarization. *Front. Mol. Biosci.* **6**, 143 (2019).

17. El Khoury, L. et al. Computationally driven discovery of SARS-CoV-2 M pro inhibitors: from design to experimental validation. *Chem. Sci.* **13**, 3674–3687 (2022).

18. Ponder, J. W. et al. Current status of the amoeba polarizable force field. *J. Phys. Chem. B* **114**, 2549–2564 (2010).

19. Zhang, C. et al. Amoeba polarizable atomic multipole force field for nucleic acids. *J. Chem. Theory Comput.* **14**, 2084–2108 (2018).

20. Gresh, N., Cisneros, G. A., Darden, T. A. & Piquemal, J.-P. Anisotropic, polarizable molecular mechanics studies of inter-and intramolecular interactions and ligand- macromolecule complexes. a bottom-up strategy. *J. Chem. Theory Comput.* **3**, 1960–1986 (2007).

21. El Hage, K., Piquemal, J.-P., Hobaika, Z., Maroun, R. G. & Gresh, N. Substituent-modulated affinities of halobenzene derivatives to the HIV-1 integrase recognition site. Analyses of the interaction energies by parallel quantum chemical and polarizable molecular mechanics. *J. Phys. Chem. A* **118**, 9772–9782 (2014).

22. Adjoua, O. et al. Tinker-hp: Accelerating molecular dynamics simulations of large complex systems with advanced point dipole polarizable force fields using GPUs and multi-GPU systems. *J. Chem. Theory Comput.* **17**, 2034–2053 (2021).

23. Lagardère, L., Aviat, F. & Piquemal, J.-P. Pushing the limits of multiple-time-step strategies for polarizable point dipole molecular dynamics. *J. Phys. Chem. Lett.* **10**, 2593–2599 (2019).

24. Jaffrelot-Inizan, T. et al. High-resolution mining of SARS-CoV-2 main protease conformational space: Supercomputer-driven unsupervised adaptive sampling. *Chem. Sci.* **12**, 4889–4907 (2021).

25. Célerse, F. et al. An efficient Gaussian-accelerated molecular dynamics (GAMD) multilevel enhanced sampling strategy: application to polarizable force fields simulations of large biological systems. *J. Chem. Theory Comput.* **18**, 968–977 (2022).

26. Célerse, F., Lagardère, L., Derat, E. & Piquemal, J.-P. Massively parallel implementation of steered molecular dynamics in tinker-hp: Comparisons of polarizable and non-polarizable simulations of realistic systems. *J. Chem. Theory Comput.* **15**, 3694–3709 (2019).

27. Lagardère, L. et al. Lambda-abf: Simplified, portable, accurate, and cost-effective alchemical free-energy computation. *J. Chem. Theory Comput.* **20**, 4481–4498 (2024).

28. Blazhynska, M. et al. Water-glycan interactions drive the SARS-CoV-2 spike dynamics: insights into glycan-gate control and camouflage mechanisms. *Chem. Sci.* **15**, 14177–14187 (2024).

29. Wang, L. et al. Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field. *J. Am. Chem. Soc.* **137**, 2695–2703 (2015).

30. Zhang, C.-H. et al. Potent noncovalent inhibitors of the main protease of sars-cov-2 from molecular sculpting of the drug perampanel guided by free energy perturbation calculations. *ACS Cent. Sci.* **7**, 467–475 (2021).

31. Chodera, J. D. et al. Alchemical free energy methods for drug discovery: progress and challenges. *Curr. Opin. Struct. Biol.* **21**, 150–160 (2011).

32. El Hage, K. et al. Targeting the major groove of the palindromic d (ggcgcc) 2 sequence by oligopeptide derivatives of anthraquinone intercalators. *J. Chem. Inf. Model.* **62**, 6649–6666 (2022).

33. Gresh, N. et al. Enforcing local DNA kinks by sequence-selective trisintercalating oligopeptides of a tricationic porphyrin: a polarizable molecular dynamics study. *ChemPhysChem* **25**, e202300776 (2024).

34. El Hage, K., Mondal, P. & Meuwly, M. Free energy simulations for protein ligand binding and stability. *Mol. Sim.* **44**, 1044–1061 (2018).

35. Rasouli, A., Pickard IV, F. C., Sur, S., Grossfield, A. & Işık Bennett, M. Essential considerations for free energy calculations of RNA-small molecule complexes: lessons from the theophylline-binding RNA aptamer. *J. Chem. Inf. Model.* **65**, 223–239 (2025).

36. Abramyan, A. M. et al. Accurate physics-based prediction of binding affinities of RNA- and DNA-targeting ligands. *J. Chem. Inf. Model.* **65**, 1392–1403 (2025).

37. Clark, F., Robb, G., Cole, D. J. & Michel, J. Comparison of receptor-ligand restraint schemes for alchemical absolute binding free energy calculations. *J. Chem. Theory Comput.* **19**, 3686–3704 (2023).

38. Salari, R., Joseph, T., Lohia, R., Hénin, J. & Brannigan, G. A streamlined, general approach for computing ligand binding free energies and its application to GPCR-bound cholesterol. *J. Chem. Theory Comput.* **14**, 6560–6573 (2018).

39. Zhang, H. et al. Accurate estimation of the standard binding free energy of netropsin with DNA. *Molecules* **23**, 228 (2018).

40. Gapsys, V. et al. Accurate absolute free energies for ligand-protein binding based on non-equilibrium approaches. *Commun. Chem.* **4**, 61 (2021).

41. Dibrov, S. M. et al. Structure of a hepatitis C virus RNA domain in complex with a translation inhibitor reveals a binding mode reminiscent of riboswitches. *Proc. Natl Acad. Sci.* **109**, 5223–5228 (2012).

42. Seth, P. P. et al. Sar by ms: discovery of a new class of RNA-binding small molecules for the hepatitis C virus: internal ribosome entry site IIA subdomain. *J. Med. Chem.* **48**, 7099–7102 (2005).

43. Wu, J. C., Chattree, G. & Ren, P. Automation of amoeba polarizable force field parameterization for small molecules. *Theor. Chem. Acc.* **131**, 1–11 (2012).

44. Shi, Y. et al. Polarizable atomic multipole-based amoeba force field for proteins. *J. Chem. Theory Comput.* **9**, 4046–4063 (2013).

45. Yang, X., Liu, C., Kuo, Y.-A., Yeh, H.-C. & Ren, P. Computational study on the binding of mango-II RNA aptamer and fluorogen using the polarizable force field amoeba. *Front. Mol. Biosci.* **9**, 946708 (2022).

46. Lagardère, L. et al. Tinker-hp: a massively parallel molecular dynamics package for multiscale simulations of large complex systems with advanced point dipole polarizable force fields. *Chem. Sci.* **9**, 956–972 (2018).

47. Jolly, L.-H. et al. Raising the performance of the tinker-hp molecular modeling package [article v1.0]. *Living J. Comput. Mol. Sci.* **1**, 10409 (2019).

48. Fiorin, G., Klein, M. L. & Hénin, J. Using collective variables to drive molecular dynamics simulations. *Mol. Phys.* **111**, 3345–3362 (2013).

49. Bonati, L., Rizzi, V. & Parrinello, M. Data-driven collective variables for enhanced sampling. *J. Phys. Chem. Lett.* **11**, 2998–3004 (2020).

50. Padroni, G., Patwardhan, N., Schapira, M. & Hargrove, A. Systematic analysis of the interactions driving small molecule–rna recognition. *RSC Med. Chem.* **11**, 802–813 (2020).

51. Chen, W. et al. Enhancing hit discovery in virtual screening through absolute protein-ligand binding free-energy calculations. *J. Chem. Inf. Model.* **63**, 3171–3185 (2023).

52. Parsons, J. et al. Conformational inhibition of the hepatitis C virus internal ribosome entry site RNA. *Nat. Chem. Biol.* **5**, 823–825 (2009).

53. Santiago-McRae, E., Ebrahimi, M., Sandberg, J. W., Brannigan, G. & Hénin, J. Computing absolute binding affinities by streamlined alchemical free energy perturbation (safep)[article v1. 0]. *Living J. Comput. Mol. Sci.* **5**, 2067–2067 (2023).

54. Invernizzi, M. & Parrinello, M. Rethinking metadynamics: from bias potentials to probability distributions. *J. Phys. Chem. Lett.* **11**, 2731–2736 (2020).
55. Invernizzi, M. & Parrinello, M. Exploration vs convergence speed in adaptive-bias enhanced sampling. *J. Chem. Theory Comput.* **18**, 3988-3996 (2022).

56. Walker, B., Liu, C., Wait, E. & Ren, P. Automation of amoeba polarizable force field for small molecules: Poltype 2. *J. Comput. Chem.* **43**, 1530-1542 (2022).

57. Turney, J. M. et al. Psi4: an open-source ab initio electronic structure program. *Wiley Interdiscip. Rev. Comput. Mol. Sci.* **2**, 556-565 (2012).

58. Stone, A. J. & Alderton, M. Distributed multipole analysis: methods and applications. *Mol. Phys.* **56**, 1047-1064 (1985).

59. Rackers, J. A. et al. Tinker 8: software tools for molecular design. *J. Chem. theory Comput.* **14**, 5273-5289 (2018).

60. Ren, P., Wu, C. & Ponder, J. W. Polarizable atomic multipole-based molecular mechanics for organic molecules. *J. Chem. Theory Comput.* **7**, 3143-3161 (2011).

61. Zhang, C., Bell, D., Harger, M. & Ren, P. Polarizable multipole-based force field for aromatic molecules and nucleobases. *J. Chem. Theory Comput.* **13**, 666-678 (2017).

62. Bannwarth, C. et al. Extended tight-binding quantum chemistry methods. *Wiley Interdiscip. Rev. Comput. Mol. Sci.* **11**, e1493 (2021).

63. Bannwarth, C., Ehlert, S. & Grimme, S. Gfn2-xtb-an accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions. *J. Chem. Theory Comput.* **15**, 1652-1671 (2019).

64. Chai, J.-D. & Head-Gordon, M. Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections. *Phys. Chem. Chem. Phys.* **10**, 6615-6620 (2008).

65. Tuckerman, M., Berne, B. J. & Martyna, G. J. Reversible multiple time scale molecular dynamics. *J. Chem. Phys.* **97**, 1990-2001 (1992).

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

67. Berendsen, H. J., Postma, J. V., Van Gunsteren, W. F., DiNola, A. & Haak, J. R. Molecular dynamics with coupling to an external bath. *J. Chem. Phys.* **81**, 3684-3690 (1984).

68. Essmann, U. et al. A smooth particle mesh ewald method. *J. Chem. Phys.* **103**, 8577-8593 (1995).

69. Lagardère, L. et al. Scalable evaluation of polarization energy and associated forces in polarizable molecular dynamics: li. toward massively parallel computations using smooth particle mesh ewald. *J. Chem. Theory Comput.* **11**, 2589-2599 (2015).

70. Laury, M. L., Wang, Z., Gordon, A. S. & Ponder, J. W. Absolute binding free energies for the sample cucurbit [8] uril host-guest challenge via the amoeba polarizable force field. *J. Comput. Aided Mol. Des.* **32**, 1087-1095 (2018).

71. Boresch, S., Tettinger, F., Leitgeb, M. & Karplus, M. Absolute binding free energies: a quantitative approach for their calculation. *J. Phys. Chem. B* **107**, 9535-9551 (2003).

72. Hénin, J., Lopes, L. J. & Fiorin, G. Human learning for molecular simulations: the collective variables dashboard in VMD. *J. Chem. Theory Comput.* **18**, 1945-1956 (2022).

73. Humphrey, W., Dalke, A. & Schulten, K. Vmd: visual molecular dynamics. *J. Mol. Graph.* **14**, 33-38 (1996).

74. Straatsma, T. P. & McCammon, J. A. Multiconfiguration thermodynamic integration. *J. Chem. Phys.* **95**, 1175-1188 (1991).

75. Zwanzig, R. W. High-temperature equation of state by a perturbation method. i. nonpolar gases. *J. Chem. Phys.* **22**, 1420-1426 (1954).

76. Jiang, W. & Roux, B. Free energy perturbation hamiltonian replica-exchange molecular dynamics (fep/h-remd) for absolute ligand binding free energy calculations. *J. Chem. Theory Comput.* **6**, 2559-2565 (2010).

77. Lyubartsev, A., Martsinovski, A., Shevkunov, S. & Vorontsov-Velyaminov, P. New approach to Monte Carlo calculation of the free

energy: method of expanded ensembles. *J. Chem. Phys.* **96**, 1776-1783 (1992).

78. Thomas, J. R. & Hergenrother, P. J. Targeting RNA with small molecules. *Chem. Rev.* **108**, 1171-1224 (2008).

79. Far, S. et al. Bis-and tris-DNA intercalating porphyrins designed to target the major groove: Synthesis of acridylbis-arginyl-porphyrins, molecular modelling of their DNA complexes, and experimental tests. *Eur. J. Org. Chem.* **2004**, 1781-1797 (2004).

80. Petrov, D., Perthold, J. W., Oostenbrink, C., de Groot, B. L. & Gapsys, V. Guidelines for free-energy calculations involving charge changes. *J. Chem. Theory Comput.* **20**, 914-925 (2024).

81. Barducci, A., Bonomi, M. & Parrinello, M. Metadynamics. *WIREs Comput. Mol. Sci.* **1**, 826-843 (2011).

82. Welling, M. Fisher Linear Discriminant Analysis. Tech. Rep., Dep. Comput. Sci. Univ. Toronto (2005).

83. Rizzi, V., Bonati, L., Ansari, N. & Parrinello, M. The role of water in host-guest interaction. *Nat. Commun.* **12**, 93 (2021).

84. Ansari, N., Rizzi, V., Carloni, P. & Parrinello, M. Water-triggered, irreversible conformational change of SARS-CoV-2 main protease on passing from the solid state to aqueous solution. *J. Am. Chem. Soc.* **143**, 12930-12934 (2021).

85. Ansari, N., Rizzi, V. & Parrinello, M. Water regulates the residence time of benzamidine in trypsin. *Nat. Commun.* **13**, 5438 (2022).

86. Bjelobrk, Z. et al. Naphthalene crystal shape prediction from molecular dynamics simulations. *Cryst. Eng. Comm.* **21**, 3280-3288 (2019).

87. Leontis, N. B., Stombaugh, J. & Westhof, E. The non-Watson-Crick base pairs and their associated isostericity matrices. *Nucleic Acids Res.* **30**, 3497-3531 (2002).

## Acknowledgements

We thank the Grand Equipement de Calcul Intensif (GENCI), Institut du Développement et des Ressources en Informatique(IDRIS), Centre Informatique de l'Enseignement Supérieur (CINES), France, for supporting this work under the grant no A0140714153 and A0160714153 (K.E.H.).

## Author contributions

N.A., C.L., L.L., K.E.H. performed research; N.A., C.L., F.H., J.W.P., P.R., J.H., J.-P.P., L.L., K.E.H. investigation; N.A., L.L., K.E.H. analyzed the results and wrote the paper with input from all authors; L.L. developed new code and theoretical tools; K.E.H. conceived and supervised the project, and acquired funding.

## Competing interests

J.W.P., P.R., J.-P.P., and L.L. are co-founders and shareholders of Qubit Pharmaceuticals. All other authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s42003-025-08809-y.

**Correspondence** and requests for materials should be addressed to Narjes Ansari, Louis Lagardère or Krystel El Hage.

**Peer review information** *Communications Biology* thanks Anna Bochicchio and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor: Dr Ophelia Bu.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2025