Thermal titration molecular dynamics (TTMD): shedding light on the stability of RNA-small molecule complexes

Ribonucleic acids are gradually becoming relevant players among putative drug targets, thanks to the increasing amount of structural data exploitable for the rational design of selective and potent binders that can modulate their activity. Mainly, this information allows employing different computational techniques for predicting how well would a ribonucleic-targeting agent fit within the active site of its target macromolecule. Due to some intrinsic peculiarities of complexes involving nucleic acids, such as structural plasticity, surface charge distribution, and solvent-mediated interactions, the application of routinely adopted methodologies like molecular docking is challenged by scoring inaccuracies, while more physically rigorous methods such as molecular dynamics require long simulation times which hamper their conformational sampling capabilities. In the present work, we present the first application of Thermal Titration Molecular Dynamics (TTMD), a recently developed method for the qualitative estimation of unbinding kinetics, to characterize RNA-ligand complexes. In this article, we explored its applicability as a post-docking refinement tool on RNA in complex with small molecules, highlighting the capability of this method to identify the native binding mode among a set of decoys across various pharmaceutically relevant test cases.

1 Introduction

Following the “Central Dogma” of molecular biology, RNA has been historically perceived as a bridging element between DNA genetic information and protein biosynthesis (Bissaro et al., 2020; Pavan et al., 2022a). However, this paradigm shifted in the last decades due to discoveries that showcased RNA’s strikingly complex genetic and catalytic functions (Sponer et al., 2018; Childs-Disney et al., 2022). Therefore, this biomolecule is no longer considered just a carrier of information. Instead, it is perceived as one of the most pluripotent actors in molecular biology (Sponer et al., 2018).

Considering that only a fraction (3%) of the genome is translated into proteins (Zafferani and Hargrove, 2021) non-coding RNAs (ncRNAs) are under the spotlight as they are involved in epigenetics (Falese et al., 2021; Childs-Disney et al., 2022), playing a pivotal role in the etiopathogenesis of both cancer and neurodegenerative diseases (Falese et al., 2021).

The identification of RNAs as possible therapeutic targets and the increased availability of experimentally solved three-dimensional structures of RNA complexes paves the way for the application of structure-based drug design (SBDD) techniques to characterize the interaction between small molecules and RNA, allowing the possibility of rationally discover new hits and steer their development into mature leads (Pavan et al., 2022a; Childs-Disney et al., 2022). Although it is an appealing perspective, applying routinely adopted molecular modeling protocols to RNA systems is not trivial. Historically, these techniques have been optimized to study the recognition process between small organic molecules as ligands and proteins as receptors (Salmaso and Moro, 2018; Bassani and Moro, 2023; Pavan and Moro, 2023). Nevertheless, structural differences between proteins and RNAs, such as the peculiar surface charge properties portrayed by the polyanionic phosphate backbone, the ions’ role in the structural stability and folding of RNA, the role of the solvent in mediating structural stability and forming bridged interactions, other than the intrinsic structural flexibility of ribonucleic acids, limited so far the possibility to repurpose these methodologies to the study of RNA complexes (Bissaro et al., 2020).

So far, different computational techniques have been developed to tackle the study of RNA-small molecule complexes from different perspectives, from binding site identification to binding mode prediction, scoring, and characterization of time-dependent properties (Zhou et al., 2022).

One of the most successful approaches in structure-based drug design is molecular docking, a fast and reliable protocol that provides possible binding hypotheses generated through a conformational search algorithm and ranked through a scoring function (Meng et al., 2012; Ferreira et al., 2015). Albeit relatively fast to perform and computationally effortless when compared to other techniques like molecular dynamics, pose scoring can suffer from inaccuracies (Chaput and Mouawad, 2017). Moreover, docking neglects the role of the solvent and RNA’s plasticity, preventing the exploration of its conformational landscape (Bissaro et al., 2020; Pavan et al., 2022a).

To overcome the limitations of a docking approach, molecular dynamics simulations can be theoretically exploited to characterize RNA-ligand interactions with a more physically rigorous methodology. However, the computational effort required to spontaneously sample both binding and unbinding processes during unbiased MD simulations makes this technique incompatible with the timings of modern-day drug discovery campaigns (Bernardi et al., 2015).

Enhanced sampling algorithms are thus exploited to cut simulation times without altering the technique’s validity. These protocols rely on energetically biasing the system to increase the frequency of observation of the desired event and have been successfully applied to the study of unbinding processes in protein-ligand complexes (Do et al., 2013).

Many popular methods, including Steered Molecular Dynamics (Izrailev et al., 1999), Random Accelerated Molecular Dynamics (Lüdemann et al., 2000), and Umbrella Sampling (Torrie and Valleau, 1977), depend on the definition of so-called “collective variables” (CV), i.e., a set of descriptors that can be used to monitor the simulation and appropriately biasing the potential energy landscape. These CVs are difficult to identify since they are heavily system-dependent and rarely generally applicable (Salmaso and Moro, 2018). On the contrary, other tempering methods like Replica Exchange (Sugita and Okamoto, 1999), and Temperature Accelerated Molecular Dynamics (So/rensen and Voter, 2000), do not rely on collective variables but still require some tinkering for optimal sampling (Zamora et al., 2016).

Within the CV-free category, Thermal Titration Molecular Dynamics (TTMD) (Pavan et al., 2022c) is a recently developed method that addresses the increasing interest in the prediction of drug-target residence time, since kinetic properties such as the dissociation rate (koff) better correlate to in vivo ligand efficacy compared to thermodynamic properties like the equilibrium dissociation constant (Kd) (Copeland et al., 2006).

Contrary to the aforementioned techniques, the biggest perk of TTMD is its simplicity concerning simulation setup and trajectory analyses, making it more accessible to medicinal chemists without a strong modeling background. Initially developed for the comparison of protein-ligand complexes based on the persistence of their native intermolecular interaction, TTMD provides a simple and robust platform for ranking binding poses upon a defined receptor binding site (Menin et al., 2023) or classifying ligands based on their dissociation rate (Pavan et al., 2022c). Specifically, TTMD operates through a series of short classic molecular dynamics simulations performed at increasing temperatures while monitoring the persistence of native receptor-ligand interactions through an interaction fingerprints-based scoring function (Pavan et al., 2022b).

Due to the encouraging success observed in the characterization of protein-ligand complex stability with TTMD and based on the successful repurposing of various computational tools designed to work on proteins to the study of RNA targets (Mlýnský and Bussi, 2018), in the present work we aim to extend the applicability domain of the technique to the world of RNA-ligand complexes. In detail, we tried to understand if TTMD can be used as a post-docking filter to steer the identification of a native-like binding pose for small organic ligands onto RNA receptor binding sites.

2 Materials and methods

2.1 Hardware overview

