F P Panei

P Gkeka

M Bonomi
email: mbonomi@pasteur.fr

Identifying small-molecules binding sites in RNA conformational ensembles with SHAMAN

The rational targeting of RNA with small molecules is hampered by our still limited understanding of RNA structural and dynamic properties. Most in silico tools for binding site identification rely on static structures and therefore cannot face the challenges posed by the dynamic nature of RNA molecules.

Here we present SHAMAN, a computational technique to identify potential small-molecule binding sites in RNA structural ensembles. SHAMAN enables exploring the conformational landscape of RNA with atomistic molecular dynamics and at the same time identifying RNA pockets in an efficient way with the aid of probes and enhanced-sampling techniques. In our benchmark composed of large, structured riboswitches as well as small, flexible viral RNAs, SHAMAN successfully identified all the experimentally resolved pockets and ranked them among the most favorite probe hotspots. Overall, SHAMAN sets a solid foundation for future drug design efforts targeting RNA with small molecules, effectively addressing the long-standing challenges in the field.

.

Introduction

RNA molecules, initially thought to be only carriers of genetic information from gene to proteins, are now known to perform a variety of biological functions, such as regulating the process of protein synthesis and defending against the entry of foreign nucleic acids into cells [START_REF] Cech | The noncoding RNA revolution -Trashing old rules to forge new ones[END_REF][START_REF] Cable | Noncoding RNAs: biology and applications-a Keystone Symposia report[END_REF][START_REF] Mortimer | Insights into RNA structure and function from genome-wide studies[END_REF][START_REF] Yao | Cellular functions of long noncoding RNAs[END_REF] . Alongside these findings, modulation of RNA functions is becoming a promising therapeutic approach for treating diseases such as cancer, viral infections, cardiovascular and muscular disorders, and neurodegenerative conditions [START_REF] Wang | RNA therapeutics on the rise[END_REF][START_REF] Damase | The Limitless Future of RNA Therapeutics[END_REF][START_REF] Halloy | Innovative developments and emerging technologies in RNA therapeutics[END_REF] .

Beside classical approaches, such as the design of antisense oligonucleotides interfering with mRNAs or directly editing RNA with CRISPR-Cas9, targeting RNA with small molecules is emerging as a promising strategy [START_REF] Rizvi | RNA as a small molecule druggable target[END_REF][START_REF] Falese | Targeting RNA with small molecules: from fundamental principles towards the clinic[END_REF][START_REF] Disney | Targeting RNA with Small Molecules To Capture Opportunities at the Intersection of Chemistry, Biology, and Medicine[END_REF][START_REF] Warner | Principles for targeting RNA with drug-like small molecules[END_REF] in terms of number of potential targets, bioavailability, and delivery [START_REF] Warner | Principles for targeting RNA with drug-like small molecules[END_REF][START_REF] Kole | RNA therapeutics: beyond RNA interference and antisense oligonucleotides[END_REF][START_REF] Kaczmarek | Advances in the delivery of RNA therapeutics: from concept to clinical reality[END_REF][START_REF] Winkle | Noncoding RNA therapeuticschallenges and potential solutions[END_REF][START_REF] Luther | Delivery approaches for CRISPR/Cas9 therapeutics in vivo : advances and challenges[END_REF] .

Although in recent years the research in this field has surged [START_REF] Morgan | Survey and summary insights into the development of chemical probes for RNA[END_REF][START_REF] Meyer | Small molecule recognition of disease-relevant RNA structures[END_REF] , the number of FDA-approved drugs is still limited and the compounds currently available on the market were identified exclusively by costly and time-consuming experimental screenings [START_REF] Howe | Selective small-molecule inhibition of an RNA structural element[END_REF][START_REF] Ratni | Discovery of Risdiplam, a Selective Survival of Motor Neuron-2 (SMN2) Gene Splicing Modifier for the Treatment of Spinal Muscular Atrophy (SMA)[END_REF][START_REF] Hashemian | Linezolid: a review of its properties, function, and use in critical care[END_REF] .

Computer-aided drug design (CADD) provides several essential tools to assist various stages of drug discovery, from druggability assessment to virtual screening for hit identification, binding affinity calculations, and generative methods for lead optimization. While these tools are well established for proteins, their application to RNA molecules is still in its infancy. The available biochemical and structural data is gradually elucidating the chemical properties of RNA binders [START_REF] Yazdani | Machine Learning Informs RNA-Binding Chemical Space[END_REF] and the structural properties of RNA binding sites [START_REF] Panei | HARIBOSS: a curated database of RNA-small molecules structures to aid rational drug design[END_REF] . This knowledge has been stimulating the development of ligand- [START_REF] Mehta | SMMRNA: a database of small molecule modulators of RNA[END_REF][START_REF] Kumar Mishra | NALDB: nucleic acid ligand database for small molecules targeting nucleic acid[END_REF] and 2D structure- [START_REF] Sun | RNALigands: a database and web server for RNA-ligand interactions[END_REF][START_REF] Donlic | R-BIND 2.0: An Updated Database of Bioactive RNA-Targeting Small Molecules and Associated RNA Secondary Structures[END_REF][START_REF] Disney | Inforna 2.0: A Platform for the Sequence-Based Design of Small Molecules Targeting Structured RNAs[END_REF] based virtual screening approaches, 3D binding-site detection tools [START_REF] Rekand | DrugPred_RNA-A Tool for Structure-Based Druggability Predictions for RNA Binding Sites[END_REF][START_REF] Zeng | Rsite2: an efficient computational method to predict the functional sites of noncoding RNAs[END_REF][START_REF] Wang | RLBind: a deep learning method to predict RNA-ligand binding sites[END_REF][START_REF] Kognole | SILCS-RNA: Toward a Structure-Based Drug Design Approach for Targeting RNAs with Small Molecules[END_REF][START_REF] Su | Recognition of small molecule-RNA binding sites using RNA sequence and structure[END_REF] , docking software [START_REF] Ruiz-Carmona | rDock: A Fast, Versatile and Open Source Program for Docking Ligands to Proteins and Nucleic Acids[END_REF][START_REF] Feng | NLDock: a Fast Nucleic Acid-Ligand Docking Algorithm for Modeling RNA/DNA-Ligand Complexes[END_REF][START_REF] Jiang | RLDOCK method for predicting RNA-small molecule binding modes[END_REF][START_REF] Guilbert | Docking to RNA via Root-Mean-Square-Deviation-Driven Energy Minimization with Flexible Ligands and Flexible Targets[END_REF] and scoring functions [START_REF] Stefaniak | AnnapuRNA: A scoring function for predicting RNA-small molecule binding poses[END_REF][START_REF] Chhabra | RNAPosers: Machine Learning Classifiers for Ribonucleic Acid-Ligand Poses[END_REF][START_REF] Pfeffer | DrugScore RNA Knowledge-Based Scoring Function To Predict RNA-Ligand Interactions[END_REF][START_REF] Philips | LigandRNA: computational predictor of RNA-ligand interactions[END_REF] specific for RNA molecules. However, our understanding of the structural and dynamic properties of RNA molecules and their interaction with small molecules still remains limited, thus ultimately hindering the rational design of novel and effective compounds [START_REF] Manigrasso | Computer-aided design of RNA-targeted small molecules: A growing need in drug discovery[END_REF] .

