Mechanism of Ligand Binding to Theophylline RNA Aptamer

Surface Plasmon Resonance (SPR).

SPR experiments were performed on a Biacore T200 (GE Healthcare) instrument at 25 °C using streptavidin-precoated SA sensor chips (GE Healthcare). The running buffer, which contained 20 mM HEPES, 100 mM NaCl, and 5 mM MgCl₂ at pH 7.3, was freshly made and used after being filtered through a 0.22 μm PVDF membrane. GenScript, Inc. synthesized the 3′-biotinylated TEP aptamer RNA (i.e., GGCGAUACCAGCCGAAAGGCCCUUGGCAGCGUCUU/3′Biotin/), which we reconstituted at 0.1 mM in molecular biology water. For immobilization of the biotinylated RNA, the sensor chip was first conditioned with 3 consecutive 1 min injections of a high salt solution (50 mM NaOH, 1 M NaCl) at a flow rate of 50 μL/min. Next, the biotinylated RNA was diluted in the running buffer (1 μM) and applied over the streptavidin sensor chip surface at a flow rate of 10 μL/min to achieve an immobilization level of about 1000 RU. Finally, alkyne-PEG-biotin (50 μM in running buffer) was injected (1 min, 10 μL/min) to block the remaining streptavidin surface binding sites. The BiaControl Software Wizard Kinetics methodology was used for the kinetic analysis. In order to titrate over the immobilized RNA (contact time: 1 min, flow rate: 50 μL/min), the TEP or CFF was dissolved in the running buffer to the necessary concentrations (0.032, 0.064, 0.32, 0.64, and 3.2 μM for the TEP, and 0.35, 0.7, 3.5, 7, and 35 mM for the CFF).

The data analysis and plotting were performed by using BiaEvaluation Software. All monitored resonance signals were subtracted with signals from a nonbinding reference channel. Equilibrium and rate constants ($K_{\mathrm{d}},k_{\mathrm{o}\mathrm{n}}$, and $k_{\text{off}}$) were calculated using the BiaEvaluation Software Binding Affinity protocol with 1:1 fitting.

Ligand Gaussian Accelerated Molecular Dynamics (LiGaMD).

Ligand Gaussian accelerated Molecular Dynamics (LiGaMD) is a computational method developed to enhance the sampling of ligand binding and dissociation in target receptors. Here, we briefly outline the algorithm of LiGaMD.⁶

We consider a scenario in which a ligand $\mathrm{L}$ binds to a receptor $\mathrm{R}$ in biological environment $\mathrm{E}.\mathrm{N}$ atoms with their locations make up the system $r\equiv \left\{{\overrightarrow{r}}_1,\cdots ,\overrightarrow{r_N}\right\}$ and momenta $p\equiv \left\{\overrightarrow{p_1},\cdots ,\overrightarrow{p_N}\right\}$. The Hamiltonian of the system can be expressed as

where $K\left(p\right)$ and $V\left(r\right)$ are the system kinetic and total potential energies, respectively. Next, the potential energy is broken down into the following terms:

where $V_{\mathrm{R},\mathrm{b}},V_{\mathrm{L},\mathrm{b}}$, and $V_{\mathrm{E},\mathrm{b}}$ are the bonded potential energies in the receptor $\mathrm{R}$, ligand $\mathrm{L}$ and environment $\mathrm{E}$, respectively. The self-nonbonded potential energies in receptor $\mathrm{R}$, Ligand $\mathrm{L}$, and Environment $\mathrm{E}$, respectively, are $V_{\mathrm{R}\mathrm{R},\mathrm{n}\mathrm{b}},V_{\mathrm{L}\mathrm{L},\mathrm{n}\mathrm{b}}$ and $V_{\mathrm{E}\mathrm{E},\mathrm{n}\mathrm{b}}.\mathrm{R}-\mathrm{L},\mathrm{R}-\mathrm{E}$, and $\mathrm{L}-\mathrm{E}$’s related nonbonded interaction energies are $V_{\mathrm{R}\mathrm{L},\mathrm{n}\mathrm{b}},V_{\mathrm{R}\mathrm{E},\mathrm{n}\mathrm{b},}$, and $V_{\mathrm{L}\mathrm{E},\mathrm{n}\mathrm{b},}$ respectively.⁶ According to classical molecular mechanics force fields,¹⁸ the nonbonded potential energies are calculated as