Modeling tasks such as the structure preparation of RNA targets and respective ligands, docking, system setup for molecular dynamics (MD) simulations, and subsequent trajectory analysis were carried out on a Linux workstation running Ubuntu 20.04 as its operating system equipped with a 20 cores Intel Core i9-9820 × 3.3 Ghz processor. MD simulations were performed exploiting an in-house GPU cluster composed of 20 NVIDIA devices ranging from GTX1080Ti to RTX3090.

2.2 Structure preparation

The three-dimensional structures of the targets presented in this work were retrieved from the Protein Data Bank (PDB (Berman et al., 2000)), and successfully processed through different built-in modules of the Molecular Operating Environment (Molecular Operating Environment MOE, 2022.02, 2022). At first, each inconsistency between the primary sequence and the tertiary structure was fixed through the “Structure preparation” module. Then, the “Protonate3D″ tool was exploited for adding missing hydrogens according to the most probable tautomeric and protonation state of titratable groups at pH = 7.4. Finally, every non-RNA and non-ligand residue was removed, except for K+ ions in the G-quadruplex portion of 5BJO complex, which were retained as they play a pivotal role in the stabilization of these non-canonical structures that shape the binding site at the monomer-monomer interface (Bhattacharyya et al., 2016).

To avoid any bias in the docking calculation that may favor a crystal-like pose, the ligands have been prepared starting from their SMILES, exploiting Open Babel (O’Boyle et al., 2011) and various QUACPAC OpenEye tools (Cadence Molecular Sciences, 2023). Open Babel was used to generate two-dimensional and three-dimensional coordinates, while “Tautomers” and “Fixpka” assigned the correct protomeric and tautomeric state at pH 7.00. Afterwards “Molcharge” set the partial charges according to the MMFF94 force field. Both the RNA and ligand were stored for docking calculations.

2.3 Docking calculation