In the cellular context, function-specific biological signals trigger complex multi-step RNA conformational changes that in turn guide a variety of RNA functions, such as ligand sensing and signaling, catalysis, or co-transcriptional folding [START_REF] Ganser | The roles of structural dynamics in the cellular functions of RNAs[END_REF][START_REF] Ken | RNA conformational propensities determine cellular activity[END_REF] . These conformational changes and the underlying dynamics are influenced both by the inherent flexibility of RNA molecules, i.e. many large-scale motional modes spanning a variety of timescales, and other cellular co-factors [START_REF] Al-Hashimi | RNA dynamics: it is about time[END_REF] . Despite the significant efforts to characterize RNA dynamics using both experimental [START_REF] Soni | Structural basis for specific RNA recognition by the alternative splicing factor RBM5[END_REF] , in-silico [START_REF] Šponer | RNA Structural Dynamics As Captured by Molecular Simulations: A Comprehensive Overview[END_REF] , and integrative approaches [START_REF] Bernetti | Integrating experimental data with molecular simulations to investigate RNA structural dynamics[END_REF] , most available tools for CADD, and in particular for the identification of small molecules binding sites, still rely on a static description of RNA structure [START_REF] Rekand | DrugPred_RNA-A Tool for Structure-Based Druggability Predictions for RNA Binding Sites[END_REF][START_REF] Zeng | Rsite2: an efficient computational method to predict the functional sites of noncoding RNAs[END_REF][START_REF] Wang | RLBind: a deep learning method to predict RNA-ligand binding sites[END_REF][START_REF] Kognole | SILCS-RNA: Toward a Structure-Based Drug Design Approach for Targeting RNAs with Small Molecules[END_REF][START_REF] Su | Recognition of small molecule-RNA binding sites using RNA sequence and structure[END_REF] . The only exception is SILCS-RNA [START_REF] Kognole | SILCS-RNA: Toward a Structure-Based Drug Design Approach for Targeting RNAs with Small Molecules[END_REF] where potential binding sites are identified by exploring the conformation of the target RNA with small cosolvent probes, similar to mixed-solvent approaches already extensively used for proteins [START_REF] Defelipe | Solvents to Fragments to Drugs: MD Applications in Drug Design[END_REF] . While SILCS-RNA can describe small structural rearrangements induced by the probes, it is not designed to capture large RNA conformational changes and, therefore, it is not able to detect binding sites present in metastable states that are marginally populated yet crucial for therapeutic applications [START_REF] Manigrasso | Computer-aided design of RNA-targeted small molecules: A growing need in drug discovery[END_REF][START_REF] Ganser | The roles of structural dynamics in the cellular functions of RNAs[END_REF][START_REF] Ken | RNA conformational propensities determine cellular activity[END_REF][START_REF] Salmon | A general method for constructing atomic-resolution RNA ensembles using NMR residual dipolar couplings: The basis for interhelical motions revealed[END_REF] .

Here, we present SHAdow Mixed solvent metAdyNamics (SHAMAN), a computational technique for binding site identification in dynamic RNA structural ensembles. Thanks to its unique parallel architecture, SHAMAN allows at the same time to: i) explore the conformational landscape of RNA with atomistic explicit-solvent molecular dynamics (MD) simulations driven by state-of-the-art forcefields and ii) identify potential small-molecules binding sites in an efficient way with the aid of probes and the metadynamics [START_REF] Laio | Escaping free-energy minima[END_REF] enhanced-sampling technique. SHAMAN was benchmarked on a set of biologically relevant target systems, including large, structured riboswitches as well as smaller highly dynamic RNAs involved in viral proliferation. Our method successfully identified all the experimentally resolved pockets present in our benchmark set and was able to rank them among the most favorite probe hotspots. Our work constitutes to date the most advanced computational pipeline for binding site identification in dynamic RNA structural ensembles, thus providing crucial information for structurebased rational design of novel compounds targeting RNA.

Results

This section is organized as follows. First, we provide a general overview of SHAMAN and illustrate its accuracy in identifying experimentally resolved binding sites in a set of biologically relevant RNA targets. Second, we focus on the probes used in our SHAMAN simulations and investigate their relation to physico-chemical features of both the RNA pockets and the small molecules bound to them in known experimental structures. Finally, we present two case studies, the FNM riboswitch and the HIV-1 TAR, to i) demonstrate how SHAMAN can be used to study well-structured as well as more flexible RNAs;

ii) highlight the main strengths of our technique in modeling both local and global flexibility of the target.

Overview of the SHAMAN approach

SHAMAN is a computational technique that uses small fragments or probes and atomistic explicitsolvent MD simulations to identify potential small-molecule binding sites in RNA structural ensembles (Fig. 1A). SHAMAN is based on a unique architecture in which multiple replicas of the system are simulated in parallel (Fig. 1B). A mother simulation, containing only RNA and possibly structural ions, explores the conformational landscape of the target and communicates the positions of the RNA atoms to the replicas. Each replica contains a different probe that explores the RNA conformation provided by the mother simulation using the metadynamics enhanced-sampling approach [START_REF] Laio | Escaping free-energy minima[END_REF] . Soft positional restraints applied to the RNA backbone atoms of the replica allow for local induce-fit effects caused by the probes, while following or "shadowing" the conformational changes of the mother RNA simulation. This parallel architecture enables an efficient exploration of the same RNA conformation by different probes and to build, for each representative cluster of RNA conformations, a set of potential small-molecule binding sites or SHAMAPs (Fig. 1C). Each SHAMAP corresponds to a region of space occupied with high probability by at least one probe and is ranked by the binding free energy Δ𝐺 of the probe(s) to a specific RNA conformation (Fig. 1D). A more detailed description of SHAMAN is provided in Materials and Methods.

Benchmark of the SHAMAN accuracy

The accuracy of SHAMAN in identifying experimentally resolved binding sites was evaluated on 7 biologically relevant systems, including riboswitches (Fig. 2A) and viral RNAs (Fig. 2B). For each system, SHAMAN simulations were initialized from both holo conformations after the removal of the ligand (holo-like) and, when available, apo conformations, resulting in a total of 12 runs (Tab. S1 and S2). The validation set was composed of 14 unique binding pockets obtained from 69 experimental structures of riboswitches (Tab. S3) and viral RNAs in complex with different ligands (Tab. S4). For each simulation, the accuracy was defined in terms of the distance between our SHAMAPs and the ligand position in the reference experimental structures (Eq. 10 and Fig. 2C).

SHAMAN was able to identify the experimentally resolved pockets in all the systems of our benchmark set, both when initializing the simulations from holo-like and apo conformations (Tab. S5 and S6). Most importantly, the experimental binding sites were ranked among the most probable SHAMAPs in each corresponding run. To quantify the rank, we defined the difference in binding free energy ΔΔ𝐺 between each SHAMAP and the one with lowest free energy (Eq. 9). When starting from the apo conformation of the target RNA, the ΔΔ𝐺 of the SHAMAPs overlapping with the ligands was in 80% of cases below 𝑘 ! 𝑇 and in the 100% of cases below 2𝑘 ! 𝑇 (Fig. 2D). When starting from holo-like conformations, these percentages dropped to 64% and 84% (Fig. 2D). Ranking the experimental binding pockets among the SHAMAPs with lowest free energy (top scored) is fundamental in the context of CADD, and in particular in virtual screening applications (Discussions).

