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

Abstract

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 simulations 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 identifies all the experimentally resolved pockets and ranks 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.

Fig. 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). For each SHAMAP, we provide the binding free energy to RNA ($\Delta G$) and the difference with respect to the lowest free energy (top scored).

Fig. 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 $\Delta \Delta G$ 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.

Fig. 3. Analysis of the SHAMAN probes.

A Violin plots representing the buriedness of the experimental pockets (y-axis) successfully identified by a given SHAMAN probe (x-axis). Buriedness values were extracted from the HARIBOSS database (Tab. S3 and S4, sample size n = 420). Each box represents the interquartile range between the first and third quartiles, with the median indicated by a horizontal black line. Outliers are marked as black diamonds. B Buriedness distribution for the RNA pockets occupied by ligands in all the structures deposited in HARIBOSS. C Total number of times that a probe explored an experimental binding site in the riboswitches of our validation set. D Cartoon representation of the 2’-deoxyguanosine (dG) riboswitch (PDB 3ski) with 2D structure of the GNG binder. In the dashed box, the 2D structures of the probes that identified the GNG binding site. E As in (C), for the viral RNAs of our validation set. F Cartoon representation of the RNA from the Hepatitis C Virus (PDB 3tzr) with 2D structure of the SS0 binder. In the dashed box, the 2D structures of the probes that identified the SS0 binding site.

Fig. 4. Comparison with other tools.

From left to right, boxplots reporting the predictive quality of different binding site prediction tools evaluated by four statistical metrics for binary classifiers (Methods). A Binding site prediction on the holo-like systems (Tab. S1, red column, sample size n = 7) validated against the single corresponding experimental structure. B, C Binding site prediction on holo-like (B) and apo (C) systems (Tab. S1, red and cyan columns) against all the validation structures (Tab. S3 and S4, Methods, sample size n = 69). Each box represents the interquartile range between the first and third quartiles, with the median indicated by a horizontal black line. Outliers are marked as black diamonds.

Fig. 5. The case of the FMN riboswitch.

A Key RNA binding site residues for the FMN ligand (PDB 2yie) and ribocil (PDB 5kx9) families. B Key RNA binding site residues for the DKM ligand (PDB 6bfb) and in the apo conformation (PDB 6wjr). C, D Cartoon representation of holo-like (C) and apo (D) starting structures used in the SHAMAN simulations of the FMN riboswitch. In the insets, the key binding site residues are overlayed with the probe densities (colors as in Tab. S7 and S8). E 2D structures of the probes that successfully identified the experimental binding sites in the FMN riboswitch. The brown and cyan dashed circles indicate the successful probes in the holo-like and apo simulations, respectively. F, G For the holo-like (F) and apo (G) simulations, the SHAMAPs with best overlap with FMN (left) and DKM (right) ligands, representing the two different binding modes of the FNM riboswitch.

Fig. 6. The case of the HIV-1 TAR.

A 2D structure of the HIV-1 TAR. The two stem regions are indicated in light gray; 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 gray 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. . 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.

Fig. 7. Identification of an alternative pocket in the THF riboswitch.

In the upper panel, cartoon representation and molecular surface of the center of the most populated RNA cluster found in the SHAMAN simulation initiated from a holo-like conformation (PDB 4lvx). The THF riboswitch presents two binding pockets (dashed circles), one in a three-way junction (HB4 ligand bound between helical domains P2, P3 and P4, right side) and the other in a pseudoknot (HB4 ligand bound in PK region, left side). The experimental ligands in PDB 4lvx are superimposed by aligning the coordinates to the RNA cluster center. Our protocol detected a low free-energy SHAMAP in the middle of the THF riboswitch between helix P2 and P3 (surfaces surrounded by orange circle, colored as defined in Tab. S7 and S8),. In the lower panel, the light gray and light orange tables report the details of the SHAMAPs that identified the two experimental and the alternative binding sites, respectively.