A self-docking calculation was conducted through the Protein-Ligand ANT System (PLANTS (Korb et al., 2006; Korb et al., 2007; Korb et al., 2009) program, one of the state-of-the-art docking tools for nucleic acids-ligand docking (Jiang et al., 2023). The binding site was defined as a sphere of radius 10.5 Å centered around the center of mass of the crystal ligand in the experimentally solved complex. The first 5 poses according to the ChemPLP scoring function were stored for refinement through the Thermal Titration Molecular Dynamics protocol.

2.4 System setup for MD simulations and equilibration protocol

Several packages of Visual Molecular Dynamics (VMD (Humphrey et al., 1996)) 1.9.3 and AmberTools22 (Case et al., 2005; Case, 2022) were used to prepare each RNA-ligand complex for MD simulations.

Each nucleic atom was parametrized according to ff14SB force with χ modification tuned for RNA (χ_OL3) field (Pérez et al., 2007; Zgarbová et al., 2011; Maier et al., 2015), while ligands parameters were assigned according to General Amber Force Field (GAFF (Wang et al., 2004)). Since the Corn Aptamer (PDB ID: 5BJO) presents a modified Uracil (I5-U17), a custom preparation was exploited following Amber’s workflow to create modified residues employing the antechamber and parmcheck2 tools, assigning parameters from ff14SB and GAFF force field.

Each RNA-ligand complex was solvated within a rectangular base prism box, with a 15 Å padding between the box border and the nearest solute atom, using the TIP3P water model (Jorgensen et al., 1983).

Sodium and chlorine monovalent ions were added to neutralize the net charge of the box, reaching the physiological salt concentration of 0.154 M. Before equilibration, each system was then subjected to 500 steps of energy minimization with the conjugate-gradient method to remove clashes and bad contacts.

Following this preparation phase, a two-step equilibration process was performed. The first phase consisted of a 0.5 ns simulation in canonical ensemble (NVT) imposing a 5 kcal mol⁻¹ Å⁻² harmonic positional restraint to each RNA and ligand atom, leaving water and ions unconstrained. In the second equilibration run, a 0.5 ns simulation was conducted in the isothermal-isobaric ensemble (NPT) constraining only the ligand and the backbone atoms of RNA through the same force field constant used in the first stage. For each equilibration stage, the temperature was kept constant at the lowest value indicated in the temperature range through a Langevin thermostat (Davidchack et al., 2009), and for the second NPT simulation, the pressure was fixed at 1 atm through a Monte Carlo barostat (Faller and de Pablo, 2002).

For each MD simulation, an integration timestep of 2 fs was used. The simulations were performed exploiting the proprietary ACEMD 3.5 engine (Harvey et al., 2009), which is based on the open-source library for molecular simulations OpenMM (Eastman et al., 2017). The M-SHAKE algorithm was used to constrain the length of bonds involving hydrogen atoms, the particle-mesh Ewald (PME (Essmann et al., 1995)) method was exploited to compute electrostatic interactions using cubic spline interpolation, and finally, a 9.0 Å cutoff was used for calculating Lennard–Jones interactions.

2.5 Thermal titration molecular dynamics (TTMD) simulations

Thermal Titration Molecular Dynamics (TTMD (Pavan et al., 2022c; Menin et al., 2023)) is an enhanced sampling molecular dynamics approach originally developed for the estimation of protein-ligand unbinding kinetics. This method consists of a series of short MD simulations (TTMD-steps) performed at progressively increasing temperatures. The length and the temperature for each step, defining the temperature ramp for the TTMD simulation, are defined by the user based on the knowledge of the system of interest, especially concerning the conservation of the native receptor fold throughout the whole simulation. In this work, two different temperature ramps were used: the first one, defined as “standard”, is the same ramp described in the original publication (starting temperature 300 K, final temperature 450 K, temperature increase 10 K, step length 10 ns), and a second one defined as “alternative” (starting temperature 73 K, final temperature 223 K, temperature increase 10 K, step length 10 ns). An interaction fingerprint-based scoring function (Pavan et al., 2022b) is exploited to monitor the conservation of the native binding mode. Specifically, the IFP_CS scoring function (Pavan et al., 2022b) is used to compare each of the binding features of each trajectory frame (encoded as protein-ligand interaction fingerprints) to the last frame of the second equilibration stage. In the context of this work, two different PLIFs were employed, specifically the InteractionFingerprint function of the Open Drug Discovery Toolkit (ODDT (Wójcikowski et al., 2015)) and the one provided by the ProLIF package (Bouysset and Fiorucci, 2021), while the cosine similarity metrics was used for the comparison, as implemented in the Scikit-learn Python library. The calculated value is then multiplied by −1 to comply with the scale of most scoring functions. In the end, the resulting score can range from −1 to 0, where −1 indicates total congruence of the PLIFs (total conservation of the binding mode) and 0 indicates that all native binding features are lost.

The simulation continues until the end of the temperature ramp is reached, or an early termination criterion is reached. Specifically, at the end of each TTMD step, the average IFPcs score is calculated for the last 10% of the step: if this value is above −0.05, the simulation is stopped as the ligand lost its original binding mode. The code for running TTMD simulations is open-source and available through the MIT license at github.com/molecularmodelingsection/TTMD.

2.6 Trajectory analysis, MS and IFF coefficients determination

TTMD trajectories were analyzed partially through the same Python code described in the previous paragraph and partially through the SuMD-analyzer (Pavan et al., 2022a) script available at github.com/molecularmodelingsection/SuMD-analyzer. Specifically, the MDAnalysis package (Michaud-Agrawal et al., 2011; Gowers et al., 2016) was exploited to calculate the RMSD of the receptor backbone, the ligand, and a defined set of binding site residues.

The “titration timeline” plot reports the time-dependent evolution of both the ligand, binding site and receptor backbone RMSD and of the IFP_CS score.

The “titration profile” plots the average IFP_CS score for each TTMD step against the temperature at which the step was executed. In this plot, the slope of the straight line linking the first and last point of the simulation (the MS coefficient) is extracted and used as a proxy measure for the estimation of the overall complex stability. The mathematical formulation of the MS coefficient is reported in Eq. (1).

$MS=\frac{{meanIFP}_{CS}^{Tend}-\left(-1\right)}{T^{end}-T^{start}}(1)$

MS formula.

The MS coefficient can range from zero (indicative of a tight and persistent binding/high conservation of the native binding mode) to 1 (indicative of high volatility of the native pose). For each investigated receptor-ligand complex, five independent TTMD simulations were performed, with the average MS coefficient being calculated across three different replicates, after discarding the highest and the lowest value.

Furthermore, to complement the MS coefficient, a second one defined as IFF was calculated as defined in Eq. (2).

$IFF=\sqrt{\frac{{\displaystyle \sum _t^T}{\left({IFP}_{CSt}-〈{IFP}_{CS}〉\right)}^2}{T}}(2)$

IFF formula.

The IFF coefficient is calculated as the root mean square fluctuation of the IFP_CS value across the simulation to the whole run mean IFP_CS value. As the MS coefficient only considers the initial and the final state of the simulation, neglecting even significant fingerprint fluctuations that may suggest poor complex stability, the IFF coefficient can overcome this issue, further improving the protocol accuracy.

3 Results

To extend the applicability domain of Thermal Titration Molecular Dynamics to the characterization of RNA-ligand complexes, six different test cases were chosen based on the chemical diversity of the ligand and structural diversity of the receptor among the experimental complexes deposited in the Protein Data Bank. Four out of the six test cases were drawn out from the work of Bissaro et al. (Bissaro et al., 2020), specifically 2LWK, 3Q50, 5BJO, and 6E1U, while 1UUD and 1UUI were added because of the challenge provided by the high conformational freedom of the P12 and P14 ligands.

Hereafter a brief overview of the six test cases.

3.1 Influenza A virus RNA promoter

Influenza A virus RNA promoter has a significant role in the modulation of transcription and replication of this group of viruses belonging to the Orthomyxoviridae family. Therefore, it is considered an interesting target for the development of antiviral drugs (Lee et al., 2013). Through a fragment screening approach, 6,7-dimethoxy-2-(1-piperazinyl)-4-quinazolinamine (DPQ) was identified by Varani’s group (Lee et al., 2013). This compound has a micromolar affinity for the promoter region (K_D of 50.5 μM), and consequently, it is a scaffold for further development of antivirals targeting the Influenza A promoter region (PDB ID: 2LWK).

3.2 Corn Aptamer

The Corn aptamer is an in vitro selected aptamer that binds a fluorophore ligand (DFHO), with a nanomolar binding affinity (K_D = 70 nM) (Warner et al., 2017). Thanks to its limited cytotoxicity, the Corn-DFHO could represent a valuable imaging tool. Despite being therapeutically less interesting than the other RNA structures selected for this study, its peculiar three-dimensional organization in a quasi-symmetric dimer with a non-cationic ligand buried at the monomer-monomer interface between two G-quadruplex structures stabilized by K+ ions makes it an appealing target to apply TTMD protocol, since its stability in molecular dynamics, and the low structural flexibility of the fluorophore ligand (3,5-difluoro-4-hydroxybenzylidene imidazolinone-2-oxime) led Bissaro et al. to optimal results when applying SuMD protocol to this system (Bissaro et al., 2020).

3.3 Pre-queuosine 1 riboswitch

PreQ1 (Pre-queosine 1, 7-aminomethyl-7-deazaguanine) is a precursor of the hypermodified guanine nucleotide Queuosine Q), a modified nucleoside that influences anticodon-codon stability in tRNAs, hence increasing translational fidelity. Prokaryotes can synthesize PreQ1 molecules from a multienzyme pathway starting from GTP. By contrast, eukaryotes do not synthesize PreQ1 de novo and need to retrieve queuoine from the diet (Eichhorn et al., 2014). For this reason, PreQ1 riboswitches are only found in prokaryotes, where, by binding PreQ1, they play a fundamental role in the biosynthesis and transport of this molecule. Therefore, riboswitches are an appealing target for antibacterial drug development. Two PreQ1 riboswitch complexes have been selected, the first one (PDB ID: 3Q50) is bound to a PreQ1 precursor with a nanomolar binding affinity (K_D = 2 nM), the second one (PDB ID6E1U) is a hybrid riboswitch aptamer bound to a synthetic compound (HMJ) that presents a sub-micromolar affinity (K_D = 0.5 μM). This complex has been obtained through the deletion of nucleobase A14 and the two vicinal ones, as a direct consequence, the binding site lacks important nucleotides and, congruently to the work of Bissaro et al. (Bissaro et al., 2020), before running TTMD simulations, only the ligand was retained from this complex, and the aptamer structure was retrieved from 3Q50 crystal.

3.4 Tar

Tat-Tar interaction modulates the transcription process in HIV-1. Therefore, targeting Tar RNA could be an appealing pathway to follow for the development of HIV-1 therapeutics. Although combination antiretroviral therapy (cART) can significantly suppress viral load, the eradication of this virus is still challenging. Moreover, chronic long-term comorbidities can arise since viral reservoirs are not eliminated by cART. In this scenario blocking HIV transcription could be beneficial (Alanazi et al., 2021). Several efforts have been made in the characterization of this target and various small-molecule complexes have been reported in the PDB. In this work, two Tar complexes have been reported, specifically selected to assay TTMD behavior with challenging flexible ligands (P12, P14) with high conformational freedom.

The investigated cases are summarized in Table 1 and their crystal structures are reported in Table 2 to showcase the structural differences of these targets.

TABLE 1

TABLE 1. Summary of the six test-cases investigated in the present study.

TABLE 2

TABLE 2. Crystal structures of the six test-cases reported in the article.

For each complex, at first, a self-docking experiment was conducted using the PLANTS docking program, since it has recently been reported as one of the best docking tools for RNA-ligand complexes (Jiang et al., 2023) and it is also free for academics. The results of the self-docking calculations are summarized in Tables 3–8.