The geometrical proximity of our SHAMAPs to the experimental binding sites present in our benchmark set was noteworthy. The average distance between the centers of the interacting sites overlapping with a ligand and its position in the experimental structure was equal to 3.8 Å and 4.4 Å in the holo-like (Fig. 2E, upper panel) and apo (Fig. 2E, lower panel) cases, respectively. Both values are relatively small when compared to the distance threshold used in our validation criterion (Eq. 10), which was defined as the sum of the radius of gyration of the SHAMAP (on average ~1.6 Å, Fig. S1A) and the ligand (on average ~3.7 Å, Fig. S1B). As expected, this proximity to the experimental binding sites was remarkably greater in the simulations initiated from holo-like conformations in which the binding sites were already present. As a matter of fact, 22% of the successful interacting sites identified in the holo-like simulations were close to the experimental pocket by half of our distance threshold, while this holds only for 1% of the apo simulations.

Analysis of the probes

Two sets of probes were used in the SHAMAN benchmark described in the previous section. The first set of 8 probes (Tab. S7) was previously used in the development of SILCS-RNA [START_REF] Kognole | SILCS-RNA: Toward a Structure-Based Drug Design Approach for Targeting RNAs with Small Molecules[END_REF] and was mostly composed of aliphatic compounds selected to represent specific types of interaction with the RNA target. This set includes: acetate (ACEY), benzene (BENX), dimethyl-ether (DMEE), formamide (FORM), imidazole (IMIA), methyl-ammonium (MAMY), methanol (MEOH), and propane (PRPX). A second set of 5 probes (Tab. S8) was generated in this work using a fragmentation protocol (Materials and Methods) applied to ligands present in i) the HARIBOSS [START_REF] Panei | HARIBOSS: a curated database of RNA-small molecules structures to aid rational drug design[END_REF] database of RNA-ligand resolved structures (https://hariboss.pasteur.cloud); and ii) the R-BIND [START_REF] Donlic | R-BIND 2.0: An Updated Database of Bioactive RNA-Targeting Small Molecules and Associated RNA Secondary Structures[END_REF] database of bioactive small molecules targeting RNA (https://rbind.chem.duke.edu). This second set includes mostly aromatic compounds: benzene (BENX), dihydro-pyrido-pyrimidinone-Imidazo-pyridine (BENF), benzothiophene (BETH), methylpyrimidine (MEPY), and the cyclic non-aromatic piperazine (PIRZ).

We first explored the relation between the probes that successfully identified experimental binding sites and some of the structural features of RNA pockets. Aromatic probes showed a preference for exploring cavities buried deep inside the RNA structure (Fig. 3A, dark green bars), with an estimated average buriedness of 0.75 ± 0.06, which is relatively high compared to known RNA-small molecules pockets (Fig. 3B). On the other hand, non-aromatic or aliphatic probes displayed two distinct patterns. FORM, MEOH, and MAMY selectively explored shallow pockets with an average buriedness of 0.59 ± 0.04 (Fig. 3A, olive green bars), while DMEE, PRPX and ACEY promiscuously explored pockets with varying solvent exposure and an average buriedness of 0.70 ± 0.08 (Fig. 3A). PIRZ exhibited an intermediate behavior, with an average buriedness of 0.65 ± 0.06 (Fig. 3A, brown bar). As a consequence, aromatic probes were particularly successful (66% of cases) in identifying riboswitches binding sites, which in our validation set typically resided in buried cavities (Fig. 3C). For example, the location of the representative riboswitch binder GNG (PDB 3ski) was exclusively identified by aromatic probes (Fig. 3D). On the other hand, aliphatic probes identified pockets with high likelihood (70%) in viral RNAs (Fig. 3E), whose inherent flexibility resulted in shallow cavities exposed to solvent. An example is the binding site of SS0, a typical viral RNA binder (PDB 3tzr), which was identified primely by aliphatic probes (Fig. 3F).

Although the main goal of SHAMAN is pocket identification, motivated by its perspective use in virtual screening and ligand optimization (Discussion) we also investigated the link between the similarity of a given probe to a ligand and its ability to identify the corresponding experimental pocket. We started by comparing standard physico-chemical properties of the entire ligand or the corresponding Murcko scaffold (Materials and Methods). Our analysis did not reveal a strong correlation between ligands and probes (Tab. S9). We then calculated the Tanimoto similarity using different fingerprints (Materials and Methods). Our analysis suggested that we cannot predict whether a probe would be successful based on its similarity with a ligand (Fig. S2). However, based on a statistical classification (Materials and Methods), we can conclude that probes that did not resemble the ligand were highly unlikely to successfully identify the corresponding binding site, with a negative predictive value (NPV) equal to 0.82 (Eq. 11 and Tab. S10).

The case of the FMN riboswitch

The Flavin MonoNucleotide (FMN) riboswitch is an RNA molecule found in bacteria that regulates FMN gene expression via binding the FMN metabolite [START_REF] Wilt | FMN riboswitch aptamer symmetry facilitates conformational switching through mutually exclusive coaxial stacking configurations[END_REF] . Being the target of ribocil [START_REF] Howe | Selective small-molecule inhibition of an RNA structural element[END_REF] , one of the few FDA-approved compounds targeting RNA, the FMN riboswitch constitutes a natural test case for our drug-design purposes. As of today, 19 X-ray structures of the FMN riboswitch are deposited in the PDB database, 3 in apo and 16 in holo conformations. The 9 unique small molecules resolved in the holo structures fall into three main families: the cognate FMN family, the synthetic ribocil family, and the tetracyclic DKM binder (Fig. S3). The ligands belonging to the FMN and ribocil families share a Ushaped conformation and occupy the same binding site, buried into the RNA structure within the junctional region of the six stems between the A-48 and A-85 bases (Fig. 4A). The DKM tetracyclic ligand exhibits instead a distinct binding mode [START_REF] Rizvi | Discovery of Selective RNA-Binding Small Molecules by Affinity-Selection Mass Spectrometry[END_REF] as it induces a flip in A-48 and stacks face-to-face between A-48 and G-62, resembling the apo form (Fig. 4B). We therefore challenged our SHAMAN approach to capture the local rearrangements of the FMN riboswitch and to identify both types of binding poses starting from a single static structure.

We tested SHAMAN starting from both holo-like (PDB 6dn3 [START_REF] Vicens | Structure-Activity Relationship of Flavin Analogues That Target the Flavin Mononucleotide Riboswitch[END_REF] ) and apo (PDB 6wjr 51 ) structures (Fig. 4CD). One major RNA cluster, including the initial conformations, was populated for 99% and 84% of the holo-like and apo trajectories. This limited conformational variability observed in our simulations is consistent with the structural variety resolved experimentally (Tab. S11), supporting the accuracy of the force field used in our SHAMAN simulations. In this predominant RNA structural cluster, our method successfully located the experimental binding sites (Fig. 4CD) with very high accuracy, in the best case with a discrepancy of only 1.5 Å and 1.7 Å in the holo-like and apo simulations, respectively (Tab. S5).

Moreover, the experimental pocket was ranked in both cases among the most probable SHAMAPs (Fig. 2D), with a ΔΔG (Eq. 9) of 0.04 kJ/mol and 0.08 kJ/mol, respectively (Tab. S5). These results are even more remarkable if we consider the buried character of the FMN riboswitch pocket, which made it difficult for the probes to access it and sample accurately. As discussed above (Fig. 3), most of the probes that successfully identified this buried pocket were aromatic, both in the holo-like (83%) and apo (75%) cases (Fig. 4G).

Notably, the two distinct binding modes of FMN and DKM ligands were identified with comparable accuracy in both runs starting from holo-like and apo conformations. Each of these starting conformations was representative of one single binding mode: in the holo-like structure, the A-48 basis faces A-85, while in the apo case it is flipped onto A-49. SHAMAN enabled the identification of both binding modes, including the one not present in the starting conformation, something not possible with algorithms based on static structures. This is highlighted by superimposing the SHAMAPs found in the holo-like and apo simulations to the corresponding starting structure (Fig. 4CD, insets). The detection of both binding modes was made possible by simulating different probes in parallel and allowing for induce-fit effects in the RNA conformation sampled by the mother simulation (Discussions). In the hololike case, the BENX and IMIA probes captured the tail of the FMN binder (left panel, Fig. 4F, black and green surfaces, respectively), while BENF and MEPY overlapped with the tetracyclic part of DKM (right panel Fig. 4F, orange and celeste surfaces, respectively). In the apo case, MEPY interacting site overlapped with both ligands, but the tetracyclic part of DKM was captured only by IMIA (Fig. 4G).

The case of HIV-1 TAR element

The HIV-1 Trans-activation response element (HIV-1 TAR) is a highly flexible, non-coding RNA molecule responsible for regulating HIV-1 gene expression through binding with Tat protein [START_REF] Harrich | A critical role for the TAR element in promoting efficient human immunodeficiency virus type 1 reverse transcription[END_REF][START_REF] Chavali | Face-time with TAR: Portraits of an HIV-1 RNA with diverse modes of effector recognition relevant for drug discovery[END_REF] .

Understanding its conformational dynamics is crucial for drug development but remains challenging due to the major structural changes occurring upon binding diverse partners [START_REF] Davidson | A Small-Molecule Probe Induces a Conformation in HIV TAR RNA Capable of Binding Drug-Like Fragments[END_REF][START_REF] Musselman | iRED Analysis of TAR RNA Reveals Motional Coupling, Long-Range Correlations, and a Dynamical Hinge[END_REF] . This conformational plasticity of HIV-1 TAR is reflected in the more than 20 resolved structures, primarily by NMR, alone or bound to different ligands in water-exposed cavities. Our validation set was composed of 5 holo structures bound to different small molecules with different binding modes (Fig. S4) in the groove between the bulge UCU and the apical loop CUGGGA (residues 23-25 and 30-35, Fig. 5A). This is a crucial region that also encodes the Tat protein binding site [START_REF] Krawczyk | Tertiary Element Interaction in HIV-1 TAR[END_REF] . One of these structures (PDB 2l8h) indicates the presence of a transient and functionally relevant pocket formed upon binding to the MV2003 small molecule [START_REF] Davidson | A Small-Molecule Probe Induces a Conformation in HIV TAR RNA Capable of Binding Drug-Like Fragments[END_REF] . Given its complex dynamics, HIV-1 TAR constitutes an important benchmark of the capabilities of SHAMAN to detect binding sites appearing upon global conformational changes of the target molecule.

We tested SHAMAN starting from two structures of HIV-1 TAR, one in holo-like (PDB 1uts [START_REF] Murchie | Structure-based Drug Design Targeting an Inactive RNA Conformation: Exploiting the Flexibility of HIV-1 TAR RNA[END_REF] ) and one in apo (PDB 1anr 60 ) conformation. Both simulations recapitulated the expected flexibility of the target by identifying multiple significantly populated structural clusters (Fig. 5BC). A significant portion of the SHAMAPs was in the major groove of HIV-1 TAR (Fig. 5BC) with a relatively high probability (ΔΔ𝐺 within 2𝑘 ! 𝑇). Among these, SHAMAN identified all the 5 experimental binding sites, even though the overall similarity of the RNA to the deposited structures was never below ~3 Å backbone RMSD (Fig. S5). The most accurate overlaps with the experimental ligands were obtained with SHAMAPs detected in conformations b and e in the holo case (Fig. 5D) and conformations a, c, and d (Fig. 5E) in the apo case, mostly by aliphatic probes (Fig. 5F). The geometric accuracy in identifying the binding sites was inferior compared to the FMN riboswitch, with an average distance between binding and interacting sites equal to 4.0 Å and 4.1 Å for the holo-like ad apo cases, respectively (Tab. S6). However, we consider this distance still acceptable given the high flexibility of the molecule and the shallow nature of the experimental binding sites.

Notably, SHAMAN was able to identify the cryptic binding pocket proposed by Davidson et al. [START_REF] Davidson | A Small-Molecule Probe Induces a Conformation in HIV TAR RNA Capable of Binding Drug-Like Fragments[END_REF] (orange residues in Fig. 3B of their paper). In our simulations, this site was detected in conformation e (orange residues in Fig. 5C) by the ACEY and MAMY probes (red and pink densities, respectively).

While in the work of Davidson et al. the cryptic pocket appeared upon MV2003 binding to HIV-1-TAR, here its detection was made possible by the ability of SHAMAN to describe large conformational changes of small RNAs and account for induce-fit effects of the probes (Discussions).

Discussion

Here we presented SHAMAN, a computational technique for small-molecule (SM) binding site identification in RNA structural ensembles based on all-atom MD simulations accelerated by metadynamics,. We benchmarked the accuracy of our approach using a set of known RNA-SM structures, which included large, stable riboswitches and smaller, highly flexible viral RNAs. SHAMAN was able to identify all the binding pockets observed in the experimental structures and rank them among the most favorable probe interacting hotspots, both when starting from holo-like and apo conformations of the target. The interacting sites found by the SHAMAN simulations initiated from holo-like conformations were closer to the experimental pockets than those found in the apo cases. However, in the latter case the SHAMAPs corresponding to experimental binding sites were still very accurate and ranked as the top scored interacting sites for the majority of systems. Furthermore, our predictions were more accurate in the case of rigid riboswitches, with the regions explored by the probes perfectly matching the experimental binding sites. The accuracy was still very satisfying also for viral RNA molecules considering their high flexibility.

SHAMAN emerges as the most advanced physics-based approach for binding site identification in RNA structural ensembles. A major limitation of existing CADD tools in this framework is the inadequate treatment of RNA flexibility. In these regards, SILCS-RNA [START_REF] Kognole | SILCS-RNA: Toward a Structure-Based Drug Design Approach for Targeting RNAs with Small Molecules[END_REF] represents the state-of-the-art computational techniques by modelling the flexibility of the target RNA using a mixed-solvent MD approach. However, the method proposed by the MacKerell group presents two important limitations.

First, it makes use of positional restraints on the RNA backbone atoms and therefore is not designed to detect cavities formed upon major conformational changes. Second, SILCS-RNA was tested only by starting the MD simulations from holo structures after the removal of the bound ligand, therefore restraining the RNA target in a conformation in which the binding site is already formed. On the contrary, SHAMAN has been designed to enable the identification of pockets in dynamic RNA conformational ensemble characterized by both local and global conformational changes. The FMN riboswitch case study highlights how the target RNA molecules simulated in the replica systems have enough freedom to undergo local rearrangements induced by the probes and ultimately to capture the two distinct binding modes observed in the experimental structures. Furthermore, the challenging case study of HIV-1 TAR demonstrates that cryptic pockets formed upon global conformational rearrangements [START_REF] Davidson | A Small-Molecule Probe Induces a Conformation in HIV TAR RNA Capable of Binding Drug-Like Fragments[END_REF] can also be successfully identified by SHAMAN.

Despite the potentialities discussed above, the current implementation of SHAMAN presents two important limitations. First, the unbiased MD simulation of the RNA target in the mother replica will hardly ever provide a comprehensive exploration of the conformational space at low computational cost. However, this might not be a severe limitation if the scope is to determine potential druggable sites in the proximity of the metastable holo-like and apo RNA conformations resolved experimentally. To achieve a more global conformational exploration, in the future we will accelerate sampling of the RNA target in the mother replica by using enhanced-sampling techniques distributed with the PLUMED library, where SHAMAN is also implemented. Another limitation of our approach resides in the accuracy of the RNA force fields used in our MD simulations. Despite tremendous progress [START_REF] Salsbury | Recent developments in empirical atomistic force fields for nucleic acids and applications to studies of folding and dynamics[END_REF] , the accuracy of molecular mechanics force fields for nucleic acids is still as high as for proteins. One way to effectively improve the underlying force field is to integrate experimental data into MD simulations.

A large variety of integrative approaches, often based on Maximum Entropy and Bayesian principles [START_REF] Bonomi | Principles of protein structural ensemble determination[END_REF] have been developed in the past 10 years to use ensemble-averaged experimental data, such as many NMR observables, to model accurate structural ensembles of dynamic proteins. These approaches have been more recently applied to the determination of RNA structural ensembles [START_REF] Salmon | A general method for constructing atomic-resolution RNA ensembles using NMR residual dipolar couplings: The basis for interhelical motions revealed[END_REF][START_REF] Bottaro | Conformational ensembles of RNA oligonucleotides from integrating NMR and molecular simulations[END_REF] and will be used in the future to improve the accuracy of the RNA ensembles determined by SHAMAN.

In the future we foresee multiple different applications of SHAMAN in the context of CADD, in particular in combination with virtual screening applications and fragment-based drug design. Here our approach was used only to identify binding sites occupied by ligands in experimentally resolved structures. In this process, we also detected potential alternative binding sites that were in many cases ranked among the top scored SHAMAPs. For example, in the case of the THF riboswitch, we identified a top scored SHAMAP at the center of the RNA molecule between helix P2 and P3 (Fig. 6). In this region, to our knowledge, no binders have been experimentally determined yet. In the future, we will attempt at experimentally validating this pocket and eventually targeting it in a virtual screening campaign. Even more exciting is the application of SHAMAN to novel targets for which a small molecule has not been found yet. In these regards, the fact that top scored SHAMAPs often corresponded to known binding sites will allow us to restrict virtual screening campaigns to a few localized regions.

Despite the fact that we did not find a strong correlation between successful probes and ligands, we believe that SHAMAN can provide some guidance to tailor the choice of small molecules for virtual screening or to optimize known ligands. For example, in the case of riboswitches characterized by buried cavities and viral RNA with shallower and more exposed cavities, the results of our analysis suggested the use of molecules rich in aromatic or aliphatic moieties, respectively. In addition, areas close to the location of known ligands identified by certain probes as strong interacting hotspots could provide insights about how to modify the ligand to improve its affinity or even clues about ligand binding pathways (Fig. S6).

In conclusion, our method provides a novel and promising foundation for future drug design efforts targeting RNA. The accuracy, reliability, and versatility of SHAMAN in identifying small-molecule binding sites across diverse RNA systems with various degree of flexibility highlight its potential value in the field. By integrating SHAMAN in virtual screening pipelines, we aim in the future at creating an advanced platform for the rational in silico design of RNA-targeting molecules, effectively addressing the longstanding challenges in the field.

Materials & Methods

Details of the SHAMAN algorithm

SHAMAN consists of four main stages, each one composed of a set of operations described in detail in the following sections. At the beginning of each stage, we provide a brief non-technical overview to facilitate the reading.

I. Input stage

The initial input of SHAMAN consists of the 3D structures of the target RNA and of a set of N probes.

Starting from this information, we generate a reference mother system, including the RNA and possibly structural ions, and N replicas, each one with the addition of a different probe.

Setup of the mother simulation. The 3D structures of all the systems (Tab. S1) were obtained from the PDB database [START_REF] Berman | The Protein Data Bank[END_REF] . In the case of RNA structures determined by NMR, the first model was selected. In case of holo structures, the ligand was removed. Furthermore, to correctly model the RNA with our forcefield, the following elements were also eliminated, if present: crystal waters, PO3 group in the 3' terminal, modified residues at both terminals, and ions not modeled by our forcefield (SO4 in PDB 3tzr, 3ski and 7kd1). The resulting model was then prepared by adding hydrogen atoms using UCSF Chimera 65 at pH=7.4 and processed by the OpenMM library 66 v. 7.7.0 to generate an initial configuration and topology files. The forcefield used for RNA was AMBER99SB-ILDN* 67 with the BSC0 correction on torsional angles [START_REF] Pérez | Refinement of the AMBER Force Field for Nucleic Acids: Improving the Description of α/γ Conformers[END_REF] and the 𝜒 "#$ correction on anti-g shifts [START_REF] Zgarbová | Refinement of the Cornell et al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles[END_REF] . Ions were modeled using the Joung and Cheatam parameters [START_REF] Joung | Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations[END_REF] with the Villa et al. correction for magnesium [START_REF] Allnér | Magnesium Ion-Water Coordination and Exchange in Biomolecular Simulations[END_REF] . Water molecules were modelled with the OPC force field [START_REF] Izadi | Building Water Models: A Different Approach[END_REF] . Forcefield parameters were obtained from https://github.com/srnas/ff.

Setup of the replica simulations. The 3D structures of the probe were generated as described in the section Details of the probes. One replica of the mother system was generated for each probe. A single probe conformer was generated using the RDkit python library v. 2022.3 and inserted in a random position and orientation, with maximum distance of its center of mass from the RNA atoms equal to 1.0 nm. The force field and topology of the probe were created with OpenFF Sage 2.0 [START_REF] Boothroyd | Development and Benchmarking of Open Force Field 2.0.0: The Sage Small Molecule Force Field[END_REF] .

General details of the MD simulations.

Both mother and replica systems were solvated in a triclinic box with dimensions chosen in such a way each edge of the box was 1.0 nm away from the closest RNA atom. K+ and Cl-were added to ensure charge neutrality at salt concentration equal to 0.15 M. In all simulations the equations of motion were integrated by a leap-frog algorithm with timestep equal to 2 fs. The smooth particle mesh Ewald 74 method was used to calculate electrostatic interactions with a cutoff equal to 0.9 nm. Van der Waals interactions were gradually switched off at 0.8 nm and cut off at 0.9 nm. All simulations were performed with GROMACS 75 v. 2021.5 equipped with a development version of PLUMED 76 (GitHub master branch).

II. Production stage

After independently equilibrating mother and replica systems, the SHAMAN simulation proceeds in parallel. The RNA in the mother simulation is freely evolving and the positions of the RNA backbone atoms are communicated to the replica systems. A restraint is added to the positions of the backbone RNA atoms in the replica systems to make sure that they follow like shadows the conformation sampled by the mother. To accelerate the exploration of the RNA surface, the sampling of the probe in the replica systems is enhanced by metadynamics.

Equilibration procedure. All systems were independently equilibrated before the production stage.

This procedure consisted of i) energy minimization with steepest descent; ii) a 10 ns-long equilibration in the NPT ensemble using the Berendsen barostat 77 at 1 atm; iii) a 10 ns-long equilibration in the NVT ensemble using the Bussi-Donadio-Parrinello thermostat 78 at 300K. During the last two steps, harmonic restraints with harmonic constant equal to 400 kJ/mol/nm 2 were applied to the positions of the RNA backbone as well as probe atoms.