where $V_{\text{elec}}$ and $V_{\text{vdW}}$ are the system electrostatic and van der Waals potential energies.⁶ Ligand binding mainly involves the nonbonded interaction energies of the ligand, $V_{\mathrm{L},\mathrm{n}\mathrm{b}}(r)=V_{\mathrm{L}\mathrm{L},\mathrm{n}\mathrm{b}}\left(r_{\mathrm{L}}\right)+V_{\mathrm{R}\mathrm{L},\mathrm{n}\mathrm{b}}\left(r_{\mathrm{P}\mathrm{L}}\right)+V_{\mathrm{L}\mathrm{E},\mathrm{n}\mathrm{b}}\left(r_{\mathrm{L}\mathrm{E}}\right)$. We selectively provide a boost potential to the ligand nonbonded potential energy in LiGaMD as

where $E_{\mathrm{L},\mathrm{n}\mathrm{b}}$ is the threshold energy for applying boost potential and $k_{\mathrm{L},\mathrm{n}\mathrm{b}}$ is the harmonic constant. The subscript of $\mathrm{\Delta }{V}_{\mathrm{L},\mathrm{n}\mathrm{b}}(r)$, $E_{\mathrm{L},\mathrm{n}\mathrm{b}}$, and $k_{\mathrm{L},\mathrm{n}\mathrm{b}}$ is dropped in the section that follows for the sake of simplicity. The simulation parameters $E$ and $k$ are determined according to the GaMD enhanced sampling principles. ⁶

When the system’s maximum potential energy $E$ is set to the lower bound $\left(E=V_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)$, the effective harmonic force constant $k_0$ can be calculated as

where $V_{\mathrm{m}\mathrm{a}\mathrm{x}},V_{\mathrm{m}\mathrm{i}\mathrm{n}},V_{\text{avg}}$ and ${\sigma }_{\mathrm{V}}$ are the maximum, minimum, average, and standard deviation of the boosted system potential energy, and ${\sigma }_0$ is the user-specified upper limit of the standard deviation of $\mathrm{\Delta }V$ (e.g., $10k_{\mathrm{B}}\mathrm{T}$) for proper reweighting.⁶ The upper limit of the standard deviation of the first potential boost and second potential boost that allows for accurate reweighting in dual-boost simulations is referred to as ${\sigma }_{0\mathrm{P}}$ and ${\sigma }_{0\mathrm{D}}$, respectively. ⁶ The harmonic constant is calculated as $k=k_0\cdot \frac{1}{V_{\mathrm{m}\mathrm{a}\mathrm{x}}-V_{\mathrm{m}\mathrm{i}\mathrm{n}}}$ with $0<k_0\le 1$. Alternatively, when the threshold energy $E$ is set to its upper bound $E=V_{\text{min}}+\frac{1}{k},k_0$ is set to

if $k_0^{\prime \prime }$ is found to be between 0 and 1. Otherwise, $k_0$ is calculated using eq 5.

In order to facilitate ligand binding to RNAs in MD simulations, additional ligand molecules might be added to the solvent.¹⁹ This is based on the fact that the time needed for ligand binding is inversely proportional to the ligand concentration.²⁰ The higher the ligand concentration, the faster the ligand binds, provided that the ligand concentration is still within its solubility limits.^6,20 In addition to selectively boosting the bound ligand to accelerate its dissociation, another boost potential is applied to the unbound ligand molecules, RNA, and solvent to facilitate ligand rebinding.⁶ The second boost potential is calculated using the total system potential energy other than the nonbonded potential energy of the bound ligand as

where $E_{\mathrm{D}}$ and $k_{\mathrm{D}}$ are the corresponding threshold energy for applying the second boost potential and the harmonic constant, respectively. This leads to dual-boost LiGaMD with the total boost potential $\mathrm{\Delta }V(r)=\mathrm{\Delta }{V}_{\mathrm{L},\mathrm{n}\mathrm{b}}(r)+\mathrm{\Delta }{V}_{\mathrm{D}}(r)$.⁶

Ligand Binding Free Energy Calculations from 3D Potential of Mean Force.

From the 3D potential of mean force (PMF) of ligand displacements from the target receptor, we calculate the ligand binding free energy as²¹

where $V_0$ is the standard volume, $V_{\mathrm{b}}={\int }_b{e}^{-\beta W(r)}\mathrm{d}r$ is the average sampled bound volume of the ligand with $\beta =1/k_{\mathrm{B}}T,k_{\mathrm{B}}$ is the Boltzmann constant, $T$ is the temperature, and $\mathrm{\Delta }{W}_{3\mathrm{D}}$ is the depth of the 3D PMF. $\mathrm{\Delta }{W}_{3\mathrm{D}}$ can be calculated by integrating Boltzmann distribution of the 3D PMF $W(r)$ overall system coordinates except the $x,y,z$ of the ligand:

where $V_{\mathrm{u}}={\int }_u\mathrm{d}r$ is the sampled unbound volume of the ligand. The exact definitions of the bound and unbound volumes $V_{\mathrm{b}}$ and $V_{\mathrm{u}}$ are not important as the exponential average cut-off contributions are far away from the PMF minima.²¹ The PyReweighting tool kit (http://miao.compbio.ku.edu/PyReweighting/) includes a free Python script called PyReweighting-3D.py that can be used to compute the 3D PMF and related ligand binding free energies. It functions for both enhanced sampling simulations using LiGaMD with energetic reweighting and cMD (without energetic reweighting).⁶

Ligand Binding Kinetics Obtained from Reweighting of LiGaMD Simulations.

One can record the time periods and compute their averages for the ligand found in the bound $\left({\tau }_{\mathrm{B}}\right)$ and unbound $\left({\tau }_{\mathrm{U}}\right)$ states from the simulation trajectories, assuming adequate sampling of repetitive ligand dissociation and binding.²⁰ The ${\tau }_{\mathrm{B}}$ corresponds to residence time in drug design.^6,22 The binding and dissociation rate constants ($k_{\text{off}}$ and $k_{\text{on}}$) for ligands were calculated as

where $\left[\mathrm{L}\right]$ is the ligand concentration in the simulation system.

The ligand kinetics from the LiGaMD simulations are reweighted using Kramers’ rate theory. According to Kramers’ rate theory, the rate of a chemical reaction in the large viscosity limit is determined by²³

where $w_{\mathrm{m}}$ and $w_{\mathrm{b}}$ are frequencies of the approximated harmonic oscillators (also referred to as curvatures of free energy surface^24,25 near the energy minimum and barrier, respectively, $\xi$ is the frictional rate constant, and $\mathrm{\Delta }F$ is the free energy barrier of transition. The friction constant $\mathrm{\xi }$ is related to the diffusion coefficient $D$ with $\xi =k_{\mathrm{B}}T/D$. The apparent diffusion coefficient $D$ is calculated by dividing the kinetic rate calculated using the transition time series directly obtained from simulations by the probability density solution of the Smoluchowski equation.²⁶ The free energy barriers of ligand binding and dissociation are calculated from the original LiGaMD simulations in order to reweight ligand kinetics using Kramer’s rate theory (reweighted, $\mathbf{\Delta }\mathit{F}$) and modified (no reweighting, $\mathbf{\Delta }\mathit{F}*$) PMF profiles, similarly for curvatures of the reweighed ($w$) and modified ($w*$, no reweighting) PMF profiles near the guest bound (“B”) and unbound (“U”) low-energy wells and the energy barrier (“Br”), and the ratio of apparent diffusion coefficients from LiGaMD simulations without reweighting (modified, $D*$) and with reweighting ($D$).²⁰ The resulting numbers are then plugged into eq 9 to estimate accelerations of the ligand binding and dissociation rates during LiGaMD simulations, which allows us to recover the original kinetic rate constants.²³

System Setup.

The NMR structure of the TEP-bound RNA complex (PDB ID: 1O15) was obtained for simulation.^27,28 To create the RNA-CFF complex, a methyl group was added to the N7 atom of TEP in the NMR structure, while keeping the coordinates of the rest of the receptor. The system was energy minimized to remove possible steric clashes. The RNA.OL3^29,30 force field was employed for RNA and GAFF2³¹ for small molecules with the tleap module in the AMBER20 package.³² To replicate the conditions of the solution structure in a TIP3P³³ water box using tleap, both systems were neutralized with 0.01 M MgCl₂³⁴ and 0.15 M NaCl.^32,35 A total of 3 ligands (one in the NMR bound conformation and 2 others placed randomly in the bulk solvent outside at a distance of >15 Å from the RNA surface) were included. The design is based on the fact that the time needed for ligand binding is inversely proportional to the ligand concentration.²⁰ The greater the ligand concentration, the less time is required for ligand binding, provided that the ligand concentration is within solubility limits.²⁰

LiGaMD Simulation Protocol.

The GPU-accelerating program pmemd.cuda in AMBER20 was used to perform LiGaMD simulations.³² The simulation system was minimized by using the steepest descent algorithm. The system was then heated from 0 to 300 K for 200 ps. It was further equilibrated using the NVT ensemble at 300 K for 800 ps and the NPT ensemble at 300 K and 1 bar for 1 ns with 1 kcal/mol/Ų constraints on the heavy atoms of the receptor and ligand, followed by a 2 ns short cMD without any constraint. The LiGaMD simulations proceeded with a 10 ns short cMD to collect the potential statistics, a 63 ns GaMD equilibration after adding the boost potential, and then three independent 2500 ns production runs for TEP. After a 63 ns equilibration run with similar parameters for cMD, five separate production runs for CFF were extended to 5000 ns.

The threshold energy for applying the ligand essential potential boost was also set to the upper bound in the LiGaMD simulation for both TEP and CFF. The selective boost potential was applied to the ligand. For the second boost potential applied to the system’s total potential energy other than the ligand’s essential potential energy, sufficient acceleration was obtained by setting the threshold energy to the lower bound. We aimed to observe ligand dissociation during equilibration while maintaining a low boost potential for accurate energetic reweighting; the $\left({\sigma }_{\mathrm{0}\mathrm{P}},{\sigma }_{0\mathrm{D}}\right)$ parameters were finally set to (4.0, 3.0 kcal/mol) for both systems (Figures S1 and S2). LiGaMD production simulation frames were saved every 0.4 ps for analysis.

Simulation Analysis.

The VMD and CPPTRAJ tools were used for simulation analysis.^36,37 The number of ligand dissociation and binding events ($N_{\mathrm{D}}$ and $N_{\mathrm{B}}$) and the ligand binding and unbinding time periods (${\tau }_{\mathrm{B}}$ and ${\tau }_{\mathrm{U}}$) were recorded from individual simulations for TEP (Tables 1, S1 and S2) and CFF (Tables 2, S3 and S4) simulations. With high fluctuations, (${\tau }_{\mathrm{B}}$ and ${\tau }_{\mathrm{U}}$) were recorded for only time periods longer than 1 ns.⁶ The 1D, 2D, and 3D PMF profiles as well as the ligand binding free energy were calculated through energetic reweighting of the LiGaMD simulations. The distance between the N4 atom of the C8 residue of the nucleic acid and the N7 atom of TEP and CFF was used as a reaction coordinate for calculating 1D PMF. The ligand heavy atom RMSD relative to the NMR bound pose was explored as a reaction coordinate for calculating the 1D PMF profile. The bin size was set to 1.0 Å for the atom distances and RMSD calculations. 2D PMF profiles were calculated for conformational changes in RNA upon ligand binding. We used the RMSD of the ligand relative to the NMR structure and the Solvent Accessible Surface Area (SASA) of nucleotide C27 on the RNA. The bin size for the latter was set to 5.0 Ų. The cutoff for the number of simulation frames in one bin was set to 500 for reweighing of 1D and 2D PMF profiles. The 3D PMF profiles of ligand displacements from the RNA host in the $X,Y$, and $Z$ directions were further calculated from the LiGaMD simulations. The bin sizes were set to 1 Å in the $X,Y$, and $Z$ directions. The cutoff of simulation frames in one bin for 3D PMF reweighting (ranging from 100 to 500 for three individual LiGaMD simulations) was set to the minimum number below which the calculated 3D PMF minimum will be shifted. Furthermore, the ligand binding free energies $(\mathrm{\Delta }G)$ were calculated using the reweighted 3D PMF profiles and binding kinetic rates by $\mathrm{\Delta }G=-RT\mathrm{l}\mathrm{n}\left(k_{\mathrm{o}\mathrm{f}\mathrm{f}}/k_{\mathrm{o}\mathrm{n}}\right)$, respectively.²⁰ The ligand dissociation and binding rate constants ($k_{\mathrm{o}\mathrm{n}}$ and $k_{\mathrm{o}\mathrm{f}\mathrm{f}}$) were calculated from the LiGaMD simulations, with their accelerations analyzed using Kramers’ rate theory.^20,24,38

Furthermore, structural clustering of the trajectories was performed for each 2.5 μs (TEP) and 5 μs simulation (CFF) using the Hierarchical Agglomerative clustering algorithm³⁹ in CPPTRAJ.³⁷ The frames were sieved at a stride of 1000 s for clustering. The remaining frames were assigned to the closest cluster afterward. The RMSD cutoff for clustering was set to 3.5 Å. The resulting structural clusters were reweighted to obtain energetically significant binding pathways for the ligand. In addition, the ligand dissociation and binding rate constants ($k_{\mathrm{o}\mathrm{n}}$ and $k_{\mathrm{o}\mathrm{f}\mathrm{f}}$) were calculated from the LiGaMD simulations, with their accelerations analyzed using the Kramers’ rate theory as described above (Table S5).

SPR Analysis of Theophylline and Caffeine Binding to the RNA Aptamer.

The kinetics and thermodynamics of interactions between the ligands and the RNA aptamer were first examined by using an experimental SPR technique. Using the established protocol, we immobilized the biotinylated TEP RNA aptamer on a streptavidin SPR chip and titrated different concentrations of ligands in solution (Figure 1A). From kinetic fitting, we determined that $k_{\mathrm{o}\mathrm{n}}=\left(1.92\pm 0.37\right)\times {10}^5{\mathrm{M}}^{-1}{\mathrm{s}}^{-1}$ and $k_{\text{off}}=(0.081\pm 0.45){\mathrm{s}}^{-1}$, resulting in a $K_{\mathrm{d}}=0.42\mu \mathrm{M}$. Notably, these values are consistent with the previously reported values.^2,40,41 For CFF, due to weak binding affinity to the TEP aptamer, $k_{\mathrm{o}\mathrm{n}}$ and $k_{\text{off}}$ could not be determined accurately from kinetic fitting (Figure 1C). However, using a steady-state binding affinity analysis, we were able to obtain $K_{\mathrm{d}}=9.1\mathrm{m}\mathrm{M}$ (Figure 1D). Using the same analysis for TEP yields $K_{\mathrm{d}}=0.31\mu \mathrm{M}$ (Figure 1B), which is quite consistent with the result from kinetic fitting. These results were also in agreement with previous observations that show the binding affinity for TEP is 10,000 times higher than it is for CFF.⁴⁰

LiGaMD Captured Repetitive Theophylline Binding and Dissociation in RNA Aptamer.

LiGaMD equilibration simulations of TEP captured the complete dissociation of the bound ligand into the bulk solvent. This was followed by rebinding to the RNA target site within 63 ns. We observed that C27 in RNA showed two distinct conformations during ligand binding. These conformations alternated between one where C27 was solvent exposed, called the “Out” state, and the other where it was buried inside the ligand pocket, called the “In” state.

Following the LiGaMD equilibration, six independent 2500 ns production simulations (“Sim1″-“Sim6”) (Table S1) were further performed on the RNA-TEP system with randomized initial atomic velocities. Three complete cycles of TEP unbinding and rebinding in three of the six LiGaMD simulations were observed (“Sim1″, “Sim2”, and “Sim3” in Table S1). Once boosted out of the pocket, the ligand-bound back into the RNA target site with a minimum RMSD of as low as ~2.4 Å relative to the native NMR structure (PDB ID: 1O15) (Figure 2A).²⁷ During “Sim1″, TEP bound to the target site during ~1500–1525 ns and then dissociated into a bulk solvent (Figure 2B). Two other rebinding and dissociation events were observed at ~2200 and ~2400 ns, respectively. During “Sim2″, TEP bound to the RNA in ~200–500 ns (Figure 2C) and dissociated quickly, followed by another rapid rebinding event during ~500–1100 ns. The ligand again approached the RNA pocket at ~1800 ns with fast dissociation. In “Sim3″, TEP was observed to leave the pocket and reassociated with the binding site at ~550, ~1480, and ~1700 ns (Figure 2D). The other simulations in the set for TEP binding to the RNA aptamer were able to capture the ligand binding events clearly, such as rebinding in “Sim4” at ~500 ns, “Sim5” for ~1000–1029 ns, and “Sim6” at ~500 ns and again at ~2000 ns (Figure S3D–F). However, these simulations were not able to capture three complete cycles of dissociation and rebinding of the ligand to the pocket. In order to further analyze the TEP and estimate the ligand kinetics, the first three simulations were used. These simulations provided reasonable statistics for the analysis of ligand binding and dissociation.

Next, we explored the correlation between ligand binding and conformational changes in the RNA. RMSD of the ligand relative to the NMR structure and solvent accessible surface area (SASA) of C27 in the aptamer were used as reaction coordinates to calculate a 2D free energy profile (Figure 3). Six low-energy conformations were identified in the 2D free energy profile of the TEP-RNA system (Figure 3A). Namely, the “Bound/Out”, “Bound/In”, “Intermediate/Out”, “Intermediate/In”, “Unbound/Out”, and “Unbound/In” states. When the C27 was solvent exposed, TEP diffused into the binding site, forming a “Bound/Out” complex with (TEP RMSD, C27 SASA) centered at (~2 Å, ~252 Ų) (Figure 3B). C27 could also be packed behind TEP in the RNA pocket in the “Bound/In” state with the (TEP RMSD, C27 SASA) centered at (~3 Å, 120 Ų) (Figure 3C). In the “Intermediate/Out” state the ligand hovered near the RNA surface, making transient local contacts at an average RMSD of ~25 Å, while the C27 remained extended and sustained a high SASA of ~240 Ų (Figure 3D). In the “Intermediate/In” state, a low energy well was identified corresponding to ~20 Å TEP RMSD and ~135 Ų SASA of C27 (Figure 3E). When C27 was solvent exposed with no ligand bound to the aptamer, a low-energy “Unbound/Out” state was identified with the (TEP RMSD, C27 SASA) centered at (~50 Å, ~250 Ų) (Figure 3F). In the “Unbound/In” state, the C27 was found to be fully buried inward with ~130 Ų SASA, and TEP was pushed out into the bulk solvent with ~50 Å RMSD (Figure 3G). As the RNA restructured itself to make room for the incoming ligand, we observed that flipping C27 from inside to outside the binding pocket played a key role in the recognition of TEP.

The Mg²⁺ ions stabilize the S-turn fold between nucleotides C22 and G26, resulting in a narrow fold with a 5 Å distance between G25-C22 (Figure S7A). This conformation allows C27 to flip outside and U24 to interact with A28 via π–π stacking, stabilizing the RNA pocket. In the ligand-unbound state, Mg²⁺ shifts toward the back of the fold, increasing the G25-C22 distance to up to ~12 Å and weakening the stacking interaction between U24 and A28 (Figure S7B), causing the RNA pocket to become larger. Mg²⁺ ions facilitate important RNA conformational changes.

LiGaMD Captured Binding and Dissociation of Caffeine to RNA.

In comparison, we performed longer LiGaMD simulations (i.e., 5000 ns) in order to capture both binding and unbinding of caffeine to the RNA.²⁸ The simulations were initiated similarly with CFF in the TEP binding pocket. The 63 ns equilibration simulation captured the dissociation and rebinding of CFF to the aptamer. CFF diffused into the TEP binding pocket during the ~5000 ns LiGaMD simulations with a binding conformation similar to TEP in the NMR structure (Figure 4A). In “Sim1″, CFF bound to the pocket only once at ~1000 ns (Figure 4B), stayed inside the pocket briefly for ~6.85 ns, and thereafter dissociated quickly. “Sim2” also captured CFF rebinding at ~2500 ns and subsequent dissociation at ~2542 ns (Figure 4C). In “Sim3″, we were able to capture 2 complete rebinding and 2 dissociation events of CFF at ~400 and ~800 ns (Figure 4D). During rebinding events, the CFF RMSD relative to the native bound pose reached a minimum of ~2.8 Å.

To further probe conformational changes in the CFF-RNA system, CFF RMSD and SASA of C27 were used to calculate 2D PMF (Figure 5A). Five low-energy states were identified, including the “Bound/Out”, “Unbound/Out”, “I1/In”, “I2/In”, and an “Unbound/In”, for which the (CFF RMSD, C27 SASA) centered at (~2 Å, ~242 Ų) (Figure 5B), (~50 Å, ~230 Ų) (Figure 5C), (~22 Å, ~130 Ų) (Figure 5D), (~10 Å, 140 Ų) (Figure 5E), and (~50 Å, ~120 Ų) (Figure 5F), respectively. A “Bound/Out” state was identified when CFF bound to the target site with high SASA of C27 in the range of ~230–260 Ų (Figure 5B). A low energy state referred to as “Unbound/Out,” sampled CFF far away from RNA in the bulk solvent with RMSD centered at ~50 Å, and C27 SASA sampled at ~230 Ų (Figure 5C). In the “I1/In” intermediate state the CFF was observed to be located at RMSD ~22 Å near the stem of the aptamer, with a low C27 SASA sampled at ~130 Ų (Figure 5D). We recorded a second intermediate “I2/In” (Figure 5E), where CFF occupied a site ~10 Å away from the target TEP binding pocket nestled in a groove made by C9, A10, G11, C21, G25, and C20 nucleotides via π-stacking interactions in the loop region. In the “Unbound/In” state centered with CFF RMSD at ~50 Å, and C27 SASA ~120 Ų (Figure 5F), the TEP binding pocket was observed to be tightly packed, preventing CFF from making its way to lodge inside. The two unbound states differed from each other in their sampled C27 conformations.

Ligand Binding Kinetic Rates and Free Energies Calculated from LiGaMD Were Consistent with Experimental Data.

LiGaMD simulations captured repetitive binding and dissociation of TEP to the RNA aptamer, which allowed us to calculate the ligand binding kinetic rate constants. The time periods for the ligand found in the bound state $\left({\tau }_{\mathrm{B}}\right)$ and unbound $\left({\tau }_{\mathrm{U}}\right)$ states during the simulation replicates were recorded for each simulation event (Tables S2 and S4). Without reweighting, the binding rate constant $\left(k_{\mathrm{o}\mathrm{n}}*\right)$ of TEP to the RNA aptamer was directly calculated from LiGaMD trajectories as (2.10 ± 0.23) × 10¹⁰ M⁻¹ s ⁻¹, whereas the dissociation rate constant $\left(k_{\mathrm{o}\mathrm{f}\mathrm{f}}*\right)$ was calculated as (1.94 ± 0.13) × 10⁶ s ⁻¹.

Next, we reweighted the LiGaMD simulations to calculate acceleration factors of TEP binding and dissociation processes (Table S5) and recovered the original kinetic rate constants using Kramers’ rate theory. The dissociation free energy barrier $\left(\mathrm{\Delta }{F}_{\text{off}}\right)$ increased to 5.65 ± 0.34 kcal/mol in the reweighted PMF profiles from 3.61 ± 0.26 kcal/mol in the modified PMF profiles (Table S5). The free energy barrier for ligand binding $\left(\mathrm{\Delta }{F}_{\text{on}}\right)$ changed to 1.23 ± 0.65 kcal/mol in the reweighted profiles from 0.96 ± 0.53 kcal/mol in the modified PMF profile (Table S5). Curvatures of the reweighed ($w$) and modified ($w*$, no reweighting) free energy profiles were calculated near the ligand Bound (“B”) and Unbound (“U”) low-energy wells and the energy barrier (“Br”), along with their ratio of apparent diffusion coefficients calculated from LiGaMD simulations with reweighting ($D$) and without reweighting (modified, $D*$) (Table S5). The reweighted $k_{\mathrm{o}\mathrm{n}}$ was calculated to be (2.72 ± 0.23) × 10⁵ M⁻¹ s⁻¹, being consistent with the corresponding experimental value of (1.92 ± 0.37) × 10⁵ M⁻¹ s ⁻¹. The reweighted $k_{\text{off}}$ was calculated to be (2.34 ± 0.07) × 10⁻³ s⁻¹ compared to our experimental value of (0.081 ± 0.45) s⁻¹. From our reweighted kinetic parameters, we can also calculate a binding affinity ($K_{\mathrm{d}}=k_{\text{off}}/k_{\text{on}}=0.009\mu \mathrm{M}$) as well as the ligand binding free energy $(\mathrm{\Delta }G\left.=-RT\mathrm{l}\mathrm{n}\left(k_{\text{off}}/k_{\text{on}}\right)=-6.35\pm 1.38\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{m}\mathrm{o}\mathrm{l}\right)$. Both of which are quite comparable with their corresponding experimental determined values ($K_{\mathrm{d}}=0.42\mu \mathrm{M}$ and $\mathrm{\Delta }G=-7.83\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{m}\mathrm{o}\mathrm{l}$) for (Table 1).

Multiple Ligand Binding and Dissociation Pathways Were Identified from LiGaMD Simulations.

We dwelled deeper into our LiGaMD simulation trajectories to understand the different pathways associated with ligand binding and dissociation of both TEP and CFF. TEP adopted three unique pathways of spontaneous binding to the RNA target site (Figure 6A–C). During all of our simulations, we observed that C27 could sample frequently between the “Out” and “In” conformations. In “Sim3″, when the aptamer was unbound with its C27 base fully extended toward the solvent (SASA ~ 250 Ų), TEP diffused from outside into the pocket through a pathway via π–π stacking interaction with the G26 nucleotide in the aptamer’s loop region (Figure 6A and Movie S1). As TEP approached the binding site, it oriented its five-membered ring toward the pocket such that the two-methyl group faced outside. This ensured that the N7 and N9 atoms of TEP were available to make interactions with C22 in the same plane (Figures 6A and 3C). In another binding pathway observed in Sim2 (“BP2”) (Figure 6B), TEP formed π -stacking interaction with the extended C27 in the “Out” conformation, as shown in the ligand trace (Figure 6B and Movie S2). In the third pathway (“BP3”), as TEP passed through the stem of the aptamer, it formed stacking interactions with nucleotides A5 and C30 (Figure 6C and Movie S3). A different nucleotide A28, located below the plane of C27, displayed a novel base-flipping mechanism as it guided the TEP to the target site through the gap between A28 and C27. In our LiGaMD simulations, the BP1 pathway had a dominant occurrence over BP2, followed by BP3, which was the least sampled by TEP.

The LiGaMD simulation trajectories were further used to study key insights into the dissociation of TEP from the aptamer. Following ligand binding, TEP started to rotate and change its orientation, causing the five-membered ring to reorient itself so that N7 and N9 faced the solvent. The ligand moved out of the binding pocket via stacking interaction with G26 (Figure 6D and Movie S4). When the ligand changed its orientation, G29 formed a non-native base pair with U6 and a hydrogen bond with the N3 atom of TEP. This caused the A5-G29 wobble base pair to undergo fluctuations, leading to the site becoming incompatible with accommodating the bound TEP. The movement of G29, accompanied by the inward movement of C27, pushed the TEP out of the binding site. This observation was also consistent with previous work in the literature, which describes the flipping mechanism of the C27 nucleotide in unbound aptamers.¹²

In our LiGaMD simulations of CFF, we observed a sequence of conformational transitions as the CFF associated with the aptamer. The extra methyl group on N7 made it difficult for CFF to bind to the RNA pocket. We observed only 1–2 binding and dissociation events in 4 out of the 8 3500–5000 ns LiGaMD simulations (Tables S3 and S4). CFF typically formed stacking interaction with the solvent-exposed C27 and diffused inward through a gap between G26 and A28 into the RNA target site (Figure 7A and Movie S5). CFF bound to the target site in an orientation similar to TEP, with its 5-member ring facing toward the pocket. This pathway was observed to be the BP2 taken up by TEP. For dissociation of CFF, we observed that the ligand made its way out into the solvent through the gap between G26 and A28 of the RNA, through which it first made its way inside (Figure 7B and Movie S6). The methyl group in the same plane of the ligand may not cause steric repulsion of the π–π stacking interaction between the planes.

We analyze the low-energy states of the ligand binding and dissociation pathways for TEP and CFF to gain a more intuitive understanding of the prerequisites for ligand binding (Figure S8). In the case of TEP binding (Figure S8A), when C27 is oriented outward in the solvent, we observed the Unbound/Out conformation. As the simulation progresses, TEP attempts to reach the binding site, leading to a decrease in the distance between RNA and TEP before binding occurs, sampling the Intermediate/Out state. Subsequently, the ligand enters the binding pocket through one of the three binding pathways, adopting the stable Bound/Out conformation; this ensures that the N7 and N9 atoms of TEP were available to make interactions with C22 in the same plane. This is followed by a conformational change in the TEP within the binding pocket as it rotates in the plane, ultimately leading to dissociation from the binding pocket. In the sampled Bound/In state, C27 transitions inward. This is followed by a conformational change in TEP within the binding pocket as it rotates in the plane, ultimately leading to dissociation from the binding pocket. In the sampled Bound/In state, C27 transitions inward. Upon binding, TEP undergoes rotation and a change in orientation, causing the five-membered ring to face the solvent. The ligand that exists in the binding pocket through stacking interactions with G26 and G29 forms a non-native base pair with U6 while hydrogen bonding with TEP, rendering the site incompatible. On the other hand, the CFF binding pathway to RNA involves the initial Unbound/Out conformation of RNA, which allows CFF to diffuse toward the binding site via binding pathway 2 (BP2) to Bound/Out conformation. As C27 shifts inward represented by the SASA values in the range of ~120–150 Ų, rendering the pocket unsuitable for accommodating the ligand, The Bound/Out state transitions to I1/In and eventually to an Unbound/In state as CFF dissociates. Additionally, we observe the binding of CFF to a different site approximately ~10 Å away from the binding pocket. With C27 positioned inward, the I2/IN state follows the dissociation pathway through the I1/In state, ultimately reaching the Unbound/In state (Figure S8B).

In our simulations, CFF spent the majority of its time prolonging stacking interactions with C27 outside of the pocket. This could be due to steric repulsions generated by the additional methyl group found on N7 of CFF, forcing it to adhere to the conserved site less frequently. When CFF acquired this pocket, it oriented itself so that the methyl group faced outward. This orientation allowed for the O2 atom of CFF to form an intermolecular hydrogen bond with the amino group on C22. TEP formed a hydrogen bond between N7 and the O4 atom of RNA U23. The nucleotides C8, A7, and G26 form a stable base triple to accommodate the ligand in the pocket (Figure S9B). Unlike TEP binding, the presence of CFF in the binding pocket disrupted the A7-C8-G26 base triple (Figure S9B). The disruption of this important base triple also aligned our simulation findings with literature explanations explaining the weaker interactions of CFF with the RNA aptamer.