TABLE 3

TABLE 3. This panel encompasses all binding modes (experimental + docking poses) investigated in this work for complex deposited in the PDB with accession code 1UUD. For each pose, the ChemPLP docking score and RMSD to the experimental pose are reported.

TABLE 4

TABLE 4. This panel encompasses all binding modes (experimental + docking poses) investigated in this work for complex deposited in the PDB with accession code 1UUI. For each pose, the ChemPLP docking score and RMSD to the experimental pose are reported.

TABLE 5

TABLE 5. This panel encompasses all binding modes (experimental + docking poses) investigated in this work for complex deposited in the PDB with accession code 2LWK. For each pose, the ChemPLP docking score and RMSD to the experimental pose are reported.

TABLE 6

TABLE 6. This panel encompasses all binding modes (experimental + docking poses) investigated in this work for complex deposited in the PDB with accession code 3Q50. For each pose, the ChemPLP docking score and RMSD to the experimental pose are reported.

TABLE 7

TABLE 7. This panel encompasses all binding modes (experimental + docking poses) investigated in this work for complex deposited in the PDB with accession code 5BJO. For each pose, the ChemPLP docking score and RMSD to the experimental pose are reported.

TABLE 8

TABLE 8. This panel encompasses all binding modes (experimental + docking poses) investigated in this work for complex deposited in the PDB with accession code 6E1U. For each pose, the ChemPLP docking score and RMSD to the experimental pose are reported.

For each complex, the top five scoring poses and the experimentally determined one were subjected to TTMD refinement using the same conditions described in the original publications (Pavan et al., 2022c; Menin et al., 2023). Furthermore, another run of TTMD calculations was conducted using an alternative temperature ramp. Two different sets of protein-ligand fingerprints were used for monitoring the binding mode’s evolution and two different metrics for the qualitative estimation of pose stability, namely, the MS and IFF coefficient. Every methodological detail is reported in the Materials and Methods section. Hereafter, each TTMD experiment is discussed separately based on execution conditions, while aggregated considerations are reported in the Discussion.

3.5 TTMD protocol 1: standard temperature ramp, ODDT fingerprints, MS coefficient

As a first attempt, we tried repurposing the same TTMD workflow described in the original publications (Pavan et al., 2022c; Menin et al., 2023) without any modifications. We therefore used a standard temperature ramp, setting up the starting simulation temperature at 300 K, the end temperature at 450 K, the temperature increase between each TTMD-step to 10 K, and the duration of each TTMD-step at 10 ns. Poses were then ranked according to the MS coefficient, using ODDT interaction fingerprints for monitoring the evolution of the binding mode. Results for this first round of TTMD simulations are summarized in Table 9.

TABLE 9

TABLE 9. This table summarizes results for TTMD simulations carried out under the protocol 1 simulation and analysis conditions. For each investigated system, the individual MS coefficient for each TTMD replicate and the average MS value are reported. The closest value to the mean is highlighted in bold.

As can be seen in Table 9, the standard TTMD protocol that was successfully applied to the characterization of protein-ligand complexes’ residence time cannot be readily repurposed to RNA-ligand complexes. What we noticed immediately, indeed, was that several replicates presented remarkably high MS values (1, for example), implying that the ligands lose most or even all the native interaction features in the first TTMD step. By visually inspecting the trajectory and subsequent analyses, we noticed that this was not true, as the ligands were still in place and in a conformation within the binding site that was compatible with a similar interaction pattern as the one of the reference frame. We carefully checked the fingerprint calculation step, and we noticed that ODDT interaction fingerprints cannot capture most of the interactions that exist between the ligand and the RNA receptor, especially concerning hydrogen bonds and pi stacking, which are dominant in this kind of system and captured just a fraction of the hydrophobic contacts. We tried manually modifying the cutoffs for the calculation of such interactions, to see if they were too strict, but the results did not dramatically change. We therefore concluded that it was an atom-typing issue, with ODDT not being able to rightfully classify all RNA atoms and depict the intermolecular interactions involving them. For this reason, we modified the TTMD code to use a different protein-ligand interaction fingerprint, provided by the ProLIF package (Bouysset and Fiorucci, 2021), since it supports DNA and RNA molecules as well.

3.6 TTMD protocol 2: standard temperature ramp, PROLIF fingerprints, MS coefficient

As a second attempt, we re-analyzed existing TTMD trajectories monitoring the evolution of the receptor-ligand interaction fingerprint through the ProLIF package instead of the ODDT one. The result of this analysis is reported in Table 10.

TABLE 10

TABLE 10. This table summarizes results for TTMD simulations carried out under the protocol 2 simulation and analysis conditions. For each investigated system, the individual MS coefficient for each TTMD replicate and the average MS value are reported. The closest value to the mean is highlighted in bold.

As can be observed in Table 10, although the sensitivity of the fingerprint-based scoring function improved compared to the traditional TTMD protocol thanks to the change of the package used for the PLIF calculation, the results of the TTMD simulations are still not satisfactory. Taking aside the first three test cases, i.e., 1UUD, 1UUI, and 2LWK for which the docking itself was not reliable thus justifying the challenging nature of the task, in the other three cases the results were promising but showed a margin of improvement. Indeed, in all these three cases, the TTMD refinement would have at best excluded a couple of wrong poses from the equation but would have attributed some other wrong poses with similar scores to the experimentally determined binding mode, resulting in an inconclusive post-docking refinement procedure. On the other hand, in the case of complex 3Q50, where the closest poses to the crystal reference are pose 1 and 4, respectively, the MS coefficient would have awarded pose 1 and 2, giving them scores similar to the experimentally determined binding pose, and higher scores to the other ones. Intriguingly, the pose ranking capabilities would have been very good in the case of complexes 5BJO and 6E1U, which were the closest poses to the experimental data (1, 4 and 3 for 5BJO, and 3, 1, and 5 for 6E1U, respectively) would have been correctly distinguished by the MS coefficient from the more incorrect ones, other than being close in MS score to the crystal reference and being correctly ranked based on the RMSD to the native binding mode. Based on the visual inspection of TTMD trajectories, and on previous experience accumulated working on similar targets (Bissaro et al., 2020; Pavan et al., 2022a), we thought that a possible solution for improving the accuracy of the prediction could be to use a different temperature ramp, specifically a less aggressive one, to reduce the amplitude of the conformational movements of the RNA receptors to favorably decouple the evolution of the three-dimensional structure of the receptor from the diffusion of the ligand from the binding site.

3.7 TTMD protocol 3: alternative temperature ramp, PROLIF fingerprints, MS coefficient