SHAMAN simulations. The systems were simulated in parallel for 1 µs each. The following settings were implemented using PLUMED. First, at each MD step the position of the atoms of the RNA backbone in the mother system were communicated to all the replicas and the corresponding atoms were restrained to have a maximum RMSD of 0.2 nm from the mother configuration. Second, to accelerate the probe exploration of the RNA surface, we used metadynamics [START_REF] Laio | Escaping free-energy minima[END_REF] . As collective variables 𝑺(𝑹), we used the xyz coordinates of the center of mass of the probe, defined after aligning the atoms of the RNA backbone to the initial reference conformation using the FIT_TO_TEMPLATE action in PLUMED. The well-tempered variant of metadynamics [START_REF] Barducci | Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method[END_REF] was used with biasfactor equal to 10. Gaussians with initial height of 1.2 kJ/mol and width of 0.1 nm were deposited every 1 ps. Finally, we restrained the position of the center of mass of the probe to be at most 1.0 nm away from the closest RNA atom using an upper harmonic wall with intensity equal to 10000 kJ/mol/nm 2 .

III. Analysis stage

For each representative cluster of RNA conformations explored by SHAMAN, we i) identified the regions with high probe occupancy; ii) defined a set of potential interacting sites for each probe; iii) clustered together the sites found by all probes to create the final SHAMAPs.

Metadynamics reweighting.

We removed the effect of the metadynamics bias potential on the probe trajectories by calculating for each frame the unbiasing weight 𝑤 % as 80 :

𝑤 % ∝ 𝑒𝑥𝑝 𝑉 & (𝑺(𝑹 % ), 𝑡 ̅ ) 𝑘 ! 𝑇 (1)
where 𝑉 & (𝑺(𝑹 % ), 𝑡 ̅ ) is the well-tempered metadynamics potential accumulated at the end of the simulation 𝑡 ̅ and evaluated on the conformation 𝑹 % . All these operations were performed independently for each simulation using the driver utility of PLUMED.

RNA clustering.

We first concatenated all the trajectories of the mother and replica simulations, after removal of the probes, and fixed the discontinuities due to the periodic boundary conditions. We then clustered all the RNA conformations with the gromos algorithm 81 implemented in GROMACS using as metrics the RMSD calculated on the RNA backbone atoms with a cutoff of 0.3 nm. To reduce memory requirements, the clustering was first performed on a subset of frames (1 every 10) and then the excluded frames were assigned to the closest cluster using a python script based on the MDAnalysis library 82 v.

2.2.0. The cluster center was taken as the representative structure for each state. The cluster populations were calculated independently for the mother and each replica simulation and clusters populated less than 10% were discarded in the subsequent analysis.