To further improve the predictive capabilities of the method, we decided to test out an alternative, less aggressive, temperature ramp. Specifically, we set up a starting temperature of 73 K and an end temperature of 223 K, while maintaining the same temperature increase and step duration. The ramp was purposefully chosen within a temperature interval that could be considered “extreme cold” based on a pair of assumptions: first, with highly thermostable macromolecules such as protein receptors investigated in previous TTMD works, working within a temperature interval that is way outside the realistic representation of the system did not impair the ability to extract meaningful data out of the simulation (i.e., a qualitative estimation of the residence time), second lowering the temperature of the simulated system will reduce the atomic velocities, thus resulting in a decreased flexibility of the RNA receptor. Although temperatures in Molecular Dynamics are not directly equivalent to the real-world ones, a similar concept of extracting meaningful data at non-physiological temperatures is commonly applied in Molecular Biology, for example, while determining crystallographic structures with X-ray spectroscopy and cryo-EM techniques. The results of this set of TTMD simulations are encompassed in Table 11.

TABLE 11

TABLE 11. This table summarizes results for TTMD simulations carried out under the protocol 3 simulation and analysis conditions. For each investigated system, the individual MS coefficient for each TTMD replicate and the average MS value are reported. The closest value to the mean is highlighted in bold.

As can be deducted from Table 11, the use of an alternative temperature ramp that operates at lower temperature values did not necessarily improve the quality of the TTMD post-docking refinement. First, the quality of the prediction decreased in those cases (1UUD, 1UUI, and 2LWK) where docking was not able to generate a good docking pose: indeed, while protocol 2 was at least able to assign the lowest scores to a pool of poses that included the experimentally determined ones, here it can be noticed how the native binding mode is always among the ones with the worst scores in term of MS coefficient. Speaking of the other three cases, instead, where protocol 2 performed quite decently, here the pose ranking capabilities of TTMD stay identical, with the only difference being that the MS score difference between various poses is decreased, flattening the difference in score between native-like poses and wrong ones. Looking at the fingerprint profile, though, we noticed how lowering the temperature ranges not only reduced the amplitude of conformational shifts in the nucleic receptor but also decreased the tendency of the ligand to diffuse from the binding site, thus making the protocol less efficient. Practically speaking, although the “titration timeline” profile clearly showed a difference in the amplitude of the changes in the fingerprint-based score throughout the TTMD simulation, the difference between the end and starting point of the simulation was not so marked. This observation led us to conclude that, maybe, the MS coefficient in its original formulation could not be the most sensitive metric to capture the subtle rearrangements of the ligand within the binding site and the conservation of the binding features, thus inducing us to introduce a different metric defined as IFF coefficient (see Materials and Methods for a detailed explanation on how the IFF coefficient is computed).

3.8 TTMD protocol 4: alternative temperature ramp, PROLIF fingerprints, IFF coefficient

To test if an alternative metric to the MS coefficient would be able to improve the sensitivity of the TTMD protocol using an alternative temperature ramp that operates at sub-freezing temperature, we re-analyzed those trajectories determining the IFF coefficient. The results are summarized in Table 12.

TABLE 12

TABLE 12. This table summarizes results for TTMD simulations carried out under the protocol 4 simulation and analysis conditions. For each investigated system, the individual IFF coefficient for each TTMD replicate and the average IFF value are reported. The closest value to the mean is highlighted in bold.

By observing Table 12, it can be seen that even the new metric is not able to improve the results of the TTMD refinement performed at much lower temperatures. Specifically, the performance of the method remains like what could be observed with the old-fashioned MS coefficient. To further prove the point, we re-analyzed those trajectories performed with the standard temperature ramp (see Supplementary Table S1) and we found a great level of coherence between the MS and IFF coefficient regarding their ability to rank poses based on the persistence of native binding features throughout the simulation.

4 Discussion

Molecular docking represents the state-of-the-art technique for the prediction of protein-ligand complexes, thanks to the good compromise between accuracy and rapidity of execution (Pavan and Moro, 2023). Although routinely used in various structure-based drug discovery campaigns, molecular docking has been specifically designed and optimized with protein-ligand receptor in mind, due to pharmaceutical relevance and availability of experimentally determined structures to use as a benchmark for method development (Salmaso and Moro, 2018). Despite the growing interest of medicinal chemists in targeting RNA macromolecules for therapeutic and diagnostic purposes in the last decade, the field of rational design of RNA-targeting ligands is still relatively unexplored, due to some intrinsic challenges portrayed by these molecular entities, such as their structural plasticity which makes it difficult to obtain their experimentally determined three-dimensional structure, a pivotal requirement for modern-days drug discovery campaigns (Bissaro et al., 2020; Pavan et al., 2022a).

Despite all these challenges, recent works have demonstrated how, despite some intrinsic limitations, routinely used molecular docking programs such as PLANTS can be quite efficiently used for investigating RNA-ligand complexes as well (Jiang et al., 2023). Furthermore, we previously showcased how Supervised Molecular Dynamics (SuMD) simulations can be successfully utilized to flank molecular docking in investigating the recognition process of complexes involving RNA molecules, both as ligands (Pavan et al., 2022a) and as receptors (Bissaro et al., 2020), retrieving useful information about the whole binding process beyond the final bound state. Critically, one of the major problems of docking in all its iterations, both classic and dynamic, is finding a good scoring metric that can discriminate native-like poses from the wrong ones (Chaput and Mouawad, 2017). The peculiarities of RNA molecules, such as the distinctive surface charge distribution, exacerbate the limitations of classical scoring functions, which have been developed specifically for protein binding sites. For all these reasons, introducing new and improved ways of investigating RNA-ligand complexes is a very pressing issue.

In the present article, we presented the first application of Thermal Titration Molecular Dynamics (TTMD), a recently developed MD-based protocol for the qualitative estimation of protein-ligand unbinding kinetics, to the characterization of RNA-ligand complexes. We investigated six different pharmaceutically relevant test cases of different molecular complexity, both on the ligand and on the receptor side. We performed two rounds of simulations, one with the standard temperature ramp already described in the original publications (Pavan et al., 2022c; Menin et al., 2023), and one with an alternative ramp operating at sub-freezing temperatures, to investigate the effect of the temperature ramp on the accuracy of prediction. We modified the original code to use the ProLIF package (Bouysset and Fiorucci, 2021) for the calculation of RNA-ligand interaction fingerprint instead of the ODDT one (Wójcikowski et al., 2015), since they exhibit a better description of these systems. We also explored a different metric, the IFF coefficient, to use alongside the MS coefficient for monitoring the conservation of the native binding determinants and calculating a proxy measure for the receptor-ligand residence time.

By analyzing the whole set of simulations and analysis performed upon them, it seems clear that TTMD can be a helpful tool to refine docking results in those cases where the geometry of the binding site is well defined and maintained throughout the simulation, in a similar way to what was already observed for protein-ligand systems. 5BJO simulations are proof of this concept as the binding site remains extremely stable throughout the simulation (Video V1, Supplementary Table S1, Supplementary Table S2). In those cases where the perturbation of the receptor conformation and or the receptor fold overcomes the tendency of the ligand to diffuse from the binding site, the method accuracy drops significantly, not improving docking results and even invalidating successful ones.