Calculation of probe free energy maps. The following analysis was performed independently for each replica and probe system as well as for each RNA cluster. We first extracted from each trajectory the frames corresponding to the selected cluster and aligned all the conformations to the RNA backbone atoms of the cluster center. We then defined a grid in the 3D space with voxel size equal to 0.1 nm and computed for each voxel 𝑖𝑗𝑘 the corresponding probe binding free energy 𝛿𝐺 '() as:

𝛿𝐺 '() = -𝑘 ! 𝑇 log 𝑁 '() 𝑁 * ( 2
)
where 𝑘 ! 𝑇 = 2.494339 kJ/mol and 𝑁 '() is the sum over all probe atoms of the (normalized) metadynamics unbiasing weights (Eq. 1) of the frames in which that atom explored the voxel 𝑖𝑗𝑘. 𝑁 * is the probe occupancy in the bulk solvent:

𝑁 * = 𝑛 +,-./ 𝑉 0-1/2 𝑉 34 (3)
where 𝑛 +,-./ is the number of probe atoms, 𝑉 0-1/2 and 𝑉 34 the volume of the voxels and simulation box, respectively. 𝛿𝐺 '() quantifies the propensity of finding a probe atom within the voxel 𝑖𝑗𝑘 rather than in the bulk solvent: voxels with low value of 𝛿𝐺 '() represent therefore potential strong binding sites to the RNA molecule. We estimated the associated error 𝜎 & by calculating the standard deviation of 𝛿𝐺 '() calculated in the first and second half of the trajectory (Fig. S7).

Voxels selection, clustering into interacting sites, and filtering.

To exclude weak affinity regions and independently for each probe, we first selected all the voxels within 10 kJ/mol from the minimum value of 𝛿𝐺 '() across all voxels. The selected voxels were then clustered into interacting sites using the DBSCAN algorithm implemented in the scikit python library 83 v. 1.8.1, with a maximum distance between points equal to 0.2 nm and a minimum number of samples equal to 5. For each interacting site, we calculated the associated binding free energy ∆𝐺 2 :

∆𝐺 2 = -𝑘 ! 𝑇 log L 𝑝 '() '() (4)
where 𝑝 '() = exp [-

5& !"# ) $ 6
] and the sum is over all the voxels belonging to the site. For each interacting site, we also defined its center 𝒈 2 as the free-energy weighted average position of the voxel centers 𝒓 '() :

𝒈 2 = ∑ 𝑝 '() 𝒓 '() '() ∑ 𝑝 '() '() (5)
and a free-energy-weighted radius of gyration 𝑅 2 as:

𝑅 2 = V ∑ W𝑝 '() ⋅ 𝑑Z𝒓 '() , 𝒈 2 [ 7 \ '() ∑ 𝑝 '() '() (6)
where 𝑑 is the Euclidean distance. Finally, we calculated the buriedness score 𝑥 .8,

2

of an interacting site to quantify its exposure to solvent. For each voxel 𝑖𝑗𝑘, we first defined the RNA density 𝑁 '() 9:; as the sum of the metadynamics unbiasing weights (Eq. 1) of the frames in which an RNA atom explored the voxel 𝑖𝑗𝑘. We then defined 𝑥 .8, 2 as:

𝑥 .8, 2 = 100 𝑁 2 L 𝑁 '() 9:; '() (7)
where the sum runs over all the 𝑁 2 voxels at the surface of the interacting site. Interacting sites with low buriedness score correspond to regions surrounded by few RNA atoms, i.e. exposed to solvent. All the sites with buriedness score lower than 0.15 were filtered out.

Calculation of the final SHAMAPs. For each representative cluster of RNA conformations, we defined a set of SHAMAPs by clustering together all the interacting sites found by all probes. To perform this operation, we used the DBSCAN algorithm applied to the centers of the interacting sites 𝒈 2 , with maximum distance between points given by 2 * [𝑅 2 ^+ 𝜎 9 ], where 𝑅 2 ^ is the average radius of gyration across all sites and 𝜎 9 their standard deviation, and a minimum number of samples equal to 1. For each SHAMAP, we defined the binding free energy ∆𝐺 < as the minimum free energy over all the interacting sites that clustered into this SHAMAP:

Δ𝐺 < = min 2∈< {∆𝐺 2 } (8)
and Δ∆𝐺 < has the difference between the binding free energy of a SHAMAP and the minimum value across all SHAMAPs (top scored):

Δ∆𝐺 < = ∆𝐺 < -min < {∆𝐺 < } (9)
IV.

Output stage

The SHAMAPs obtained at the end of the previous stage constituted the final set of hotspots associated to a given conformational state of the RNA target. The SHAMAPs are reported in a table and ordered by Δ𝐺 < . Along with this information, each SHAMAP is annotated with the properties of its constituent interacting sites: a list of probes that explored the region, their correspondent ∆𝐺 2 , the population of the RNA cluster in which the site has been visited, the coordinates of the centers 𝒈 2 and the radius of gyration 𝑅 2 .

Details of the SHAMAN benchmark

Details of the target RNAs. For our SHAMAN simulations, we selected 7 RNA systems, whose structures in complex with at least one ligand were deposited in the PDB databank [START_REF] Berman | The Protein Data Bank[END_REF] (Tab. S1). To initiate the simulations, we selected 1 holo structure per system and, when available, an apo structure of the same RNA molecule. In total we performed 12 SHAMAN simulations. A summary of all simulations performed along with details about the systems are reported in Tab. S2.

Details of the PDB structures used for validation.

To benchmark the accuracy of our approach, we first retrieved for each system all the holo structures deposited in the PDB with different ligands and binding poses. We then visually inspected each structure and identified 14 structures with unique binding poses and pockets. All the structures used for validation along with details about the RNA, the ligand, and the experimental method and resolution are reported in Tab. S3 and S4.

Details of the probes. In our SHAMAN simulations we used two set of probes. The first set (Tab. S7) was taken from the SILCS-RNA study [START_REF] Kognole | SILCS-RNA: Toward a Structure-Based Drug Design Approach for Targeting RNAs with Small Molecules[END_REF] and for the most part originally developed for proteins. It is (https://hariboss.pasteur.cloud) and RBind 26 (https://rbind.chem.duke.edu) databases, respectively. The molecules were then prepared using Epik 85 at pH 7.4. All molecules were converted into canonical SMILES using RDkit v. 2022.3 and then processed using the RDkit Murcko Scaffolds KNIME node.

The two sets were compared using the 'Compare Ligand Sets' node provided by Schrodinger v. 2022.3 and 6 Murcko scaffolds were identified as common. Finally, the RECAP fragmentation method [START_REF] Lewell | RECAPRetrosynthetic Combinatorial Analysis Procedure: A Powerful New Technique for Identifying Privileged Molecular Fragments with Useful Applications in Combinatorial Chemistry[END_REF] (implemented as the 'Fragments from Molecules' node provided by Schrodinger) was employed and resulted in 6 fragments (Tab. S8). From these fragments, a selection of 5 representative ones was used for the SHAMAN simulations: benzene (BENX), dihydro-pyrido-pyrimidinone-imidazo-pyridine (BENF), benzothiophene (BETH), methyl-pyrimidine (MEPY), and the cyclic non-aromatic piperazine (PIRZ). BETH was intentionally modeled in a protonated state, as it appears as part of the molecules extracted from RBind and HARIBOSS rather than as a fragment alone.

Details of the validation procedure. To benchmark the accuracy of our approach in identifying binding sites occupied by a ligand in known experimental structures, we used the following procedure:

i.

Multiple sequence alignment

For each simulated system, we aligned the sequence of our target RNA with the sequences of all the validation PDBs using CLUSTALW 87 v. 2.0.

ii.

Structural alignment of validation PDBs to SHAMAN cluster centers

For each validation PDB, we defined the binding site as the set of nucleotides with at least one atom within 0.6 nm of a ligand atom. The backbone atoms of the validation PDB belonging to this region were then structurally aligned to the corresponding nucleotides in each RNA cluster center, based on the sequence alignment defined above.

iii.

Definition of success for a probe interacting site

For each validation PDB, we defined an experimental sphere centered on the center of mass of the heavy atoms of the ligand 𝒈 /1+ and with a radius given by its radius of gyration 𝑅 /1+ . For each probe interacting site, we defined a validation sphere centered on the free-energy weighted center of the interacting site 𝒈 2 and with radius given by its free-energy weighted radius of gyration 𝑅 2 . We then considered a probe interacting site as successful if the validation sphere was overlapping with the experimental sphere:

𝑑Z𝒈 2 , 𝒈 /1+ [ ≤ 𝑅 2 + 𝑅 /1+ (10)
In case of match with multiple validation structures, we retained only the one corresponding to the interacting site with lower Δ∆𝐺 from the top scored SHAMAP.

iv.

Definition of success for a SHAMAP

A SHAMAP was considered successful in identifying a known ligand binding site if at least one of the probe interacting sites that compose the SHAMAP was successful according to the criterion defined above.

Probes-ligands comparison

For probes and ligands in the SHAMAN simulations initiated from holo structures, we first calculated the following set of descriptors with RDKit v. 2022.3: molecular weight, number of aromatic rings, number of H-bond donors/acceptors, topological polar surface area (TPSA), and number of heterocycles. The correlation between probes and ligands descriptors was then computed with scipy v.

1.8.1 using the Pearson correlation coefficient. The analysis was performed using either the entire ligand or its Murcko scaffold. We also quantified the similarity between ligands and successful probes using different types of fingerprints (FPs) implemented in RDKit. In particular, we used Morgan (radius = 2, 2048 bits), RDKit (2048 bits), and MACCS FPs . Using these FPs and the Tanimoto distance, we calculated the similarity between successful probes and reference ligands, considered either as entire ligands or using their corresponding Murcko scaffold.

To further investigate a possible correlation between ligand and successful probes, we formulated the following hypothesis: the ability of a probe to identify a binding site is related to its similarity to the corresponding ligand. We then compared each of the 13 probes (Tab. S7 and S8) with all the 8 ligands resolved in the experimental pockets (Tab. S1) and considered a probe to be similar (dissimilar) to a ligand if the Tanimoto distance calculated with MACCS FP was greater (lower) than 0.4 (0.2). Based on the SHAMAN results in our benchmark, we built a confusion matrix of the four possible outcomes (Tab. S10) and defined the SHAMAN negative predictive value 𝑁𝑃𝑉 as the ratio between true negatives TN and total number of negatives TN+FN:

𝑁𝑃𝑉 = 6:
6:>?: (Eq. 11)

Software and data availability

SHAMAN simulations can be run with the development version (GitHub master branch) of PLUMED (https://github.com/plumed/plumed.github.io). The GROMACS topology files and PLUMED input files used in our benchmark are available on PLUMED-NEST (www.plumed-nest.org), the public repository of the PLUMED consortium, as plumID:23.031. Scripts to facilitate the preparation of the input files and the analysis of the results as well as a complete tutorial will be released soon under a license "free for academics, not for commercial use". For each SHAMAP, we provide the binding free energy to RNA (𝛥𝐺) and the difference with respect to the lowest free energy (top scored) SHAMAP (𝛥𝛥𝐺) along with a list of probes that explored the corresponding regions.

Z 3 L A z c n q O 1 v r 1 r I w T o 6 R k 6 A s E = " > A A A C D X i c b V D L S s N A F J 3 U V 6 2 v q E s 3 g 1 V w Y 0 n E 1 0 Y o V d B V q d A X N C V M p p N 2 6 C Q T Z i Z C C f k B N / 6 K G x e K u H X v z r 9 x 2 m a h r Q c u H M 6 5 l 3 v v 8 S J G p b K s b y O 3 s L i 0 v J J f L a y t b 2 x u m d s 7 T c l j g U k D c 8 Z F 2 0 O S M B q S h q K K k X Y k C A o 8 R l r e 8 H r s t x 6 I k J S H d T W K S D d A / Z D 6 F C O l J d c 8 c G 4 I U w j e w q v j o V u B d Q g d x v v Q 8 Q X C S T V N q q 6 V u m b R K l k T w H l i Z 6 Q I M t R c 8 8 v p c R w H J F S Y I S k 7 t h W p b o K E o p i R t O D E k k Q I D 1 G f d D Q N U U B k N 5 l 8 k 8 J D r f S g z 4 W u U M G J + n s i Q Y G U o 8 D T n Q F S A z n r j c X / v E 6 s / M t u Q s M o V i T E 0 0 V + z K D i c B w N 7 F F B s G I j T R A W V N 8 K 8 Q D p H J Q O s K B D s G d f n i f N k 5 J 9 X j q 7 P y 2 W K 1 k c e b A H 9 s E R s M E F K I M 7 U A M N g M E j e A a v

  composed of: acetate (ACEY), benzene (BENX), dimethyl-ether (DMEE), formamide (FORM), imidazole (IMIA), methyl-ammonium (MAMY), methanol (MEOH), and propane (PRPX). Starting from the 2D structure of each probe, we prepared a minimized 3D structure at the protonation state corresponding to pH 7.4 using the LigPrep module 84 of Schrodinger Suite. A second set of RNA-specific probes was generated using the following KNIME 4.6 protocol. We first extracted 265

Figure 1 .

 1 Figure 1. Overview of the SHAMAN approach. A) Input stage: Selection of the RNA target structure and of the probes to initialize the mother and replica systems, each one with a different probe. B) Production stage: the unbiased/unrestrained MD simulation of the mother system communicates the positions of the RNA backbone atoms to the replicas, which are restrained to follow the mother like shadows. The probe exploration of the RNA conformation is accelerated by metadynamics. C) Analysis stage (from top to bottom): i) the sampled RNA ensemble is clustered into a set of representative conformations; ii)for each cluster and probe, a free-energy map is calculated from the probe occupancy during the course of the simulation; iii) voxels in the free-energy maps are clustered together into interacting sites; iv) for each interacting site, free energy and buriedness score are calculated and sites too exposed to solvent are discarded; v) for each RNA cluster, all interacting sites obtained from all probes are clustered together into SHAMAPs. D) Output stage: two RNA representative clusters with population equal to 32% (light brown, left panel) and 28% (pink, right panel) with the corresponding SHAMAPs (green circles).

  kJ/mol ] < l a t e x i t s h a 1 _ b a s e 6 4 = " l 0