Concerning this aspect, the choice of the temperature ramp is pivotal for the success of the TTMD rescoring process. Despite the high flexibility of RNA molecules, especially at high simulation temperatures, the standard temperature ramp still outperformed the alternative, sub-freezing one. A plausible explanation for this evidence is that the electrostatic component dominates the protein-ligand interaction energy, making it very unlikely for the ligand to spontaneously detach and unbind without any external aid in the shape of a temperature increase. Regardless of the choice of the metric used for estimating the pose residence time (the MS or the IFF coefficient), the ranking provided by the standard ramp was more useful than the one provided by the alternative one, due to the flattening of the score differences between good and bad poses at lower temperatures. Despite this, the exploration of the alternative ramp allowed us to determine that simulating in those conditions will not result in discarding a good pose, at least for those cases where TTMD (and docking too) perform well. The possibility to tune the temperature interval with a certain degree of freedom based on the knowledge of the system makes it quite appealing, since in a prospective use case one could simply start with a titration in the low-temperature range, evaluate the separation in score between different poses, and eventually re-execute the titration at higher temperatures based on the results provided by the first round of simulations.

Another useful expediency that the user could implement would be to carry out a TTMD “dry run” on just the RNA molecule itself, a functionality that has been introduced purposefully in the latest version of the TTMD code used for this article, and evaluate its behavior in terms of structural plasticity, especially at the binding site level, before running the titration on systems including the ligand as well. It is important to stress that probably exists an intermediate ramp between the two that were utilized in the present work that might be a better compromise and provide slightly better results, but at the same time, we firmly believe that to be adopted by the largest number of user possible a method has to be easy to implement and setup, without spending too much time on optimizing parameters for its execution. For this reason, we wanted to showcase that the method can perform reasonably well almost regardless of the choice of the ramp, although the original ramp that has already proven to be successful in the case of protein-ligand complexes is still the best one in terms of accuracy. Finally, it is worth mentioning that working with macromolecules that have an intrinsic conformational flexibility such as RNA strongly limits the applicability of any given MD technique in a high-throughput fashion, due to its limited sampling capabilities, making docking still the most efficient choice for investigating these systems. On the other hand, the ever-increasing computational power available for MD simulations will make and more appealing the use of pipelines like TTMD in the future, thus justifying the interest in the further development of this and other MD-based methods for the investigation of systems involving nucleic acids (Sponer et al., 2018).

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

AD: Formal Analysis, Investigation, Writing–original draft. MP: Investigation, Methodology, Writing–original draft. SMe: Data curation, Formal Analysis, Software, Visualization, Writing–original draft. VS: Methodology, Supervision, Validation, Writing–review and editing. MS: Supervision, Supervision, Writing–review and editing, Formal Analysis. SMo: Supervision, Conceptualization, Writing–review and editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work has been funded by the European Union - Next-Generation EU (“PNRR M4C2-Investimento 1.4- CN00000041”). MMS lab is very grateful to Chemical Computing Group, OpenEye, and Acellera for their scientific and technical partnership.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2023.1294543/full#supplementary-material

References

Alanazi, A., Ivanov, A., Kumari, N., Lin, X., Wang, S., Kovalskyy, D., et al. (2021). Targeting tat–tar rna interaction for hiv-1 inhibition. Viruses 13, 2004. doi:10.3390/v13102004

PubMed Abstract | CrossRef Full Text | Google Scholar

Bassani, D., and Moro, S. (2023). Past, present, and future perspectives on computer-aided drug design methodologies. Molecules 28, 3906. doi:10.3390/MOLECULES28093906

PubMed Abstract | CrossRef Full Text | Google Scholar

Berman, H. M., Westbrook, J., Feng, Z., Gilliland, G., Bhat, T. N., Weissig, H., et al. (2000). The protein Data Bank. Nucleic Acids Res. 28, 235–242. doi:10.1093/NAR/28.1.235

PubMed Abstract | CrossRef Full Text | Google Scholar

Bernardi, R. C., Melo, M. C. R., and Schulten, K. (2015). Enhanced sampling techniques in molecular dynamics simulations of biological systems. Biochimica Biophysica Acta (BBA) - General Subj. 1850, 872–877. doi:10.1016/J.BBAGEN.2014.10.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacharyya, D., Mirihana Arachchilage, G., and Basu, S. (2016). Metal cations in G-quadruplex folding and stability. Front. Chem. 4, 38. doi:10.3389/fchem.2016.00038

PubMed Abstract | CrossRef Full Text | Google Scholar

Bissaro, M., Sturlese, M., and Moro, S. (2020). Exploring the RNA-recognition mechanism using supervised molecular dynamics (SuMD) simulations: toward a rational design for ribonucleic-targeting molecules? Front. Chem. 8, 107. doi:10.3389/fchem.2020.00107

PubMed Abstract | CrossRef Full Text | Google Scholar

Bouysset, C., and Fiorucci, S. (2021). ProLIF: a library to encode molecular interactions as fingerprints. J. Cheminformatics 13, 72–79. doi:10.1186/s13321-021-00548-6

CrossRef Full Text | Google Scholar

Cadence Molecular Sciences (2023). OpenEye. Cadence Molecular Sciences, Santa Fe, NM. http://www.eyesopen.com. Available at:http://www.eyesopen.com.

Google Scholar

Case, D. A. (2022). AmberTools22.

Google Scholar

Case, D. A., Cheatham, T. E., Darden, T., Gohlke, H., Luo, R., Merz, K. M., et al. (2005). The Amber biomolecular simulation programs. J. Comput. Chem. 26, 1668–1688. doi:10.1002/JCC.20290

PubMed Abstract | CrossRef Full Text | Google Scholar

Chaput, L., and Mouawad, L. (2017). Efficient conformational sampling and weak scoring in docking programs? Strategy of the wisdom of crowds. J. Cheminformatics 9, 37. doi:10.1186/s13321-017-0227-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Childs-Disney, J. L., Yang, X., Gibaut, Q. M. R., Tong, Y., Batey, R. T., and Disney, M. D. (2022). Targeting RNA structures with small molecules. Nat. Rev. Drug Discov., 736–762. 2022 21:10 21. doi:10.1038/s41573-022-00521-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Connelly, C. M., Numata, T., Boer, R. E., Moon, M. H., Sinniah, R. S., Barchi, J. J., et al. (2019). Synthetic ligands for PreQ1 riboswitches provide structural and mechanistic insights into targeting RNA tertiary structure. Nat. Commun., 1501–1512. 2019 10:1 10. doi:10.1038/s41467-019-09493-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Copeland, R. A., Pompliano, D. L., and Meek, T. D. (2006). Drug–target residence time and its implications for lead optimization. Nat. Rev. Drug Discov. 5, 730–739. doi:10.1038/nrd2082