0 Figure 2 .Figure 3 .Figure 4 .

 0234 Figure 2. Assessment of the SHAMAN accuracy. A)A cartoon-surface representation of the four riboswitches in our benchmark set (Tab. S1), with the corresponding name in the upper left of each panel. In the lower part, the PDB id of the starting structure used in our SHAMAN simulations is reported in a brown and cyan box for the holo-like and apo case (when available), respectively. The cartoon representations correspond to the holo-like structures. B) As in panel A), for the three viral RNAs of our benchmark set (Tab. S1). C) Definition of the validation distance (Eq. 10) as the distance between the free-energy weighted center of an interacting site and the center of mass of the experimental ligand. D) 𝛥𝛥𝐺 distribution of the probes that correctly identified known experimental pockets for holo-like (brown) and apo simulations (cyan). E) Scatter plots of the validation distance (x axis) and cutoff defined by Eq. 10 (y axis) for holo-like (brown, upper panel) and apo (cyan, lower panel) simulations. The dashed line indicates validation distances equal to the validation cutoff, while the dotted line corresponds to half the validation cutoff. Each system is identified by a different marker shape, as defined in the legend.

Figure 5 .FFigure 6 .

 56 Figure 5. The case of the HIV-1 TAR. A) 2D structure of the HIV-1 TAR. The two stem regions are indicated in light grey; the bulge (residues 23-25) and the apical loop (residues 30-35) in black. B-C) Representative RNA clusters determined by the SHAMAN simulations initiated from the holo-like (B) and apo (C) conformations. SHAMAPs are visualized as solid surfaces with the color code defined in Tab. S7 and S8. The RNA state labeled as "conf e" in panel C is represented as a grey surface to highlight the orange region explored by ACEY (red density) and MAMY (rose density). This area corresponds to the cryptic binding site identified by Davidson et al. 56 . D-E) Representative RNA conformations and SHAMAPs with best overlap with the experimental binding sites found in the simulations initiated from the holo-like (D) and apo (E) conformations. In the insets, SHAMAPs that best identified the 5 ligands present in our validation set (Tab. S4): clockwise from top left, ARG in PDB 1arj, PMZ in PDB 1lvj, P13 in PDB 1uts, P12 in PDB 1uui, MV2003 in PDB 2l8h. F) 2D structures of the probes that successfully identified the experimental binding sites. The brown and cyan dashed circles indicate the successful probes in the holo-like and apo simulations, respectively.

Acknowledgement F. P. Panei was funded by Sanofi and the Association Nationale de la Recherche et de la Technologie (ANRT) contract 2020/1259. This work was granted access to the HPC resources of IDRIS under the allocation 2022-AD01101371 made by GENCI.

Ethics declarations

F. P. Panei and P. Gkeka are or were Sanofi employees and may own stocks in Sanofi. M. Bonomi declares no competing interests.