PubMed Abstract | CrossRef Full Text | Google Scholar

Davidchack, R. L., Handel, R., and Tretyakov, M. V. (2009). Langevin thermostat for rigid body dynamics. J. Chem. Phys. 130, 234101. doi:10.1063/1.3149788

PubMed Abstract | CrossRef Full Text | Google Scholar

Davis, B., Afshar, M., Varani, G., Murchie, A. I. H., Karn, J., Lentzen, G., et al. (2004). Rational design of inhibitors of HIV-1 TAR RNA through the stabilisation of electrostatic “hot spots.”. J. Mol. Biol. 336, 343–356. doi:10.1016/J.JMB.2003.12.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Do, T. N., Carloni, P., Varani, G., and Bussi, G. (2013). RNA/peptide binding driven by electrostatics - insight from bidirectional pulling simulations. J. Chem. Theory Comput. 9, 1720–1730. doi:10.1021/ct3009914

PubMed Abstract | CrossRef Full Text | Google Scholar

Eastman, P., Swails, J., Chodera, J. D., McGibbon, R. T., Zhao, Y., Beauchamp, K. A., et al. (2017). OpenMM 7: rapid development of high performance algorithms for molecular dynamics. PLOS Comput. Biol. 13, e1005659. doi:10.1371/journal.pcbi.1005659

PubMed Abstract | CrossRef Full Text | Google Scholar

Eichhorn, C. D., Kang, M., and Feigon, J. (2014). Structure and function of preQ1 riboswitches. Biochimica Biophysica Acta (BBA) - Gene Regul. Mech. 1839, 939–950. doi:10.1016/J.BBAGRM.2014.04.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Essmann, U., Perera, L., Berkowitz, M. L., Darden, T., Lee, H., and Pedersen, L. G. (1995). A smooth particle mesh Ewald method. J. Chem. Phys. 103, 8577–8593. doi:10.1063/1.470117

CrossRef Full Text | Google Scholar

Falese, J. P., Donlic, A., and Hargrove, A. E. (2021). Targeting RNA with small molecules: from fundamental principles towards the clinic. Chem. Soc. Rev. 50, 2224–2243. doi:10.1039/D0CS01261K

PubMed Abstract | CrossRef Full Text | Google Scholar

Faller, R., and de Pablo, J. J. (2002). Constant pressure hybrid Molecular Dynamics–Monte Carlo simulations. J. Chem. Phys. 116, 55–59. doi:10.1063/1.1420460

CrossRef Full Text | Google Scholar

Ferreira, L. G., Santos, R. N. D., Oliva, G., and Andricopulo, A. D. (2015). Molecular docking and structure-based drug design strategies. Molecules 20, 13384–13421. doi:10.3390/MOLECULES200713384

PubMed Abstract | CrossRef Full Text | Google Scholar

Gowers, R. J., Linke, M., Barnoud, J., Reddy, T. J. E., Melo, M. N., Seyler, S. L., et al. (2016). MDAnalysis: a Python package for the rapid analysis of molecular dynamics simulations. Proc. 15th Python Sci. Conf., 98–105. doi:10.25080/Majora-629e541a-00e

CrossRef Full Text | Google Scholar

Harvey, M. J., Giupponi, G., and Fabritiis, G. D. (2009). ACEMD: accelerating biomolecular dynamics in the microsecond time scale. J. Chem. Theory Comput. 5, 1632–1639. doi:10.1021/CT9000685/ASSET/IMAGES/LARGE/CT-2009-000685_0004.JPEG

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphrey, W., Dalke, A., and Schulten, K. (1996). VMD: visual molecular dynamics. J. Mol. Graph. 14, 33–38. doi:10.1016/0263-7855(96)00018-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Izrailev, S., Stepaniants, S., Isralewitz, B., Kosztin, D., Lu, H., Molnar, F., et al. (1999). “Steered molecular dynamics,” in Computational molecular dynamics: challenges, methods, ideas lecture notes in computational science and engineering. Editors P. Deuflhard, J. Hermans, B. Leimkuhler, A. E. Mark, S. Reich, and R. D. Skeel (Berlin, Heidelberg: Springer), 39–65. doi:10.1007/978-3-642-58360-5_2

CrossRef Full Text | Google Scholar

Jenkins, J. L., Krucinska, J., McCarty, R. M., Bandarian, V., and Wedekind, J. E. (2011). Comparison of a PreQ1 riboswitch aptamer in metabolite-bound and free states with implications for gene regulation. J. Biol. Chem. 286, 24626–24637. doi:10.1074/JBC.M111.230375

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, D., Zhao, H., Du, H., Deng, Y., Wu, Z., Wang, J., et al. (2023). How good are current docking programs at nucleic acid-ligand docking? A comprehensive evaluation. J. Chem. Theory Comput. 19, 5633–5647. doi:10.1021/ACS.JCTC.3C00507/ASSET/IMAGES/LARGE/CT3C00507_0006.JPEG

PubMed Abstract | CrossRef Full Text | Google Scholar

Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983). Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. doi:10.1063/1.445869

CrossRef Full Text | Google Scholar

Korb, O., Stützle, T., and Exner, T. E. (2006). “PLANTS: application of ant colony optimization to structure-based drug design,” in Ant colony Optimization and swarm intelligence lecture notes in computer science. Editors M. Dorigo, L. M. Gambardella, M. Birattari, A. Martinoli, R. Poli, and T. Stützle (Berlin, Heidelberg: Springer), 247–258. doi:10.1007/11839088_22

CrossRef Full Text | Google Scholar

Korb, O., Stützle, T., and Exner, T. E. (2007). An ant colony optimization approach to flexible protein–ligand docking. Swarm Intell. 1, 115–134. doi:10.1007/s11721-007-0006-9

CrossRef Full Text | Google Scholar

Korb, O., Stützle, T., and Exner, T. E. (2009). Empirical scoring functions for advanced Protein−Ligand docking with PLANTS. J. Chem. Inf. Model. 49, 84–96. doi:10.1021/ci800298z

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, M. K., Bottini, A., Kim, M., Bardaro, M. F., Zhang, Z., Pellecchia, M., et al. (2013). A novel small-molecule binds to the influenza A virus RNA promoter and inhibits viral replication. Chem. Commun. 50, 368–370. doi:10.1039/C3CC46973E

PubMed Abstract | CrossRef Full Text | Google Scholar

Lüdemann, S. K., Lounnas, V., and Wade, R. C. (2000). How do substrates enter and products exit the buried active site of cytochrome P450cam? 1. Random expulsion molecular dynamics investigation of ligand access channels and mechanisms11Edited by J. Thornton. J. Mol. Biol. 303, 797–811. doi:10.1006/jmbi.2000.4154

PubMed Abstract | CrossRef Full Text | Google Scholar

Maier, J. A., Martinez, C., Kasavajhala, K., Wickstrom, L., Hauser, K. E., and Simmerling, C. (2015). ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 11, 3696–3713. doi:10.1021/ACS.JCTC.5B00255/SUPPL_FILE/CT5B00255_SI_001.PDF

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, X.-Y., Zhang, H.-X., Mezei, M., and Cui, M. (2012). Molecular docking: a powerful approach for structure-based drug discovery. Curr. Comput. Aided-Drug Des. 7, 146–157. doi:10.2174/157340911795677602

PubMed Abstract | CrossRef Full Text | Google Scholar

Menin, S., Pavan, M., Salmaso, V., Sturlese, M., and Moro, S. (2023). Thermal titration molecular dynamics (TTMD): not your usual post-docking refinement. Int. J. Mol. Sci. 24, 3596. doi:10.3390/ijms24043596

PubMed Abstract | CrossRef Full Text | Google Scholar

Michaud-Agrawal, N., Denning, E. J., Woolf, T. B., and Beckstein, O. (2011). MDAnalysis: a toolkit for the analysis of molecular dynamics simulations. J. Comput. Chem. 32, 2319–2327. doi:10.1002/jcc.21787

PubMed Abstract | CrossRef Full Text | Google Scholar

Mlýnský, V., and Bussi, G. (2018). Exploring RNA structure and dynamics through enhanced sampling simulations. Curr. Opin. Struct. Biol. 49, 63–71. doi:10.1016/j.sbi.2018.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Molecular Operating Environment (MOE), 2022.02, C. C. G. (2022). Molecular operating environment (MOE), 2022.02. Available at: https://www.chemcomp.com/Research-Citing_MOE.htm.

Google Scholar

O’Boyle, N. M., Banck, M., James, C. A., Morley, C., Vandermeersch, T., and Hutchison, G. R. (2011). Open Babel: an open chemical toolbox. J. Cheminformatics 3, 33–14. doi:10.1186/1758-2946-3-33

CrossRef Full Text | Google Scholar

Pavan, M., Bassani, D., Sturlese, M., and Moro, S. (2022a). Investigating RNA–protein recognition mechanisms through supervised molecular dynamics (SuMD) simulations. NAR Genomics Bioinforma. 4, lqac088. doi:10.1093/NARGAB/LQAC088

CrossRef Full Text | Google Scholar

Pavan, M., Menin, S., Bassani, D., Sturlese, M., and Moro, S. (2022b). Implementing a scoring function based on interaction fingerprint for Autogrow4: protein kinase CK1δ as a case study. Front. Mol. Biosci. 9, 909499. doi:10.3389/fmolb.2022.909499

PubMed Abstract | CrossRef Full Text | Google Scholar

Pavan, M., Menin, S., Bassani, D., Sturlese, M., and Moro, S. (2022c). Qualitative estimation of protein–ligand complex stability through thermal titration molecular dynamics simulations. J. Chem. Inf. Model. 62, 5715–5728. doi:10.1021/acs.jcim.2c00995

PubMed Abstract | CrossRef Full Text | Google Scholar

Pavan, M., and Moro, S. (2023). Lessons learnt from COVID-19: computational strategies for facing present and future pandemics. Int. J. Mol. Sci. 24, 4401. doi:10.3390/IJMS24054401

PubMed Abstract | CrossRef Full Text | Google Scholar

Pérez, A., Marchán, I., Svozil, D., Sponer, J., Cheatham, T. E., Laughton, C. A., et al. (2007). Refinement of the AMBER force field for nucleic acids: improving the description of alpha/gamma conformers. Biophysical J. 92, 3817–3829. doi:10.1529/BIOPHYSJ.106.097782

CrossRef Full Text | Google Scholar

Salmaso, V., and Moro, S. (2018). Bridging molecular docking to molecular dynamics in exploring ligand-protein recognition process: an overview. Front. Pharmacol. 9, 393738. doi:10.3389/fphar.2018.00923

CrossRef Full Text | Google Scholar

So/rensen, M. R., and Voter, A. F. (2000). Temperature-accelerated dynamics for simulation of infrequent events. J. Chem. Phys. 112, 9599–9606. doi:10.1063/1.481576

CrossRef Full Text | Google Scholar

Sponer, J., Bussi, G., Krepl, M., Banas, P., Bottaro, S., Cunha, R. A., et al. (2018). RNA structural dynamics as captured by molecular simulations: a comprehensive overview. Chem. Rev. 118, 4177–4338. doi:10.1021/ACS.CHEMREV.7B00427

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugita, Y., and Okamoto, Y. (1999). Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 314, 141–151. doi:10.1016/S0009-2614(99)01123-9

CrossRef Full Text | Google Scholar

Torrie, G. M., and Valleau, J. P. (1977). Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 23, 187–199. doi:10.1016/0021-9991(77)90121-8

CrossRef Full Text | Google Scholar

Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004). Development and testing of a general amber force field. J. Comput. Chem. 25, 1157–1174. doi:10.1002/JCC.20035

PubMed Abstract | CrossRef Full Text | Google Scholar

Warner, K. D., Sjekloa, L., Song, W., Filonov, G. S., Jaffrey, S. R., and Ferré-D’Amaré, A. R. (2017). A homodimer interface without base pairs in an RNA mimic of red fluorescent protein. Nat. Chem. Biol., 1195–1201. 2017 13:11 13. doi:10.1038/nchembio.2475

PubMed Abstract | CrossRef Full Text | Google Scholar

Wójcikowski, M., Zielenkiewicz, P., and Siedlecki, P. (2015). Open Drug Discovery Toolkit (ODDT): a new open-source player in the drug discovery field. J. Cheminformatics 7, 26. doi:10.1186/S13321-015-0078-2

CrossRef Full Text | Google Scholar

Zafferani, M., and Hargrove, A. E. (2021). Small molecule targeting of biologically relevant RNA tertiary and quaternary structures. Cell Chem. Biol. 28, 594–609. doi:10.1016/J.CHEMBIOL.2021.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Zamora, R. J., Uberuaga, B. P., Perez, D., and Voter, A. F. (2016). The modern temperature-accelerated dynamics approach. Annu. Rev. Chem. Biomol. Eng. 7, 87–110. doi:10.1146/annurev-chembioeng-080615-033608

PubMed Abstract | CrossRef Full Text | Google Scholar

Zgarbová, M., Otyepka, M., Šponer, J., Mládek, A., Banáš, P., Cheatham, T. E., et al. (2011). Refinement of the Cornell et al. Nucleic acids force field based on reference quantum chemical calculations of glycosidic torsion profiles. J. Chem. Theory Comput. 7, 2886–2902. doi:10.1021/ct200162x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Jiang, Y., and Chen, S.-J. (2022). RNA–ligand molecular docking: advances and challenges. WIREs Comput. Mol. Sci. 12, e1571. doi:10.1002/wcms.1571

CrossRef Full Text | Google Scholar