Binding free-energy landscapes of small molecule binder and non-binder to FMN riboswitch: All-atom molecular dynamics

Abstract

A small and flexible molecule, ribocil A (non-binder) or B (binder), binds to the deep pocket of the aptamer domain of the FMN riboswitch, which is an RNA molecule. This binding was studied by mD-VcMD, which is a generalized-ensemble simulation method. Ribocil A and B are structurally similar because they are optical isomers to each other. In the initial conformation of simulation, the ligands and the aptamer were completely dissociated in explicit solvent. The aptamer–ribocil B binding was stronger than the aptamer–ribocil A binding, which agrees with experiments. The computed free-energy landscape for the aptamer–ribocil B binding was funnel-like, whereas that for the aptamer–ribocil A binding was rugged. When passing through the gate (named “front gate”) of the binding pocket, each ligand interacted with bases of the riboswitch by non-native π-π stackings, and the stackings restrained the ligand’s orientation to be advantageous to reach the binding site smoothly. When the ligands reached the binding site in the pocket, the non-native stackings were replaced by the native stackings. The ligand’s orientation restriction is discussed referring to a selection mechanism reported in an earlier work on a drug–GPCR interaction. The present simulation showed another pathway leading the ligands to the binding site. The gate (“rear gate”) for this pathway was located completely opposite to the front gate on the aptamer’s surface. However, the approach from the rear gate required overcoming a free-energy barrier regarding ligand’s rotation before reaching the binding site.

 Introduction

Computational drug discovery of RNA has been one of the emerging topics [1–6] and riboswitches are popular targets. Riboswitches are in the 5' or 3' untranslated regions of mRNA and regulate translation of mRNAs using small molecules that bind to riboswitches [7–10]. A riboswitch consists of two regions: An pptamer domain and an expression platform. Upon binding of a small ligand to the binding site of the aptamer domain, the structure of the expression platform changes, and then gene–expression process is switched on or off by the effect of the structural change induced in the expression platform. In general, the ligand-binding rate to a receptor (the riboswitch in this study) depends on the ligand concentration around the receptor. Therefore, the riboswitch has a role of switching translation on or off depending on the ligand concentration [11]. Because riboswitches have been found in pathogenic bacteria, a compound that inhibits binding of the natural ligand to the riboswitch can be a drug candidate [12].

Riboflavin (also known as vitamin B₂) is metabolized and flavin mononucleotide (FMN) is yielded in cell. The FMN riboswitch [8,10] (also known as the RFN element [13]) is a riboswitch of prokaryote and FMN binds to the aptamer domain. This binding results in suppression of riboflavin’s biosynthesis by inducing a structural change in its expression platform [9,10]. When the concentration of FMN rises in cell, the probability of complex formation of FMN and the aptamer domain increases. Thus, the FMN riboswitch works as a sensor providing feed-back control of the biosynthesis of riboflavin [13,14].

The structure of the aptamer domain of the Fusobacterium nucleatum FMN riboswitch has been studied experimentally [15–18]. The complex structure of ribocil B and the aptamer domain was solved by X-ray crystallography (PDB ID: 5c45; resolution 2.93 $\mathring{\mathrm{A}}$) [19]. Supplementary Figure S1 illustrates the ligand bound to the deep pocket of the aptamer domain. Portions of the aptamer domain are referred to as P1–P6 [20], as shown in Supplementary Figure S1a–c. Binding experiments of several compounds to the aptamer domain showed that ribocil B binds to the aptamer domain more strongly than ribocil A does, although ribocil A and B are very similar compounds (optical isomers) to each other (Figures 1a and b) [21]. Interestingly, while the binding modes of ribocil A and B to the aptamer domain are almost identical [21], an X-ray complex structure was not reported for ribocil A.

Figure 1 

Examined small ligands for mD-VcMD simulations: (a) Ribocil A (non-binder) and (b) ribocil B (binder), which are optical isomers mutually. Ligands are flexible in simulation because all the dihedral bonds are rotatable. Red-colored hydrogen atom labeled “H” in panels (a) and (b) are referred to as “$H_{diff}$ atom” for convenience and mentioned in main text. (c) Definition of three RCs: ${\lambda }^{\left(\alpha \right)}$, ${\lambda }^{\left(\beta \right)}$, and ${\lambda }^{\left(\gamma \right)}$. Method to set an RC (${\lambda }^{\left(h\right)};h=\alpha$, $\beta$, or $\gamma$) with using two atom groups $G_A^{\left(h\right)}$ and $G_B^{\left(h\right)}$ is explained in the main text and Supplementary Section 2. See Supplementary Table S1 for details of atom groups actually used in the present study. Atom groups $G_A^{\left(\alpha \right)}$, $G_A^{\left(\beta \right)}$, $G_A^{\left(\gamma \right)}$ and $G_B^{\left(\gamma \right)}$ are set in the aptamer domain, and $G_B^{\left(\alpha \right)}$ and $G_B^{\left(\beta \right)}$ are done in ribocil. The atom group $G_A^{\left(\alpha \right)}$ presented in green is split into two regions of the aptamer domain.

We define three directions to view the aptamer domain of the FNM riboswitch: the “front view” (Supplementary Figures S1a and d), the “side view” (Supplementary Figures S1b and e), and the “rear view” (Supplementary Figures S1c and f). Interestingly, the ligand, ribocil B, is visible in both the front and rear views (the magenta molecule in Supplementary Figures S1d and f), although the aptamer domain is usually illustrated by the front view in many papers. Supplementary Figure S2 displays the apo and holo forms of the aptamer domain, where ribocil B is not presented in the apo form (Supplementary Figure S2b). This figure indicates that the binding pocket forms a tunnel that is excavating the aptamer domain between the front and rear sides of the aptamer’s surface in both forms. The existence of the tunnel suggests that the ligand may reach the binding site from both the front and rear sides. Although these X-ray structures arouse researchers’ interest regarding the ligand binding process, the binding process is not understood well. The binding process is a temporary phenomenon where the ligand approaches the aptamer from a distant position, enters the deep pocket of the aptamer, and reaches the binding site in the pocket.

A molecular dynamics simulation can elucidate dynamic processes occurring in a biomolecular system. Especially, enhanced sampling (generalized-ensemble methods) reproduces rare and large conformational motions in the molecular system [22–24]. The multi-dimensional virtual-system coupled molecular dynamics (mD-VcMD) [25] is one of the enhanced sampling methods designed for investigating molecular binding processes. To perform an mD-VcMD simulation, Multiple reaction coordinates (RCs) are first set in advance. Then, the method enhances the motions in a space constructed by the multiple RCs. Importantly, a statistical weight (thermodynamic weight) equilibrated at a simulation temperature is assigned to all the sampled conformations (snapshots). Therefore, the ensemble of snapshots can be regarded a thermally equilibrated ensemble (canonical ensemble) of the snapshots. This method has been applied to various ligand–receptor systems [26–31]. We have proposed three variants of mD-VcMD so far [25]: The original, a subzone-based version, and genetic algorithm-guided mD-VcMD simulations.

In this paper we computationally investigate the ribocil–aptamer complex formation process focusing on two issues: i) whether ribocil B binds to the aptamer more strongly than ribocil A does and ii) whether ribocil approaches the binding site from the front side or the rear side of the aptamer surface. For this purpose, we first generate two computational systems consisting of the aptamer domain and each of the ribocil molecules in an explicit solvent. Because we focus on the binding process, the expression platform is excluded from the system. Then, we perform mD-VcMD starting the simulation from a completely unbound conformation, where ribocil is located far from the aptamer. The simulations sampled widely the conformational space for the two systems providing a structural ensemble of both systems: Then, we calculated the free-energy landscapes for the two systems, from which the two issues mentioned above are analyzed and discussed.

 Materials and Methods

In this work, mD-VcMD simulation is performed to obtain the conformational ensembles of the two systems. Each system consists of a receptor (the aptamer domain of the F. nucleatum FMF riboswitch) and a ligand (either ribocil A or B), and the ensembles were used to obtain the free-energy landscape. When the ligand is ribocil A, the system is referred to as “Ribo-A” and when it is ribocil B, it is done to as “Ribo-B”.

For conformational sampling, we use the subzone-based mD-VcMD, which is one of the three variants of the mD-VcMD methods [25], although the method is simply referred to as mD-VcMD in this paper. A statistical weight (i.e., thermodynamic weight at the simulation temperature) is assigned to each sampled conformation (i.e., snapshot). Then the ensemble of the weighted snapshots (i.e., canonical ensemble) is used to calculate various physical quantities at equilibrium.

Below, we first explain the simulation systems and next introduce multiple reaction coordinates (RCs). Sampling is done in the conformational space constructed by the multiple RCs. Then, we explain briefly the mD-VcMD sampling method, and lastly describe the procedures to calculate some physical quantities.

 Molecular System and the Initial Conformation of Simulation

The structure of the aptamer domain was taken from an X-ray crystallography (PDB ID: 6wjr; 2.7 Å resolution), which is the apo form (i.e., unbound form) of the aptamer domain [20]. Because we start the simulation from a completely dissociated conformation, the choice of the apo form is mandatory. We applied no structural restraints to the system. Therefore, the structures of the aptamer and each ribocil vary freely during the simulation.

The aptamer domain was put at the center of a cubic periodic boundary box filled by solvent ($97.9582\mathring{\mathrm{A}}\times 97.9582\mathring{\mathrm{A}}\times 97.9582\mathring{\mathrm{A}}$), and ribocil A or B was put at the corner of the box (Supplementary Figure S3a). Fifteen Mg ions were placed randomly in the bulk region of the box (Supplementary Section 1 for details). Additional ions, K and Cl, were introduced to set the ionic concentration to the experimental condition [20] (0.1 M), and to neutralize the net charge of the entire systems. The two systems consisted of the same atoms 92,628: 3,572 RNA atoms, 49 ligand atoms, 29,602 water molecules, 15 Mg ions, 133 K ions and 53 Cl ions. Then, an NPT simulation (temperature=300 K; pressure=1 atm) was performed to relax the system. The resultant box size was $97.3511\mathring{\mathrm{A}}\times 97.3511\mathring{\mathrm{A}}\times 97.3511\mathring{\mathrm{A}}$ for the Ribo-A system and $97.3304\mathring{\mathrm{A}}\times 97.3304\mathring{\mathrm{A}}\times 97.3304\mathring{\mathrm{A}}$ for the Ribo-B system, and the final structures were used for the initial ones of the mD-VcMD simulation. Supplementary Figure S3b shows that ribocil and the aptamer were completely dissociated. Details for generating the initial structures are explained in Supplementary Section 1.

The force field for the aptamer domain was from the ff99bsc0χOL3 (or χ_OL3) [32], while those for water molecules were from the 3-point optimal point charge (OPC3) water model [33]. The force field parameters for Mg ion were obtained from a recent study on divalent ions [34] and those for K and Cl ions were obtained from another recent study on monovalent ions [35].

The chemical structures for ribocil A and B are shown in Figures 1a and b. The force field parameters for ribocil A and B were modeled by us. First, the atomic partial charges of ribocil A or B were calculated quantum-chemically using Gaussian09 [36] at the HF/6-31G* level, followed by RESP fitting [37]. Then, the obtained atomic partial charges were incorporated with the general amber force field file 2 (GAFF2) [38], which was designed to be compatible with conventional AMBER force-fields.

 Three Reaction Coordinates (RCs) Introduced to Control System’s Motions

The mD-VcMD simulation enhances the conformational sampling in a space constructed by multiple RCs: ${\lambda }^{\left(h\right)}$ ($h=\alpha ,\beta ,\gamma ,\cdots$). Thus, the RCs should be set beforehand. In this study, a single RC, ${\lambda }^{\left(h\right)}$, is defined by the inter-centroid distance between two atom groups, $G_A^{\left(h\right)}$ and $G_B^{\left(h\right)}$, as shown in Supplementary Section 2, in which Supplementary Figure S4 presents a single RC schematically. Selection of RCs is important to increase the sampling efficiency, although RCs can be set arbitrarily in theory. Here, we imposed three conditions on the RCs: Variations of the RCs should control (i) ribocil approaching to (or departing from) the aptamer domain, (ii) ribocil’s orientation relative to the aptamer, and (iii) the opening (or closing) of the binding-pocket gate of the aptamer domain. Although various RCs can satisfy the three specific conditions, we introduced three RCs (${\lambda }^{\left(\alpha \right)}$, ${\lambda }^{\left(\beta \right)}$ and ${\lambda }^{\left(\gamma \right)}$) in a straightforward manner as shown in Figure 1c for both the Ribo-A and Ribo-B systems. The reason why we set the same RCs for the two systems is because the structures of ribocil A and B overlap very well to each other in the complex structures [21] as mentioned in the Introduction section.

The atom groups currently used (Figure 1c) are described in Supplementary Section 3 and Supplementary Table S1. Because three RCs were introduced in this study, “multi-dimensional RC (mD RC)” becomes “three-dimensional RC (3D-RC)” actually. We briefly explain the roles of the RCs. Two atom groups $G_A^{\left(\alpha \right)}$ and $G_B^{\left(\alpha \right)}$ are introduced to define the first RC ${\lambda }^{\left(\alpha \right)}$, and increment or decrement of ${\lambda }^{\left(\alpha \right)}$ corresponds to ribocil approaching to or departing from the aptamer, respectively. This relates to condition (i). The atom groups $G_A^{\left(\beta \right)}$ and $G_B^{\left(\beta \right)}$ define ${\lambda }^{\left(\beta \right)}$ and motions of ${\lambda }^{\left(\beta \right)}$ also related to condition (i). However, if ${\lambda }^{\left(\alpha \right)}$ decreases and ${\lambda }^{\left(\beta \right)}$ increases simultaneously, or if ${\lambda }^{\left(\alpha \right)}$ increases and ${\lambda }^{\left(\beta \right)}$ decreases, then the ligand rotates: Condition (ii) is satisfied. Last, $G_{\gamma }^A$ and $G_{\gamma }^B$ define the third RC ${\lambda }^{\left(\gamma \right)}$, and increment (or decrement) of ${\lambda }^{\left(\gamma \right)}$ induces opening (or closing) of the gate of the binding pocket: Condition (iii) is satisfied. Similar RC setting was used in our previous studies [27–30].

 mD-VcMD and Conformational Ensemble

The mD-VcMD samples the conformation in a region of 3D-RC space at 300 K, and then the 3D-RC region should be defined in advance by user. For this purpose, we set the minimum and maximum boundaries (denoted as ${\lambda }_{min}^{\left(h\right)}$ and ${\lambda }_{max}^{\left(h\right)}$, respectively) for each RC. Thus, the conformation moves in the range of $[{\lambda }_{min}^{\left(h\right)},{\lambda }_{max}^{\left(h\right)}]$. Importantly, the conformation should adopt both bound and unbound conformations in the range. The actual values of the boundaries are presented in Supplementary Table S2. We set ${\lambda }_{min}^{\left(h\right)}$ to $0.0\mathring{\mathrm{A}}$ for $h=\alpha$ and $\beta$. This value of zero is the smallest in theory. The ${\lambda }_{max}^{\left(h\right)}$ ($83.1178\mathring{\mathrm{A}}$) for $h=\alpha$ and $\beta$ is large enough to produce unbound conformations. Recall that ${\lambda }^{\left(\gamma \right)}$ is the gate width of the binding pocket of the aptamer domain. The value of ${\lambda }_{min}^{\left(\gamma \right)}$ ($5.0\mathring{\mathrm{A}}$) is small enough because it is likely that ribocil cannot not pass through the gate when ${\lambda }_{min}^{\left(\gamma \right)}=5.0\mathring{\mathrm{A}}$. The radius of a heavy atom is $2\mathring{\mathrm{A}}$ approximately. When ${\lambda }_{min}^{\left(\gamma \right)}=5.0\mathring{\mathrm{A}}$, only a space of $1\mathring{\mathrm{A}}$ remains, which is smaller than the diameter of a heavy atom ($5\mathring{\mathrm{A}}=2\mathring{\mathrm{A}}+2\mathring{\mathrm{A}}+1\mathring{\mathrm{A}}$). The value of ${\lambda }_{max}^{\left(\gamma \right)}$ (15.4248 $\mathring{\mathrm{A}}$), which corresponds to the maximum gate width, allowed ribocil to pass through the gate as shown later.

Next, we divided the ${\lambda }^{\left(h\right)}$ axis into small $n_{vs}^{\left(h\right)}$ zones: Supplementary Table S2 lists the actual value of $n_{vs}^{\left(h\right)}$. The lower and upper boundaries for the $i$-th zone along ${\lambda }^{\left(h\right)}$ are denoted as ${\lambda }_{low}^{\left(h\right)}\left(i\right)$ and ${\lambda }_{up}^{\left(h\right)}\left(i\right)$, respectively, whose values are listed in Supplementary Table S3. We refer to a zone along a single RC axis as a 1D zone. A relation of $[{\lambda }_{min}^{\left(h\right)},{\lambda }_{max}^{\left(h\right)}]=[{\lambda }_{low}^{\left(h\right)}\left(1\right),{\lambda }_{up}^{\left(h\right)}\left(n_{vs}^h\right)]$ is satisfied because the 1D zones were introduced by dividing the range of $[{\lambda }_{min}^{\left(h\right)},{\lambda }_{max}^{\left(h\right)}]$ into $n_{vs}^{\left(h\right)}$ zones. In the 3D-RC space, a zone is defined three-dimensionally: A set of $i$-th, $j$-th, and $k$-th 1D zones along the ${\lambda }^{\left(\alpha \right)}$-, ${\lambda }^{\left(\beta \right)}$- and ${\lambda }^{\left(\gamma \right)}$-axis, respectively, corresponds to a “3D zone” whose index is $(i,j,k)$. The ranges for this 3D zone are $[{\lambda }_{low}^{\left(\alpha \right)}\left(i\right),{\lambda }_{up}^{\left(\alpha \right)}\left(i\right)]$, $[{\lambda }_{low}^{\left(\beta \right)}\left(j\right),{\lambda }_{up}^{\left(\beta \right)}\left(j\right)]$ and $[{\lambda }_{low}^{\left(\gamma \right)}\left(k\right),{\lambda }_{up}^{\left(\gamma \right)}\left(k\right)]$ along the ${\lambda }^{\left(\alpha \right)}$-, ${\lambda }^{\left(\beta \right)}$-, and ${\lambda }^{\left(\gamma \right)}$-axis, respectively.

Suppose that the conformation of the system is currently in a 3D zone $(i,j,k)$ at time $=t$. We refer to this 3D zone as the “current zone”. In mD-VcMD, the conformation is confined in the current zone for a short interval, $\mathrm{\Delta }{\tau }_{trns}$, of simulation: A restoring force towards the current zone is applied to the system only when the conformation is flying to outside the current zone, and no bias force is applied to the conformation when the conformation is inside the current zone [25]. Then, at a time passage of $\mathrm{\Delta }{\tau }_{trns}$ (time $=t+\mathrm{\Delta }{\tau }_{trns}$), the conformation may transition from $(i,j,k)$ to another zone $(i^{\prime },j^{\prime },k^{\prime })$, which is one of the zones adjacent to $(i,j,k)$, with using a transition matrix $B$ [25]. This transition is called an “inter-zone transition”. Then the current zone is updated to $(i^{\prime },j^{\prime },k^{\prime })$. A matrix element $B_{(i,j,k)\to (i^{\prime },j^{\prime },k^{\prime })}$ is the transition probability from zone $(i,j,k)$ to zone $(i^{\prime },j^{\prime },k^{\prime })$. In the present study we set $\mathrm{\Delta }{\tau }_{trns}=20ps$, which corresponds to $1\times {10}^4$ simulation steps (a time step of simulation=$2fs$). We note that $\mathrm{\Delta }{\tau }_{trns}$ is not the simulation length but a part of it. Given a simulation run, the system experiences inter-zone transitions $L/\mathrm{\Delta }{\tau }_{trns}$ times during the run ($L$ is the simulation length): $L/\mathrm{\Delta }{\tau }_{trns}=100$ in the present study ($L=2ns$). As explained later, we executed many runs in parallel to increase the sampling efficiency.

Repeating the inter-zone transitions, the system’s conformation moves in the 3D-RC space. However, because the conformational space is too large to be sampled by a single simulation run, the mD-VcMD simulation is performed iteratively. When iteration $M$ has finished, a quantity $Q_{cano}^M(i,j,k)$, which is a probability of finding the system in the zone $(i,j,k)$, is calculated from the first to $M$-th iterations using the method introduced in our previous study [25]. Precisely speaking, $Q_{cano}^M(i,j,k)$ is a system’s canonical distribution (i.e., a thermally equilibrated distribution at the simulation temperature $T$) in a 3D zone $(i,j,k)$. The transition matrix $B$ for the $(M+1)$-th iteration is calculated from $Q_{cano}^M(i,j,k)$ whose protocol is given in our earlier studies [25,27–30]. Then the $(M+1)$-th iteration is started. The transition matrix for the first iteration is set exceptionally as: $B_{(i,j,k)\to (i^{\prime },j^{\prime },k^{\prime })}=1.0$, which means that the transition from the current zone to any adjacent zones occurs equally.

The simulation is terminated when a convergence criterion is satisfied for $Q_{cano}^M(i,j,k)$, and then $Q_{cano}^M(i,j,k)$ is expressed simply as $Q_{cano}(i,j,k)$. The convergence check of the resultant distribution is done using a function $E_{local}(i,j,k)$, which was introduced originally in an earlier report [25]. We briefly explain the meaning of $E_{local}(i,j,k)$. Suppose a cubic region consisting of $27(=3\times 3\times 3)$ zones in the 3D-RC space, where the indices $(i,j,k)$ of the cubic region are given as: $i=i_c-1,i_c,i_c+1$ for the axis ${\lambda }^{\left(\alpha \right)}$, $j=j_c-1,j_c,j_c+1$ for the axis ${\lambda }^{\left(\beta \right)}$, and $k=k_c-1,k_c,k_c+1$ for the axis ${\lambda }^{\left(\gamma \right)}$. The zone $(i_c,j_c,k_c)$ is the central zone of the cubic region, and thus the other 26 zones surround the central zone. During the mD-VcMD simulation, transitions among zones occur, and the transition involves those among the 27 zones in the cubic region. $E_{local}(i,j,k)$ has a property that the more the detailed balance condition is satisfied in the transitions around the zone $(i_c,j_c,k_c)$, the smaller the $E_{local}(i,j,k)$ becomes around the zone $(i_c,j_c,k_c)$. Then, we judged that if an inequality $E_{local}(i,j,k)<E_{local}^0$ is satisfied in a volume that involves important regions of the 3D-RC space (i.e., high $Q_{cano}(i,j,k)$ regions), we judged that the simulation is converged. We have set $E_{local}^0=0.25$ in studies on molecular binding [26–28,30] and found that the native complex structure corresponded to the lowest free-energy region computed.

A snapshot of the system is stored once every $1\times {10}^5$ steps ($0.2$ ns) in all iterations executed. After convergence, a thermodynamic weight at the simulation temperature $T$ is assigned to each snapshot using the probability $Q_{cano}(i,j,k)$ [25]. Therefore, the snapshots construct a thermally equilibrated conformational ensemble (canonical ensemble) at $T$ ($300K)$. Remember that no force is applied to the system’s conformation in each period of $\mathrm{\Delta }{\tau }_{trns}$ of the simulation. Therefore, all snapshots detected in a zone from iterations 1 to $M$ were treated with an equal weight [25].

To increase further the sampling efficiency, we performed $N_{run}$ (=1,024) runs of the mD-VcMD simulation in parallel for each iteration, where the runs started from different conformations. If the $M$-th iteration has finished, we have snapshots sampled from iterations 1 to $M$. Then, the initial conformations of the $N_{run}$ runs for the $(M+1)$-th iteration were selected from those snapshots so that they are distributed as evenly as possible in the 3D-RC space. Exceptionally, all the runs for the first iteration start from a single conformation obtained from the NPT simulation (Supplementary Figure S3b).

Our previous studies [39,40] have shown that an ensemble consisting of multiple short trajectories converges to a canonical ensemble when the initial conformations of the short trajectories distribute widely in the conformational space. This means that a trajectory, generated by connecting the short trajectories in an arbitrary order, can be a single simulation trajectory because the detailed balance condition is satisfied at the trajectory connection points [40]. Here, we refer to a short trajectory from the $i$-th short run as ${TR}_{short}^i$ and its length as $L\left({TR}_{short}^i\right)$. The trajectory generated from the short trajectories is denoted as ${TR}_{sum}$: Its length is $L\left({TR}_{sum}\right)={\sum }_i^{}L\left({TR}_{short}^i\right)=N_{run}\times L\left({trj}_{short}^i\right)$. For the present study, $L\left({TR}_{sum}\right)=92.16\mu s=2ns\times 1024runs\times 45iterations$, where 45 is the number of iterations done in the present study for each system as explained later. Now, imagine another single trajectory, referred to as ${TR}_{single}$, that was generated from an actual single run. We assume that $L\left({TR}_{single}\right)=L\left({TR}_{sum}\right)$. Importantly, the volume covered by ${TR}_{sum}$ in the conformational space is larger than that by ${TR}_{single}$ because ${TR}_{sum}$ involves conformation jumps at the connection points between two short trajectories. Therefore, the ensemble of short runs, ${TR}_{sum}$, has a better sampling efficiency than the long single simulation, ${TR}_{single}$, has.

For both the Ribo-A and Ribo-B systems, we set the simulation length of each of $N_{run}$ runs to $1\times {10}^6$ steps: $L\left({TR}_{short}^i\right)=1\times {10}^6\times 2fs=2.0ns$. Consequently, the total length for an iteration was $2.048\mu s$ (=$2.0ns\times 1,024$). We performed 45 iterations for the present study, which corresponds to $92.16\mu s$ ($2.048\mu s\times 45$) in all for each system: $L\left({TR}_{sum}\right)=92.16\mu s$. As mentioned above, a snapshot was stored every $1\times {10}^5$ steps ($0.2ns$) in each run. Therefore, 460,800 ($92.16\mu s/0.2ns$) snapshots were stored for each system. We calculated the distribution functions of various quantities at the simulation temperature ($T=300K$) from the resultant ensemble of snapshots.

The mD-VcMD is a sampling method in that a single zone is sampled by many runs performed in parallel and by many iterations. This multiple sampling increases the statistics. Suppose that some stable conformations, which are distinguishable one another, exist in a zone (denoted as $z_x$) and that energy barriers exist among the stable conformations in the zone $z_x$. Then, if a simulation is confined in $z_x$, the simulation samples only one of the stable conformations. In contrast, the mD-VcMD simulation moves among different zones (inter-zone transitions). Then this simulation may sample those stable conformations by taking a route that avoids the potential barriers in the outside of the zone $z_x$. Furthermore, different runs can sample different stable conformations in $z_x$. All snapshots from those runs are combined in a single ensemble to analyze the system.

We also emphasize that the mD-VcMD simulation is a global search method: Any thermodynamically probable conformation is sampled by varying ${\lambda }^{\left(\alpha \right)}$, ${\lambda }^{\left(\beta \right)}$ and ${\lambda }^{\left(\gamma \right)}$. This is because any relative positioning between the ligand and the receptor can be assigned to a 3D zone.

The mD-VcMD algorithm was first implemented to a MD simulation program omegagene/myPresto [23] and next to GROMACS [41] with the PLUMED plug-in [42], which was used for the current study. The simulation conditions were the following: LINCS [43] to constrain the covalent-bond lengths related to hydrogen atoms, the Nosé–Hoover thermostat [44,45] to control the simulation temperature, a time-step of 2 fs ($\mathrm{\Delta }t=2$ fs), and simulation temperature of $300K$. The zero-dipole summation method [46–48] was used to calculate quickly and accurately the long-range electrostatic interactions. This method calculates the electrostatic interactions without an artifact caused by imposing a periodicity to a non-periodic system originally (a biological system for instance) [49].

Comparison of the execution time between the zero-dipole summation method and the Particle Mesh Ewald Sum method was done in TIP3P water systems [47], which indicated that the zero-dipole summation is comparable or faster than the Particle Mesh Ewald Sum. For instance, the time of the zero-dipole summation with $\alpha =0$ and $r_c=12\mathring{\mathrm{A}}$ ($\alpha$ is the damping factor and $r_c$ is the cutoff length) and that of the Particle Mesh Ewald Sum with the real part $r_c=8\mathring{\mathrm{A}}$ without the Fourier part evaluation become comparable with increasing the number of parallel processors. See supplement file of Ref. 47 for details of this comparison. Comparison of the execution time between the zero-multipole summation method, which is a developed version of the zero-dipole summation method, and the PME using several protein systems in Gromacs was done [50], which also showed the superiority of our method.

All simulations were performed on the TSUBAME3.0 supercomputer at the Tokyo Institute of Technology using GP-GPU. A single run ($1\times {10}^6$ steps) took about 40 min using a single GPU. One iteration of mD-vcMD consisted of 1,024 runs in the present study. So, if 1,024 GPUs are available at the same time, one iteration ends by 40 min.

 Some Physical Quantities to Analyze Binding Mechanism

The mD-VcMD simulation produces a canonical ensemble of conformations (snapshots) at $300K$. This method assigns a thermodynamic weight at $300K$ to any snapshot obtained from the sampling [25]. Therefore, using this ensemble, one can calculate a spatial density of the ligand’s centroid around the receptor at 300 K using the ensemble. Besides we present some physical quantities, which are used to analyze the ligand binding process in this paper. Details for these quantities are explained in Supplementary Section 4.

Here we briefly explain the notations used to express the spatial density: See Supplementary Subsection 4.1 for details. The 3D real space (not the 3D-RC space) is divided into cubes, whose size is $L_c\times L_c\times L_c$, where $L_c$ is the length of sides of the cube, and the center of a cube is denoted as ${\mathit{r}}_{cube}$. $L_c$ is usually set to $2\mathring{\mathrm{A}}$ and rarely to $0.5\mathring{\mathrm{A}}$, and if no explanation is given, $L_c$ is $2\mathring{\mathrm{A}}$ in this paper. The spatial density ${\rho }_{CM}\left({\mathit{r}}_{cube}\right)$ is the density of the ligand’s centroid ${\mathit{r}}_{CM}$ in the cube centered at ${\mathit{r}}_{cube}$.

Next, we defined two molecular orientation vectors ${\mathit{e}}_{\uparrow }$ and ${\mathit{e}}_{\leftarrow }$ for ribocil (Supplementary Subsection 4.2). The ${\mathit{e}}_{\uparrow }$ is a unit vector parallel to a vector ${𝐯}_{\uparrow }$, which is defined from the centroid of Ring 2 and Ring 3 to that of Ring 1 and Ring 4. The other vector ${\mathit{e}}_{\leftarrow }$ is a unit vector parallel to a vector ${𝐯}_{\leftarrow }$, which is defined from the centroid of Ring 3 to that of Ring 2. Supplementary Figure S5 shows visually ${\mathit{e}}_{\uparrow }$ and ${\mathit{e}}_{\leftarrow }$. In general, an expectation value of a quantity $\Omega$ in the cube centered at ${\mathit{r}}_{cube}$ is denoted as $<\Omega \left({\mathit{r}}_{cube}\right)>$. Expression for $<\Omega \left({\mathit{r}}_{cube}\right)>$ is given in Supplementary Equation S4. Then, setting as $\Omega ={\mathit{e}}_{\leftarrow }$ or ${\mathit{e}}_{\uparrow }$, vector fields for the ribocil’s molecular orientation, $<{\mathit{e}}_{\uparrow }\left({\mathit{r}}_{cube}\right)>$ and $<{\mathit{e}}_{\leftarrow }\left({\mathit{r}}_{cube}\right)>$, are calculated.

The native complex structure (PDB ID: 5c45) shows three $\pi$-$\pi$ stackings (native stacking) between aptamer’s bases and ribocil’s rings (see Supplementary Figure S6 and Supplementary Table S4). To analyze the stacking, we introduced a quantitative procedure to assign the stacking to snapshots: See Supplementary Subsection 4.3. We checked snapshots and found that eight stackings are frequently formed in the complex formation process (Supplementary Table S4), three of which were the native stackings and the other five were non-native stacking. Supplementary Figure S7a–c illustrates snapshots with native stacking from the front-view and Figure S8a–c does the side-view of Figure S7a–c. Supplementary Figure S7d–f illustrates those with non-native stacking from the front view, and Figure S8d–f are the side-view of Figure S7d–f. We introduced two quantities: $N_{native}^{\left(stack\right)}$ and $N_{non-native}^{\left(stack\right)}$, which are respectively the numbers of native and non-native stacking patterns in a snapshot (see Supplementary Subsection 4.3).

We note that the non-native stackings are commonly formed before the ligand binds to the binding site experimentally determined. This does not contradict the fact that the ligand approaches the binding site gradually from outside of the binding pocket.

 Results

We performed 45 iterations of mD-VcMD and saved 460,800 snapshots as mentioned in the Methods and Materials section. In the Results section, we first show that the system’s distribution (the free-energy landscape) in the 3D-RC space covers both bound and unbound conformations. Next, we compute the distribution of ribocil in the real 3D space (i.e., not in the 3D-RC space) and demonstrate that robocil B at the highest-density spot (the lowest free-energy position) is in the deep pocket of the aptamer domain, whereas that of ribocil A is not. This means that ribocil B binds to the aptamer more strongly than ribocil A does. We then analyze the binding mechanism from the obtained conformational ensemble.

 Conformational Distribution in 3D-RC Space

Figure 2 illustrates $Q_{cano}(i,j,k)$, which is the equilibrated distribution of the system at 300 K in the 3D-RC space (not in the 3D real space). As explained in the Materials and Methods section, the three indices $i$, $j$ and $k$ specify the position of a zone in the 3D-RC space: See Supplementary Table S3 for the correspondence between the zoon index and the value of ${\lambda }^{\left(h\right)}$. Because $Q_{cano}(i,j,k)$ is the probability at the zone, we normalized it for both systems as:

Figure 2 

System’s canonical distribution $Q_{cano}^M(i,j,k)$ presented in the 3D-RC space. $Q_{cano}^M(i,j,k)$ is computed from the 1’st to $M$-th iterations ($M\le 45$). Panels (a) and (b) are, respectively, distributions for Ribo-A and Ribo-B systems. (c) The close-up of the highest density region of panel (b). Three coordinate axes are indices $i$, $j$, and $k$, which specify a 3D zone. The position of the zone in the 3D-RC space is: ${\lambda }^{\left(\alpha \right)}=\{{\lambda }_{low}^{\left(\alpha \right)}\left(i\right)+{\lambda }_{up}^{\left(\alpha \right)}\left(i\right)\}/2$, ${\lambda }^{\left(\beta \right)}=\{{\lambda }_{low}^{\left(\beta \right)}\left(j\right)+{\lambda }_{up}^{\left(\beta \right)}\left(j\right)\}/2$ and ${\lambda }^{\left(\gamma \right)}=\{{\lambda }_{low}^{\left(\gamma \right)}\left(k\right)+{\lambda }_{up}^{\left(\gamma \right)}\left(k\right)\}/2$. See Table S3 of SI for values of ${\lambda }_{low}^{\left(h\right)}\left(m\right)$ and ${\lambda }_{up}^{\left(\alpha \right)}\left(m\right)$ ($m=i$, $j$, or $k$). Snapshots from four sites $c_1–c_4$ of panels (a) and (b) are shown in Figure 3.

Figure 3 

Conformations picked from sites $c_i$ ($i=1,\dots ,4$) whose positions in the 3D-RC space are shown in Figure 2. Ribocils and aptamer are presented in magenta and light cyan, respectively. Panels (a), (b), (c) and (d) display snapshots from the Ribo-A system, and panels (e), (f), (g) and (h) so those from the Ribo-B system. Label “$c_i$” is presented near the panel identifier. Two conformations are picked from either site $c_2$ or $c_3$. Red arrows indicate direction of the front view (see Supplementary Figure S1).

${\sum }_{i=1}^{n_{vs}^{\left(\alpha \right)}}{\sum }_{j=1}^{n_{vs}^{\left(\beta \right)}}{\sum }_{k=1}^{n_{vs}^{\left(\gamma \right)}}{Q}_{cano}(i,j,k)=1$. (1)

This normalization is important to compare $Q_{cano}(i,j,k)$ between the Ribo-A and Ribo-B systems.

Figures 2b and c demonstrate that the Ribo-B system had a remarkable high-density spot (red-contour region involving the site $c_1$ in Figure 2b) in the 3D-RC space. We note that the red-contour region was surrounded by the blue-contour region and the blue-contour region was done by the magenta-contour region. In contrast, no such a remarkable high-density spot was found in the Ribo-A system (Figure 2a) with the unified density normalization (Equation 1), and the magenta-counter region occupied a large volume of the 3D-RC space. Therefore, the free-energy landscape for the Ribo-B system follows a funnel-like shape in the 3D-RC space, while that for the Bibo-A system is more rugged.

Conformations taken from sites $c_1–c_4$ in Figure 2 are displayed in Figure 3a–d for the Ribo-A system and in Figure 3e–h for the Ribo-B system. Note that the ribocil molecule from $c_1$ is in the binding pocket of the aptamer domain for both systems (Figures 3a and 3e). This is because the structures from $c_1$ had similar RC values for both systems (Figure 2). Because a high density was not assigned to $c_1$ for the Ribo-A system, the structures from $c_1$ of ribocil A were thermodynamically less stable than those of ribocil B. In the next subsection, we analyze details of the structures in the binding pocket.

Figures 2 and 3 demonstrates that mD-VcMD sampled widely the 3D-RC space for both systems: ribocil from $c_1$ was in the binding pocket, and ribocil from $c_4$ was completely dissociated from the aptamer domain (Figures 3d and 3h). Ribocil from $c_2$ was around the entrance of the binding pocket (Figures 3b and 3f), while the structures of ribocil from $c_3$ was outside of the binding pocket with contacting the aptamer domain. Interestingly, we found ribocil near the front and rear surfaces of the aptamer domain. This suggests that ribocil may enter the binding pocket from either side of the tunnel (Supplementary Figure S2).

The convergence of the resultant distribution $Q_{cano}^M(i,j,k)$ was checked using $E_{local}(i,j,k)$. See also Supplementary Section 5 and Supplementary Figure S9, which demonstrates that $Q_{cano}^M(i,j,k)$ converged in both the bound and unbound regions of the 3D-RC space during iterations 1–40 ($M=1,\dots ,40$) for both systems. We proceeded the simulation up to iteration 45 to increase the number of snapshots to be saved in the conformational ensemble and set as ${Q_{cano}(i,j,k)=Q}_{cano}^M(i,j,k)$.

 Density of Ribocil Around Aptamer Domain and Density Tunnel

Figure 4 demonstrates the spatial density ${\rho }_{CM}\left({\mathit{r}}_{cube}\right)$ of the ribocil’s centroid around the aptamer domain in the real 3D space, where the cube size was set to $L_c=2\mathring{\mathrm{A}}$ (see Supplementary Equation S1). A high-density region of ${\rho }_{CM}=0.0025{\mathring{\mathrm{A}}}^{-3}$ (red contour region) was only found in the Ribo-B system. This region is referred to as “the highest-density spot” to indicate clearly that this density is the highest, and then we denote this density as ${\rho }_{max}$ (i.e., ${\rho }_{max}=0.0025{\mathring{\mathrm{A}}}^{-3}$). Figure 4c demonstrates that the highest density spot is in the binding pocket for the Ribo-B system. The region of ${\rho }_{CM}={\rho }_{max}/10$ (blue contour region) surrounded the highest density spot, and the regions of ${\rho }_{CM}={\rho }_{max}/100$ (magenta contour region) surrounded the blue contour region. This spatial pattern of ${\rho }_{CM}$ relates to the funnel-like nature of Figure 2b. We note that the highest-density spot deviated by about $4\mathring{\mathrm{A}}$ from the X-ray position (small black sphere) to a direction parallel to the red arrow in Figure 4d. We discuss this shift in the next subsection.

Figure 4 

(a) Front and (b) side views of spatial density ${\rho }_{CM}\left({\mathit{r}}_{cube}\right)$ of ribocil A around the aptamer domain. Label $Pi$ ($i=1,\dots ,6$) is an identifier assigned to $Pi$ helix. Label “$Pi-Pi+1$” indicates junction region between $Pi$ and $Pi+1$. Similarly, front (c) and (d) side views of ${\rho }_{CM}\left({\mathit{r}}_{cube}\right)$ of the ribocil B. Iso-density levels are colored as indicated in inset, and ${\rho }_{max}=0.0025{\mathring{\mathrm{A}}}^{-3}$. The density is normalized for both systems: ${\sum }_{{\mathit{r}}_{cube}}^{}{\rho }_{CM}\left({\mathit{r}}_{cube}\right)\times {L_c}^3=1$. Light gray model represents aptamer domain, and green model does ribocil B in the X-ray structure (PDB ID: 5c45; P4 loop is modeled). Centroid of ribocil B in the X-ray structure is shown by black sphere. The structure of ribocil A is not determined experimentally. Black broken-line rectangle is density tunnel excavated from the front to rear surfaces of aptamer domain. Front and rear gates, from which ribocil enter the binding pocket, are shown in panels (a) and (b). See Supplementary Figure S1 for P4–P6.

For the Ribo-A system, no spot of ${{\rho }_{CM}=\rho }_{max}$ was found, and the binding pocket was partly occupied by the blue contours, although the volume for the blue contours was smaller than that for the Ribo-B system (Figure 4b). These features of ${\rho }_{CM}$ shows that the free-energy landscape of the Ribo-A system is more rugged than that for the Ribo-B system. Therefore, Figure 4 is consistent to Figure 2.

Here we emphasize that ${\rho }_{CM}$ was normalized as shown in Supplementary Equation S3. Therefore, we conclude that probability of ribocil B in the binding pocket was larger than that of ribocil A in the pocket. In other words, binding of ribocil B to the aptamer stronger than that of ribocil A.

Remember that the experimental structure of the aptamer domain has a tunnel (Supplementary Figure S2) in either the apo or holo form. In fact, we found a density tunnel as illustrated in Figures 4b and d: The magenta contours (${\rho }_{CM}={\rho }_{max}/100$) excavated the aptamer domain between the front and rear surfaces (the broken-line rectangle in Figures 4b and d). Therefore, ribocil can reach the binding position (the position in the X-ray complex) by passing gates on the front and rear surfaces of the aptamer domain. The gate (named “front gate”) on the front surface is a cleft formed by the P4–P5 junction and the P4–P5 junction of the aptamer (Figure 4a). The gate (“rear gate”) on the rear surface is a cleft opposite to the front gate (Figure 4b).

Because ${\rho }_{CM}\left({\mathit{r}}_{cube}\right)$ in Figure 4 was presented with $L_c=2\mathring{\mathrm{A}}$, this figure did not show a fine structure of ${\rho }_{CM}$ smaller than $2\mathring{\mathrm{A}}$. Then, to view the fine structure, we computed ${\rho }_{CM}$ with $L_c=0.5\mathring{\mathrm{A}}$. Figure 5 demonstrates a free-energy barrier (broken line) in both the Ribo-A and Ribo-B systems, whereas no clear free-energy barrier was seen in Figure 4. On the other hand, the number of snapshots detected in cubes decreases with decreasing $L_c$ (i.e., statistics decreases). Although we examined $L_c$ smaller than $0.5\mathring{\mathrm{A}}$, the spatial patterns of ${\rho }_{CM}\left({\mathit{r}}_{cube}\right)$ became bumpy (data not shown).

Figure 5 

${\rho }_{CM}\left({\mathit{r}}_{cube}\right)$ with setting $L_c=0.5\mathring{\mathrm{A}}$ focusing on the density tunnel (magenta contours of Figure 4) (a) Ribo-A system and (b) Ribo-B system.

 Orientation of Ribocil Around Aptamer Domain

Figure 6 illustrates the spatial patterns of four quantities regarding the molecular orientation of ribocil B: $<{\mathit{e}}_{\uparrow }\left({\mathit{r}}_{cube}\right)>$, $<{\mathit{e}}_{\leftarrow }\left({\mathit{r}}_{cube}\right)>$, ${SP}_{\uparrow }\left({\mathit{r}}_{cube}\right)$ and ${SP}_{\leftarrow }\left({\mathit{r}}_{cube}\right)$. See Supplementary Subsection 4.2 for definition of these quantities as well as ${\mathit{e}}_{\uparrow }$, ${\mathit{e}}_{\uparrow }^{(X-ray)}$, ${\mathit{e}}_{\leftarrow }$ and ${\mathit{e}}_{\leftarrow }^{(X-ray)}$. The quantities for ribocil A are presented in Supplementary Figure S10. Both of Figure 6 and Supplementally Figure S10 indicate that the molecular orientation of ligand was randomized outside the binding pocket because the number of black arrows, which satisfied the inequality $|<{\mathit{e}}_{\alpha }\left({\mathit{r}}_{cube}\right)>|\ge 0.6$ ($\alpha =\uparrow$ or $\leftarrow$), was small outside the binding pocket. In contrast, significantly large number of arrows was found in the binding pocket. These results indicate that for both systems, the ligand’s molecular orientation ordered when the ligand entered the gate of the binding pocket from the outside.

(a) $<{\mathit{e}}_{\uparrow }\left({\mathit{r}}_{cube}\right)>$ and ${SP}_{\uparrow }\left({\mathit{r}}_{cube}\right)$ for Ribo-B system. (b) $<{\mathit{e}}_{\leftarrow }\left({\mathit{r}}_{cube}\right)>$ and ${SP}_{\leftarrow }\left({\mathit{r}}_{cube}\right)$ for the Ribo-B system. In both panels, vectors $<{\mathit{e}}_{\alpha }\left({\mathit{r}}_{cube}\right)>$ ($\alpha =\uparrow$ or $\leftarrow$) with $|<{\mathit{e}}_{\alpha }\left({\mathit{r}}_{cube}\right)>|\ge 0.6$ are shown by black arrows. Similarly, sites with ${SP}_{\alpha }\left({\mathit{r}}_{cube}\right)\ge 0.6$ and ${SP}_{\alpha }\left({\mathit{r}}_{cube}\right)\le -0.6$ are shown, respectively, by magenta- and cyan-colored contours. Density tunnel is divided into the front and rear portions, which are indicated by broken-line rectangles. See also Supplementary Figure S10 for $<{\mathit{e}}_{\alpha }\left({\mathit{r}}_{cube}\right)>$ and ${SP}_{\alpha }\left({\mathit{r}}_{cube}\right)$ of Ribo-A system.

Figure 6 

Now, we divide the density tunnel (broken-line rectangle in Figure 4) into two portions, front and rear portions, as shown in Figure 6a. For the Ribo-B system, apparently, ${\mathit{e}}_{\uparrow }$ in the front portion of the tunnel tended to be parallel to ${\mathit{e}}_{\uparrow }^{(X-ray)}$ as indicated by arrows of $<{\mathit{e}}_{\uparrow }\left({\mathit{r}}_{cube}\right)>$ and the magenta contours (${SP}_{\uparrow }\left({\mathit{r}}_{cube}\right)\ge 0.6$) of Figure 6a. This tendency was also found in the Ribo-A system (Supplementary Figure S10a). In contrast, for both systems, ${\mathit{e}}_{\uparrow }$ in the rear portion tended to be antiparallel to ${\mathit{e}}_{\uparrow }^{(X-ray)}$ as shown by arrows of $<{\mathit{e}}_{\uparrow }\left({\mathit{r}}_{cube}\right)>$ and the cyan contours (${SP}_{\uparrow }\left({\mathit{r}}_{cube}\right)\le -0.6$). Interestingly, the magenta and blue contour regions switched from one to the other sharply at the boundary of the front and rear portions. Importantly, this boundary corresponds to the free-energy barrier of ${\rho }_{CM}\left({\mathit{r}}_{cube}\right)$ (the broken line in Figure 5). It is likely that the ligand rotation at the boundary accompanies a free-energy cost.

Next, we discuss the spatial patterns of ${\mathit{e}}_{\leftarrow }$. Figure 6b and Supplementary Figure S10b show that the regions of ${SP}_{\leftarrow }\left({\mathit{r}}_{cube}\right)\ge 0.6$ dominated the density tunnel, whereas some small regions of ${SP}_{\leftarrow }\left({\mathit{r}}_{cube}\right)\le -0.6$ were found in the tunnel for both systems. Therefore, the spatial patterns for ${\mathit{e}}_{\leftarrow }$ did not show a clear transitional behavior at the boundary between the front and rear portions of the density tunnel. We concluded that the transitional motion of the ribocil’s molecular orientation at the boundary is related to flipping motions of ${\mathit{e}}_{\uparrow }$, which is represented schematically in Supplementary Figure S11.

Last of this subsection, we note that both vectors ${\mathit{e}}_{\uparrow }$ and ${\mathit{e}}_{\leftarrow }$ tended to be parallel to ${\mathit{e}}_{\uparrow }^{(X-ray)}$ and ${\mathit{e}}_{\leftarrow }^{(X-ray)}$, respectively, in the front portion of the density tunnel (magenta contours in both Figure 6 and Supplementary Figure S10). Therefore, we presume that there is a factor to maintain the molecular orientation. We discuss this point in the next subsection.

 Stacking Between Aptamer Domain and Ribocil

In this and next subsections, we analyze intermolecular interactions that restrains the molecular orientation of ribocil in the binding pocket of the aptamer. A study on the X-ray complex structures of the aptamer domain and some ligands has reported that $\pi$- $\pi$ stackings between the ligands and the aptamer are formed to stabilize the complex structure [17]. Interestingly, three bases (A48, G62 and A85) of the aptamer are conserved in the stacking with different ligands. Our simulation showed that eight $\pi$-$\pi$ stackings were formed frequently when the ligand was in the binding pocket (Supplementary Table S4) in both Ribo-A and Ribo-B systems: Three of them were native stackings (Supplementary Figures S7a–c) and the other five were non-native stackings (Supplementary Figures S7d–f). Below, we analyzed the stacking quantitatively.

To pick up snapshots in that ribocil was in the binding pocket, we calculated root-mean-square-deviation of ribocil ($rms{d}_{rib}$) between a snapshot and the X-ray complex structure (PDB ID: 5c45). Remember that the aptamer part of snapshots had been superimposed to that of the X-ray structure (Supplementary Subsection 4.1). The $rms{d}_{rib}$ was calculated simply as: $rms{d}_{rib}=[{N_{rib}}^{-1}{\Sigma }_k{({\mathit{r}}_k-{\mathit{r}}_k^{(X-ray)})}^2{]}^{1/2}$, where ${\mathit{r}}_k$ and ${\mathit{r}}_k^{(X-ray)}$ are the positions of the $k$-th heavy atom of ribocil in the snapshot and the X-ray complex structure, respectively, and $N_{rib}$ the number of heavy atoms in ribocil ($N_{rib}=27$ for both ribocil A and B). Then, we collected snapshots satisfying $rms{d}_{rib}\le 5.0\mathring{\mathrm{A}}$ because we confirmed that ribocil was completely in the binding pocket when satisfying $rms{d}_{rib}\le 5.0\mathring{\mathrm{A}}$. Next, we assigned the stacking pairs to the collected snapshots using the method introduced in Supplementary Subsection 4.3.

We defined the numbers of native and non-native stackings in each snapshot, which are, respectively, denoted as $N_{native}^{\left(stack\right)}$ and $N_{non-native}^{\left(stack\right)}$ (see Supplementary Subsection 4.3 for details). Then, calculating $N_{native}^{\left(stack\right)}$ and $N_{non-native}^{\left(stack\right)}$ for all the collected snapshots, we computed the probability distribution assigned to $N_{native}^{\left(stack\right)}$ (i.e., $p_{\alpha }\left(N_{native}^{\left(stack\right)}\right)$) and that assigned to $N_{non-native}^{\left(stack\right)}$ (i.e., $p_{\alpha }\left(N_{non-native}^{\left(stack\right)}\right)$), where $\alpha$ is the system’s indicator ($\alpha =$ Ribo-A or Ribo-B), and the statistical weight $w_i$ assigned to each snapshot was used to calculate the probability (see Supplementary Subsection 4.1). Figure 7a shows that the distribution of $N_{native}^{\left(stack\right)}$ was similar between the two systems, although the probability of $N_{native}^{\left(stack\right)}=3$ for the Ribo-B system was slightly larger than that for the Ribi-A system: $p_{Ribo-B}(N_{native}^{\left(stack\right)}=3)/p_{Ribo-A}(N_{native}^{\left(stack\right)}=3)=1.8$. In contrast, the non-native stacking was formed more in the Ribo-B system than in Ribo-A (Figure 7b).

Figure 7 

(a) Probabilities $p_{Ribo-A}\left(N_{native}^{\left(stack\right)}\right)$ and $p_{Ribo-B}\left(N_{native}^{\left(stack\right)}\right)$. (b) Probabilities $p_{Ribo-A}\left(N_{non-native}^{\left(stack\right)}\right)$ and $p_{Ribo-B}\left(N_{non-native}^{\left(stack\right)}\right)$. See Supplementary Subsection 4.3 for definition of $N_{native}^{\left(stack\right)}$ and $N_{non-native}^{\left(stack\right)}$. (c) $<rms{d}_{rib}>$ as a function of $N_{native}^{\left(stack\right)}$ and (d) $<rms{d}_{rib}>$ as a function of $N_{non-native}^{\left(stack\right)}$. Snapshots with $N_{non-native}^{\left(stack\right)}=5$ were not observed for both Ribo-A and Ribo-B systems.

Next, we investigated the stacking further by calculating the average of $rms{d}_{rib}$ (i.e., $<rms{d}_{rib}>$) as a function of $N_{native}^{\left(stack\right)}$ or $N_{non-native}^{\left(stack\right)}$, where the statistical weight $w_i$ was also used. Figure 7c demonstrates that $<rms{d}_{rib}>$ increases with decreasing $N_{native}^{\left(stack\right)}$: The more the number of native stackings, the closer the conformation to the X-ray position. In contrast, Figure 7d shows that $<rms{d}_{rib}>$ increased with increasing $N_{non-native}^{\left(stack\right)}$.

Remember that the highest-density spot of ribocil B deviated by about $4\mathring{\mathrm{A}}$ from the X-ray position toward the front gate of the bonding pocket of the aptamer (red arrow in Figure 4d). Similarly, ribocil B deviated from the X-ray position when non-native stacking was formed (cyan arrows of Supplementary Figure S8d–f). In contrast, ribocil B with the native stacking did not deviate largely from the X-ray position (binding site) (Supplementary Figures S8a–c). These results propose a binding scenario as follows. Suppose that ribocil is entering the binding pocket from the front gate. The non-native stackings are formed first before ribocil reaching the binding site, and when ribocil binds to the binding site, the non-native stackings are replaced by the native stackings. We consider that this binding scenario occurs only when ribocil enters the binding pocket from the front gate. Probability for ribocil entering from the rear gate was less probable because of the free-energy barrier (Figure 5). The binding process from the rear gate is argued in the next subsection.

In general, many factors contribute to the conformational stability of a biological molecule in a complicated manner. We presume that the non-native stacking is one of the main factors for stabilizing the highest-density spot in the Ribo-B system as shown in Figure 4d. On the other hand, it is likely that the stability of ribocil B at the binding site was less contributed by the native stacking. That is, ribocil B at the X-ray position is more stabilized if the native stacking is corrected to be stronger. We argue this point again in the Discussion section.

Figure 5 showed visually that ${\rho }_{CM}\left({\mathit{r}}_{cube}\right)$ for ribocil B in the vicinity of the X-ray position was larger than that for ribocil A. To check this visual impression, we calculated the probability assigned to a small volume around the X-ray position by:

$Q_{\alpha }\left(r_{lim}\right)={\sum }_i^{}{w}_i{\delta }_i\left(r_{lim}\right)$, (2)

where ${\delta }_i\left(r_{lim}\right)=1$ when ${\mathit{r}}_{CM,i}$ is in a sphere (radius=$r_{lim}$) centered at the X-ray position and ${\delta }_i\left(r_{lim}\right)=0$ otherwise. Remember that the binding modes of ribocil A and B are almost identical [21]. The parameter $\alpha$ is an indicator of system: $\alpha =A$ for Ribo-A and $\alpha =B$ for Ribo-B. The resultant ratio is: $\frac{Q_B\left(1.0\mathring{\mathrm{A}}\right)}{Q_A\left(1.0\mathring{\mathrm{A}}\right)}=2.8$ and $\frac{Q_B\left(3.0\mathring{\mathrm{A}}\right)}{Q_A\left(3.0\mathring{\mathrm{A}}\right)}=3.1$. Therefore, the X-ray complex structure (PDB ID: 5c45) for the Ribo-B system is more stable than that for the Ribo-A system. This result agrees with the difference between $p_{Ribo-B}(N_{native}^{\left(stack\right)}=3)$ and $p_{Ribo-A}(N_{native}^{\left(stack\right)}=3)$ shown above.

Next, we analyze stacking between A48 and A49 of the aptamer domain. This stacking is formed only in the apo form (PDB ID: 6wjr and 6wjs) and disappears in the holo form (PDB ID: 5c45 and 5kx9) because A49 is replaced by Ring R4 of ribocil B in the holo form. Then, to check this property of the A48-A49 stacking from the simulation data, we calculated the ratio of the A48–A49 stacking for native-like complex conformations and unbound conformations: The native-like complex conformations and the unbound conformations are referred to as “B-state conformations” and “UB-state conformations”, respectively, in this subsection. We defined that B-state conformations satisfy $rms{d}_{rib}\le 3.0\mathring{\mathrm{A}}$ and that the UB-state conformations do $r_{min}\ge 20\mathring{\mathrm{A}}$, where $r_{min}$ is the minimum atomic distances from the aptamer domain to ribocil. The B-state and UB-state probabilities are referred to as $P_B$ and $P_{UB}$, respectively, and defined as:

$P_{\alpha }=\frac{{\sum }_i^{}{\delta }_{48-49,i}^{\left(\alpha \right)}{u}_i^{\left(\alpha \right)}{\omega }_i}{{\sum }_i^{}{u}_i^{\left(\alpha \right)}{\omega }_i}$, (3)

where $\alpha$ specifies the state ($\alpha =B$ or $UB$), and ${\delta }_{48-49,i}^{\left(\alpha \right)}$ works indicate the stacking, defined as: ${\delta }_{48-49,i}^{\left(\alpha \right)}=1$ when the A48–A49 stacking is formed in snapshot $i$ and 0 otherwise. The identification of stacking is from Supplementary Subsection 4.3. The function $u_i^{\left(\alpha \right)}$ works to restrict snapshots only in the $\alpha$-state, defined as: $u_i^{\left(\alpha \right)}=1$ when the snapshot $i$ belongs to the $\alpha$-state and 0 otherwise.

The resultant $P_B$ and $P_{UB}$ for the Ribo-A system were: $P_B=28\%$ and $P_{UB}=83\%$ ($P_B+P_{UB}>100\%$ because each $P_B$ and $P_{UB}$ were normalized in each the B and $UB$ states, respectively, in Equation 3). Those for the Ribo-B system were: $P_B=0\%$ and $P_{UB}=79\%$. Therefore, we conclude that the switching of stacking from A48-A49 to A48–R4 was reproduced in the simulation for both systems. One may consider that $P_B=28\%$ for the Ribo-A system may be large comparing to $P_B=0.0\%$ for the Ribo-B system. This might be a reason for the less stability of ribocil A in the binding pocket than that of ribocil B. However, as discussed above, it is more probable that the less stability assigned to ribocil A was induced by the less stability of the non-native stacking regarding the Ribo-A system (Figure 7b).

 Ribocil in the Rear Portion of the Density Tunnel

Figures 5 and 8a indicate that the ligand’s approach from the rear gate of the aptamer is less important for both ribocil A and ribocil B because of the free-energy barrier, which requires the large molecular orientation ${\mathit{e}}_{\uparrow }$ of ribocil in the tunnel. This approach makes the molecular binding slow. However, this result raises a question: Why does ${\mathit{e}}_{\uparrow }$ tend to be antiparallel to ${\mathit{e}}_{\uparrow }^{(X-ray)}$ in the rear portion of the tunnel? If ${\mathit{e}}_{\uparrow }$ is parallel to ${\mathit{e}}_{\uparrow }^{(X-ray)}$, then ribocil can reach the X-ray position without the large orientational rearrangement, and then the free-energy barrier vanishes probably. We consider this question in this subsection.

We picked snapshots in that ribocil was in the rear portion of the density tunnel with condition of ${SP}_{\uparrow }\le -0.9$ and ${SP}_{\leftarrow }\ge 0.5$, and obtained 218 conformations. The first inequality specifies the snapshots to be almost antiparallel to ${\mathit{e}}_{\uparrow }^{(X-ray)}$. Remember that ${\mathit{e}}_{\uparrow }$ of ribocil in the rear portion tended to be antiparallel to ${\mathit{e}}_{\uparrow }^{(X-ray)}$ (Figure 6). The second inequality imposes the snapshots so that the angle formed by ${\mathit{e}}_{\leftarrow }$ and ${\mathit{e}}_{\leftarrow }^{(X-ray)}$ is $60{}^{\circ }$ or smaller. We found that the majority (90 %) of them formed two hydrogen bonds between R3 of ribocil and G98 or A99 of the aptamer domain as well as between R1 and A102. We also found that intermolecular stacking was rarely formed between ribocil and the aptamer (data not shown). Supplementary Figure S12 illustrates a typical snapshot with those hydrogen bonds. We conclude that snapshots with ${\mathit{e}}_{\uparrow }$ to be antiparallel to ${\mathit{e}}_{\uparrow }^{(X-ray)}$ were stabilized by those hydrogen bonds. In other words, the free-energy barrier is a consequence of the hydrogen bonds when ribocil is in the rear portion of the density tunnel.

This inhibition mechanism for ribocil, which is in the rear portion, is interesting. We, however, could not argue only from the current simulation results if this mechanism is prepared evolutionally. We expect that the current simulation stimulates not only the computational field but also experimental or evolutionary biology.

 Physical Quantities in the Density Tunnel

Here we analyze properties of the ribocil in the density tunnel. For this purpose, we introduced a cylinder representing the density tunnel. Supplementary Figure S13 visualizes the cylinder. The cylinder axis was introduced so that it perforates the center of the density tunnel and divided into small bins (see Supplementary Section 6 for details). The position of each bin was specified by a parameter $x_{cy}$ defined along the cylinder axis. We analyzed only the snapshots whose ligand’s centroids were involved in the cylinder: A physical quantity of those snapshots was expressed as a function of $x_{cy}$.

Supplementary Equations S11 defines a free-energy (potential of mean force) profile, $G_{cy}\left(x_{cy}\right)$, of the system along the cylinder’s axis $x_{cy}$, and Figure 8a demonstrates $G_{cy}\left(x_{cy}\right)$ for both systems. The Ribo-B system (the red line) had a free-energy minimum at $x_{cy}=4\mathring{\mathrm{A}}$. This minimum corresponds to the high-density spot (red-contour region) of Figure 4d for the Ribo-B system. Also Figure 8a shows the existence of the free-energy barrier at $x_{cy}=-2\mathring{\mathrm{A}}$, which corresponds to the broken line of Figure 5b. The Ribo-A system (the blue line of Figure 8a) also had a free-energy minimum at $x_{cy}=8\mathring{\mathrm{A}}$. One may consider from Figure 8a that the Ribo-A system had a sharp high-density spot around the entrance. However, there was no such a high-density spot in Figure 4b. Therefore, this free-energy minimum 8b for the Ribo-A system was contributed by the broad distribution of smaller density sites (i.e., blue-colored contours in Figure 4b), which were summed at $x_{cy}=8\mathring{\mathrm{A}}$ in computing $G_{cy}\left(x_{cy}\right)$. The Ribo-A system also exhibited the free energy barrier at $x_{cy}=-4\mathring{\mathrm{A}}$, which corresponds to the broken line of Figure 5a.

Figure 8 

Profile of physical quantities in cylinder axis $x_{cy}$. Definition of the cylinder is given in Supplementary Section 6 and Supplementary Figure S13. Blue- and red-lines are, respectively, from Ribo-A and Ribo-B systems. (a) Profile of free energy (potential of mean force) $G_{cy}\left(x_{cy}\right)$. (b) Profile of $<rms{d}_{rib}\left(x_{cy}\right)>$. (c) Profile of $<rms{d}_{rib}^{\left(shift\right)}\left(x_{cy}\right)>$. (d) Profile of $<{Rg}_{3base}\left(x_{cy}\right)>$. Filled triangle labelled “apo” indicates averages of ${Rg}_{3base}\left(x_{cy}\right)$ for two experimental apo structures (PDB ID: 5c45 and 5kx9), and that labelled “holo” does average of two experimental holo form structures (PDB ID: 6wjr and 6wjs).

Figure 8b plots $<rms{d}_{rib}\left(x_{cy}\right)>$ with setting ${\Omega }_{cy,i}=rms{d}_{rib,i}$in Supplementary Equation S12, where $rms{d}_{rib,i}$ is $rms{d}_{rib}$ of snapshot $i$. Note that $rms{d}_{rib,i}$, defined in the Result section, was computed after superimposing the aptamer domain between snapshot $i$ and the initial conformation of the simulation (see Supplementary Subsection 4.1). The profile $<rms{d}_{rib}\left(x_{cy}\right)>$ increased monotonically with $x_{cy}$ departing from the origin ($x_{cy}=0$) for both systems. This increment is because $rms{d}_{rib}$ was contributed mainly by the ligand’s translation from the X-ray position to the snapshot’s position. Then, we removed this translation by shifting the ligand’s centroid of the snapshot $i$ to the X-ray position (PDB ID: 5c45). Here, we refer to the $rms{d}_{rib}$ after the shift as $rms{d}_{rib}^{\left(shift\right)}$. Figure 8c plots $<rms{d}_{rib}^{\left(shift\right)}\left(x_{cy}\right)>$, which suggests that the molecular orientation of ribocil A tended to be more similar with that of the native complex than that of ribocil B did in the density tunnel. This tendency was not clear when we compared the molecular orientational quantities ($<{\mathit{e}}_{\alpha }>$ and ${SP}_{\alpha }$; $\alpha =\uparrow$ or $\leftarrow$) between Figure 6 (the Ribo-B system) and Figure S11 (the Ribo-A system) in the density tunnel.

The result that $<rms{d}_{rib}^{\left(shift\right)}\left(x_{cy}\right)>$ of ribocil A was smaller than that of ribocil B does not necessarily mean that the Ribo-B system produced the native-like complex more than the Ribo-B system. In fact, the native-like complex was produced more from the Ribo-B system than from the Ribo-A system (Figure 4). As we discussed regarding Figures 7b and d, the more stability of Ribo-B than Ribo-A was probably brought by the non-native $\pi$-$\pi$ stacking formed in the binding pocket.

Now, we analyze the size of the binding pocket as the function of $x_{cy}$. Remember that A48, G62, and A85 contributed to the native stackings (Supplementary Subsection 4.3) and that the bases surround the ligand from different directions of the binding pocket as shown in Supplementary Figure S6. Thus, we calculated the radius of gyration, $R{g}_{3base}$, of a clump comprised of the three bases as an index to quantify the pocket size.

$R{g}_{3base}$ was calculated using heavy-atomic positions of the three bases, and the atomic masses for all heavy atoms was set to a constant because the size of the binding pocket is spatial expanse independent of the atomic masses. Figure 8d plots $<{Rg}_{3base}\left(x_{cy}\right)>$, which indicates that the binding pocket shrunk when the ligand was outside the binding pocket ($x_{cy}\ge 20\mathring{\mathrm{A}}$ or $x_{cy}\le -20\mathring{\mathrm{A}}$) for both systems. When the ligand entered the binding pocket, $<{Rg}_{3base}>$ swelled ($20\mathring{\mathrm{A}}<x_{cy}<8\mathring{\mathrm{A}}$ or $-20\mathring{\mathrm{A}}<x_{cy}<-5\mathring{\mathrm{A}}$). Then, the binding pocket shrunk again when the ligand approached the position ($x_{cy}\approx 0\mathring{\mathrm{A}}$) in the density tunnel. This behavior of $<{Rg}_{3base}>$, which was similar for both systems, is reasonable because the entrance of the binding pocket should open wide enough to pass the ligand through the tunnel. It is likely that the decrement of the binding pocket with the ligand reaching the binding site was resulted from increasing of the ligand–aptamer packing near the native complex conformation.

On the other hands, the computed $R{g}_{3base}$ provided inconsistency with experimental values: $R{g}_{3base}$ for the holo form (PDB ID: 6wjr and 6wjs) are 7.47 $\mathring{\mathrm{A}}$ and 7.44 $\mathring{\mathrm{A}}$, respectively: Average=7.46 $\mathring{\mathrm{A}}$. $R{g}_{3base}$ for the apo form (PDB ID: 5c45 and 5kx9) are 7.27 $\mathring{\mathrm{A}}$ and 7.29 $\mathring{\mathrm{A}}$, respectively: Average=7.28 $\mathring{\mathrm{A}}$, which are shown in Figure 8d. Those experimental values were computed by us with using the crystallographic structures. Remember that we used only heavy atoms to calculate $R{g}_{3base}$. Apparently, $R{g}_{3base}$ of apo form was larger than the computational value ($x_{cy}\sim 0\mathring{\mathrm{A}}$) by 0.25 $\mathring{\mathrm{A}}$. We consider the following two reasons for this inconsistency. First, the solution condition is different between the experiment and simulation. In the experimental structures of the apo form, the binding pocket involves an SO₄ ion (PDB ID: 6wjr) and a PO₄ ion (PDB ID: 6wjs). This means that the binding pocket was solvated in the apo state in the experiments. Contrarily, in the simulation, no water molecules or ions were involved in the binding pocket in the unbound conformations. This is likely to be the main reason for the shrinking of the binding pocket in the computation. The second reason, which may be an indirect reason, is the difference of the system’s composition: The experimental apo structure was obtained from a crystal, although a single aptamer was put in solution in computation. Because an RNA molecule is more flexible than a globular protein, this second reason may affect seriously to the inconsistency of $R{g}_{3base}$ between the experiment and computation. So far, we cannot specify which reason of the two is the most probable for the inconsistency. However, we consider that the computed behavior of $R{g}_{3base}$ as the function of $x_{cy}$ is useful to consider the binding process of this system because an experiment provides ${Rg}_{3base}$ only in the apo and holo forms: No experimental data is given for the intermediate state between the apo and holo forms.

Figure 6 and Supplementary Figure S10 showed that the ribocil’s molecular orientation (${\mathit{e}}_{\uparrow }$ and ${\mathit{e}}_{\leftarrow }$) was ordered well in the binding pocket. As discussed, this orientation ordering is important to lead both ribocil A and B to the binding position (the position in the X-ray complex structure) smoothly. We named this mechanism the orientation selection. To assess this molecular orientation ordering, a relatively compact structure (Supplementary Figure S5) should be formed. In other words, Rings R1 and R4 should be relatively close to each other. If ribocil adopts an extended conformation, then its shape is not compact although ${\mathit{e}}_{\uparrow }$ is computable formally for such an extended conformation.

Here, we denote the inter-centroid distance between Rings R1 and R4 as $r_{R1-R4}$. Then, to check the compactness of ribocil, we calculated two profiles regarding $r_{R1-R2}$ as a function of $x_{cy}$: $<r_{R1-R2}\left(x_{cy}\right)>$ and ${SD}_{r_{R1-R2}}\left(x_{cy}\right)$. The latter is the amount of fluctuation of $r_{R1-R2}$ defined as: ${SD}_{r_{R1-R2}}\left(x_{cy}\right)=<r_{R1-R2}{\left(x_{cy}\right)}^2>-<r_{R1-R2}\left(x_{cy}\right){>}^2$. Figure. 8e demonstrates that ribocil B became compact when ribocil B entered the binding pocket. On the other hand, ribocil A was always compact. Because ${SD}_{r_{R1-R2}}\left(x_{cy}\right)$ was relatively small for both systems, the shape of ribocil was compact in the binding pocket for both systems. Because the compactness of ribocil B was induced when entering the binding pocket, the induced-fit mechanism exerted in the binding process for the Ribo-B system.

 Contact of a Hydrogen Atom of Ribocil to A48 and A85 in the Native-Like Complex Structure

Ribocil A and B differ mutually in the position of a hydrogen, which is shown by the red-colored hydrogen atom in Figures 1a and b. We refer to this atom as “$H_{diff}$ atom” for convenience. It is likely that the $H_{diff}$ atom tends to contact A48 for ribocil A, and that it tends to do A85 for ribocil B. To check this conjecture, we calculated radial distribution functions: ${\rho }_{rd}\left(r_{H_{diff}-A48}^{\left(min\right)}\right)$ and ${\rho }_{rd}\left(r_{H_{diff}-A85}^{\left(min\right)}\right)$, where $r_{H-A48}^{\left(min\right)}$ and $r_{H-A85}^{\left(min\right)}$ are the minimum distances from the $H_{diff}$ atom to the sidechain heavy atoms of A48 and A85, respectively. See Supplementary Subsection 7 for details of the radial distribution functions. In fact, Figure 9 demonstrates that the $H_{diff}$ atom of ribocil A tended to contact to A48 and that the $H_{diff}$ atom of ribocil B did to contact to A85.

Figure 9 

Radial distribution functions ${\rho }_H\left(r_{H-A48}^{\left(min\right)}\right)$ and ${\rho }_H\left(r_{H-A85}^{\left(min\right)}\right)$ of two distances $r_{H-A48}^{\left(min\right)}$ and $r_{H-A85}^{\left(min\right)}$, respectively. The red-colored hydrogen atom in Figure 1a and b, referred to as “$H_{diff}$ atom” in this paper, differs between ribocil A and B.

We also checked if the position of the $H_{diff}$ atom affected two native stackings A48–R4 and A85–R3. For this analysis, we used snapshots with $rms{d}_{rib}\le 2.5\mathring{\mathrm{A}}$ for both systems, which means that only native-like complex conformations were used. Then, two ratios of the stacking formation, $ratio(A48-R4)$ and $ratio\left(A85–R3\right)$, were calculated: The former is the ratio of formation of the A48-R4 stacking and the latter is that for A85–R3. The resultant ratios were: For ribocil A, $ratio(A48-R4)=16\%$ and $ratio\left(A85–R3\right)=44\%$. For ribocil B, $ratio(A48-R4)=21\%$ and $ratio\left(A85–R3\right)=77\%$. These results indicates that ribocil B involved more the native stackings than ribocil A did, which suggests that ribocil B binds to the aptamer domain more strongly to the aptamer domain than ribocil A does. However, the difference of the $H_{diff}$-atom positioning in ribocil did not vary the magnitude relation between $ratio(A48-R4)$ and $ratio(A48-R4)$. I.e., the inequality of $ratio(A48-R4)<ratio(A85-R3)$ was seen for both systems.

 Discussion

Existence of encounter complex has been proposed in protein–ligand or protein–protein binding [51–54]. The present study has suggested that the ligand experiences encounter complexes before reaching the binding site of the aptamer (Figure 3). Starting from the unbound state, the ligand could contact the whole surface of riboswitch (green contours of Figure 4). Note that the green-contour region has a lower free energy or potential of mean force ($PMF$) than the unbound region outside the green-contour region because $PMF$ is defined by: $PMF=-RTln\left[{\rho }_{CM}\left({\mathit{r}}_{cube}\right)\right]$, where $R$ is gas constant. The ligand found some stable sites (magenta contour sites) when moving on the aptamer’s surface. The density ${\rho }_{CM}\left({\mathit{r}}_{cube}\right)$ increased further in the vicinity of the front gate (blue contours in Figures 4b and d), where the entrance of the binding pocket opened with the ligand entering the front or rear gate ($20\mathring{\mathrm{A}}<x_{cy}<8\mathring{\mathrm{A}}$ or $-5\mathring{\mathrm{A}}>x_{cy}<-5\mathring{\mathrm{A}}$ in Figure 8d). Last, the ligand was encapsulated in the binding pocket and the density tunnel shrunk when the ligand moved to the binding position for both systems ($x_{cy}\to 0$ in Figure 8d).

The density of ribocil B in the binding pocket was higher than that assigned to ribocil A (Figure 4). This indicates that the aptamer–ribocil B complex was more stable than the aptamer–ribocil A complex, which agrees qualitatively with the experimental result [21]. However, Ref. 21 reported that the difference of the binding free-energy between the Ribo-A and ribo-B systems was 4.3 kcal/mol. This means that in the experiment a probability assigned to the binding state for the Ribo-B system is about ${10}^3$ times larger than that for the Ribo-A system, and that our computation underestimated the free-energy difference. Remember that the Results section showed that the native-stacking interactions are counted insufficiently for the Ribo-B system. We presume that this insufficiency is the reason for the underestimation of the free-energy difference.

We propose two possibilities for this insufficiency: One is an inaccuracy of the force field, and the other is inefficiency of sampling. It has been reported that the currently used force field χ_OL3 is appropriate for treating an RNA system by molecular simulation [55]. However, a quantum-mechanical computation has shown that the potential-energy surface for the $\pi$-$\pi$ stacking involves multiple energy minima and is sensitive to the type of interacting base rings [56]. This means that the $\pi$-$\pi$ stacking interaction is a difficult energy term to be approximated classically for molecular dynamics simulations.

We discuss the second possibility, sampling inefficiency. Figures 8a–c and 9a–c showed that ribocil B with the native-stackings was superimposed well to that of the X-ray complex, which means that the native-like complex structure was searched during the simulation. Therefore, it is likely that the probability (or stability) assigned to the native-like complex increases if the stacking interaction is improved. Therefore, we consider that the inaccuracy regarding the stacking interaction is a major cause for the instability of the native stackings.

The ligand orientational ordering in the binding pocket (Figure 6 and Supplementary Figure S10) found in the present study has a similarity with the binding mechanism found in another computational study: Binding of a GPCR (human endothelin receptor type B) and a drug molecule (bosentan) [30]. Namely, the molecular orientation of bosentan ordered to enter the deep binding pocket of GPCR. We presume that the molecular orientational ordering is important when a ligand binds to the deep pocket because the ligand (not a small compound consisting of a few atoms) cannot change the orientation readily in the deep pocket. Therefore, this orientational ordering works as a selection mechanism (“orientational selection mechanism”) for the ligand–receptor binding. This mechanism can be regarded as a type of the conformation selection mechanism [57]. The difference between two studies is: In the present study, the orientational ordering was supported by the ligand–receptor $\pi$-$\pi$ stacking, whereas the ordering in the earlier study [30] was done by the interaction between the ligand and the long-disordered N-terminal tail of GPCR.

Both the apo and holo forms of the aptamer domain (PDB ID: 6wjr [20] and 5c45 [19], respectively) have a tunnel, which excavates between the front and rear gates of the aptamer (Supplementary Figure S2). In fact, the special density ${\rho }_{CM}\left({\mathit{r}}_{cube}\right)$ showed that both ribocil A and B could pass through the tunnel (Figures 4b and 4d), and this tunnel was denoted as the density tunnel in the present study. The density tunnel was divided into the front and rear portions, and the front portion involved the binding pocket (Figure 6a). Interestingly, ${\rho }_{CM}\left({\mathit{r}}_{cube}\right)$ showed a free-energy barrier at the boundary of the front and rear portions (Figure 5b and Figure 8a). Further analysis showed that the molecular orientation $<{\mathit{e}}_{\uparrow }>$ in the rear portion tended to be opposite to that in the front portion (Figure 6a) and ${\mathit{e}}_{\uparrow }^{(X-ray)}$. This means that the molecular orientation should turn completely and sharply at the boundary when the ligand passes the boundary from the rear to front portion (Supplementary Figure S11). Therefore, we presume that the main binding pathway is the approach from the front gate to the binding pocket.

The mD-VcMD method was performed to obtain the overall free-energy landscape, by which binding process can be discussed. For this purpose, sampling should be performed widely in the conformational space. Therefore, the mD-VcMD method belongs to the “global sampling” [24]. On the other hand, the detailed free-energy differences among free-energy basins and the heights of free-energy barriers in the conformational space may not be estimated accurately when the system is large or when the simulation length is short. Here, we denote the volume to be sampled by a global sampling method $V_{gloval}$.

If $PMF$ is calculated along a pathway (1D line or a narrow tube around the 1D line) in the 3D real space, the volume to be sampled, $V_{local}$, decreases drastically comparing with $V_{gloval}$: $V_{local}\ll V_{gloval}$. This sampling method belongs to “local sampling” [24]. Denoting the starting and ending points of the pathway as ${\mathit{r}}_{st}$ and ${\mathit{r}}_{en}$, respectively, the difference of $PMF$ between these two points are given by $PMF\left({\mathit{r}}_{en}\right)-PMF\left({\mathit{r}}_{st}\right)$ independent of the shape of the pathway. By setting ${\mathit{r}}_{st}$ and ${\mathit{r}}_{en}$ to an unbound conformation and the native complex structure, respectively, one can compute the binding free energy after some compensation to the obtained values of $PMF\left({\mathit{r}}_{en}\right)$ and $PMF\left({\mathit{r}}_{st}\right)$, which we do not explain in this paper.

We presume that in general the local sampling yields a binding free energy more accurately than the global sampling does because of $V_{local}\ll V_{gloval}$. Therefore, there has been many methods readily applicable [58–64] or having been applied [65–73] to the local sampling. The global sampling and local sampling should be used depending on the purpose of the study. When one investigates the overall free-energy landscape and thermodynamically probable binding process (or scenario), the global sampling is appropriate. The overall landscape is not obtained from the local sampling because its sampling is restricted in the small volume along the pathway. On the other hand, because the small volume can be sampled readily in general, the local sampling can calculate the free-energy difference between two conformations more accurately than the global sampling does.

Last, we discuss a positive role of the non-native stacking in the molecular binding. In general, the stacking is formed when two rings have appropriate positioning mutually, and it is likely that the non-native stacking is temporally formed before ribocil B reaches the binding site (i.e., X-ray position), as we showed in this paper. Remember that U61 of the aptamer domain is located near the gate of the binding pocket, and then ribocil B may be trapped shortly by stacking with U61 as shown in Supplementary Figures S8d–f. These figures demonstrate that the non-native stacking between U61 and R4 of ribocil supports the molecular orientation ${\mathit{e}}_{\uparrow }$ to be parallel to ${\mathit{e}}_{\uparrow }^{(X-ry)}$. Then, riboil can reach the X-ray position by translational motions without rotation. In this scenario, the positioning of U61 is important to arrange the ribocil’s orientation to be advantageous for the molecular binding. Last, we emphasize that this scenario is less influenced by the accuracy of p-p stacking force field. Therefore, the present study is useful to understand the RNA-ligand binding process.

We uploaded important snapshots obtained from the current simulation to the Biological Structure Model Archive (BSM-Arc) [74].

 Conflict of Interest

All authors declare that they have no conflict of interest.

 Author Contributions

J.H., N.K., I.F., and Y.F. designed the research plan; J.H. developed the mD-VcMD algorithm; G.J.B. and N.K. generated the molecular system for simulation and input files for gromacs, and then performed preparatory simulations; J.H. performed the main part of the simulation and analyses; I.F. and Y.F. performed some important parts of analyses; J.H. wrote the paper.

 Data Availability

We uploaded important snapshots obtained from the current simulation to the Biological Structure Model Archive (BSM-Arc). The entry for this study is https://bsma.pdbj.org/entry/49.

A preliminary version of this work, DOI: https://biorxiv.org/cgi/content/short/2023.07.01.547313v1, was deposited in the bioRxiv on July 03, 2023.

 Acknowledgements

We are grateful to Dr. Kota Kasahara from Ritsumeikan University for implementation of the mD-VcMD algorithm to gromacs using the PLUMED plug-in. This work was supported by JSPS KAKENHI Grants No. 21K06052 (J.H.) and 20H03229 (N.K.), and performed in part under the Cooperative Research Program of the Institute for Protein Research, Osaka University, CR22-02 and CR-23-02. J.H. and K.N. acknowledge support by the HPCI System Research Project (Project IDs: hp220002, hp220015, hp220022, hp230003, and hp230011). J.H., I.F., N.K., and Y.F. are supported by the Project Focused on Developing Key Technology for Discovering and Manufacturing Drugs for Next-Generation Treatment and Diagnosis (project ID: 23ae0121028h0003; 2018–2021 and 2021–) from AMED and the Japan Biological Informatics Consortium (JBiC). GA-mD-VcMD was performed on TSUBAME3.0 supercomputers at the Tokyo Institute of Technology. We used UCSF Chimera ver. 1.15 for drawing molecular structures and gnuplot for graphs. Other analyses were done using homemade programs.

References

• [1]   Bernetti,  M.,  Aguti,  R.,  Bosio,  S.,  Recanatini,  M.,  Masetti,  M.,  Cavalli,  A. Computational drug discovery under RNA times. QRB Discovery 3, E22 (2022). https://doi.org/10.1017/qrd.2022.20
• [2]   Manigrasso,  J.,  Marcia,  M.,  De Vivo,  M. Computer-aided design of RNA-targeted small molecules: A growing need in drug discovery. Chem 7, 2965–2988 (2021). https://doi.org/10.1016/j.chempr.2021.05.021
• [3]   Morishita,  E. C. Discovery of RNA-targeted small molecules through the merging of experimental and computational technologies. Expert Opin. Drug Discov. 18, 207–226 (2023). https://doi.org/10.1080/17460441.2022.2134852
• [4]   Childs-Disney,  J. L.,  Yang,  X.,  Gibaut,  Q. M. R.,  Tong,  Y.,  Batey,  R. T.,  Disney,  M. D. Targeting RNA structures with small molecules. Nat. Rev. Drug Discov. 21, 736–762 (2022). https://doi.org/10.1038/s41573-022-00521-4
• [5]   Bagnolini,  G.,  Luu,  T. B.,  Hargrove,  A. E. Recognizing the power of machine learning and other computational methods to accelerate progress in small molecule targeting of RNA. RNA 29, 473–488 (2023). https://doi.org/10.1261/rna.079497.122
• [6]   Kognole,  A. A.,  Hazel,  A.,  MacKerell,  A. D. Jr. SILCS-RNA: Toward a structure-based drug design approach for targeting RNAs with small molecules. J. Chem. Theory Comput. 18, 5672–5691 (2022). https://doi.org/10.1021/acs.jctc.2c00381
• [7]   Nahvi,  A.,  Sudarsan,  N.,  Ebert,  M. S.,  Zou,  X.,  Brown,  K. L.,  Breaker,  R. R. Genetic control by a metabolite binding mRNA. Chem. Biol. 9, 1043–1049 (2002). https://doi.org/10.1016/s1074-5521(02)00224-7
• [8]   Mironov,  A. S.,  Gusarov,  I.,  Rafikov,  R.,  Lopez,  L. E.,  Shatalin,  K.,  Kreneva,  R. A., et al. Sensing small molecules by nascent RNA: A mechanism to control transcription in bacteria. Cell 111, 747–756 (2002). https://doi.org/10.1016/s0092-8674(02)01134-0
• [9]   Winkler,  W.,  Nahvi,  A.,  Breaker,  R. R. Thiamine derivatives bind messenger RNAs directly to regulate bacterial gene expression. Nature 419, 952–956 (2002). https://doi.org/10.1038/nature01145
• [10]   Winkler,  W. C.,  Cohen-Chalamish,  S.,  Breaker,  R. R. An mRNA structure that controls gene expression by binding FMN. Proc. Natl. Acad. Sci. U.S.A. 99, 15908–15913 (2002). https://doi.org/10.1073/pnas.212628899
• [11]   Serganov,  A.,  Nudler,  E. A decade of riboswitches. Cell 152, 17–24 (2013). https://doi.org/10.1016/j.cell.2012.12.024
• [12]   Warner,  K. D.,  Hajdin,  C. E.,  Weeks,  K. M. Principles for targeting RNA with drug-like small molecules. Nat. Rev. Drug Discov. 17, 547–558 (2018). https://doi.org/10.1038/nrd.2018.93
• [13]   Gelfand,  M. S.,  Mironov,  A. A.,  Jomantas,  J.,  Kozlov,  Y. I.,  Perumov,  D. A. A conserved RNA structure element involved in the regulation of bacterial riboflavin synthesis genes. Trends Genet. 15, 439–442 (1999). https://doi.org/10.1016/s0168-9525(99)01856-9
• [14]   Vitreschak,  A. G.,  Rodionov,  D. A.,  Mironov,  A. A.,  Gelfand,  M. S. Regulation of RF biosynthesis and transport genes in bacteria by transcriptional and translational attenuation. Nucleic Acids Res. 30, 3141–3151 (2002). https://doi.org/10.1093/nar/gkf433
• [15]   Vicens,  Q.,  Mondragón,  E.,  Batey  R. T. Molecular sensing by the aptamer domain of the FMN riboswitch: A general model for ligand binding by conformational selection. Nucleic Acids Res. 39, 8586–8598 (2011). https://doi.org/10.1093/nar/gkr565
• [16]   Serganov,  A.,  Huang,  L.,  Patel,  D. J. Coenzyme recognition and gene regulation by a flavin mononucleotide riboswitch. Nature 458, 233–237 (2009). https://doi.org/10.1038/nature07642
• [17]   Rizvi,  N. F.,  Howe,  J. A.,  Nahvi,  A.,  Klein,  D. J.,  Fischmann,  T. O.,  Kim,  H.-Y., et al. Discovery of selective RNA-binding small molecules by affinity-selection mass spectrometry. ACS Chem. Biol. 13, 820–831 (2018). https://doi.org/10.1021/acschembio.7b01013
• [18]   Vicens,  Q.,  Mondragón,  E.,  Reyes,  F. E.,  Coish,  P.,  Aristoff,  P.,  Berman,  J., et al. Structure-activity relationship of flavin analogues that target the flavin mononucleotide riboswitch. ACS Chem. Biol. 13, 2908–2919 (2018). https://doi.org/10.1021/acschembio.8b00533
• [19]   Howe,  J. A.,  Wang,  H.,  Fischmann,  T. O.,  Balibar,  C. J.,  Xiao,  L.,  Galgoci,  A. M., et al. Selective small-molecule inhibition of an RNA structural element. Nature 526, 672–677 (2015). https://doi.org/10.1038/nature15542
• [20]   Wilt,  H. M.,  Yu,  P.,  Tan,  K.,  Wang,  Y. X.,  Stagno,  J. R. FMN riboswitch aptamer symmetry facilitates conformational switching through mutually exclusive coaxial stacking configurations. J. Struct. Biol. X 4, 100035 (2020). https://doi.org/10.1016/j.yjsbx.2020.100035
• [21]   Howe,  J. A.,  Xiao,  L.,  Fischmann,  T. O.,  Wang,  H.,  Tang,  H.,  Villafania,  A., et al. Atomic resolution mechanistic studies of ribocil: A highly selective unnatural ligand mimic of the E. coli FMN riboswitch. RNA Biol. 13, 946–954 (2016). https://doi.org/10.1080/15476286.2016.1216304
• [22]   Higo,  J.,  Ikebe,  J.,  Kamiya,  N.,  Nakamura,  H. Enhanced and effective conformational sampling of protein molecular systems for their free energy landscapes. Biophys. Rev. 4, 27–44 (2012). https://doi.org/10.1007/s12551-011-0063-6
• [23]   Kasahara,  K.,  Terazawa,  H.,  Takahashi,  T.,  Higo,  J. Studies on molecular dynamics of intrinsically disordered proteins and their fuzzy complexes: a mini-review. Comput. Struct. Biotechnol. J. 17, 712–720 (2019). https://doi.org/10.1016/j.csbj.2019.06.009
• [24]   Fukunishi,  Y.,  Higo,  J.,  Kasahara,  K. Computer simulation of molecular recognition in biomolecular system: From in-silico screening to generalized ensembles. Biophys. Rev. 14, 1423–1447 (2022). https://doi.org/10.1007/s12551-022-01015-8
• [25]   Higo,  J.,  Kusaka,  A.,  Kasahara,  K.,  Kamiya,  N.,  Hayato,  I.,  Qilin,  X., et al. GA-guided mD-VcMD: A genetic-algorithm-guided method for multi-dimensional virtual-system coupled molecular dynamics. Biophys. Physicobiol. 17, 161–176 (2020). https://doi.org/10.2142/biophysico.BSJ-2020008
• [26]   Higo,  J.,  Kawabata,  T.,  Kusaka,  A.,  Kasahara,  K.,  Kamiya,  N.,  Fukuda,  I., et al. Molecular interaction mechanism of a 14-3-3 protein with a phosphorylated peptide elucidated by enhanced conformational sampling. J. Chem. Inf. Model. 60, 4867–4880 (2020). https://doi.org/10.1021/acs.jcim.0c00551
• [27]   Higo,  J.,  Takashima,  H.,  Fukunishi,  Y.,  Yoshimori,  A. Generalized-ensemble method study: A helix-mimetic compound inhibits protein–protein interaction by long-range and short-range intermolecular interactions. J. Comput. Chem. 42, 956–969 (2021). https://doi.org/10.1002/jcc.26516
• [28]  Hayami, T., Kamiya, N., Kasahara, K., Kawabata, T., Kurita, J.-I., Fukunishi, Y. Difference of binding modes among three ligands to a receptor mSin3B corresponding to their inhibitory activities. Sci. Rep. 11, 6178 (2021). https://doi.org/10.1038/s41598-021-85612-9
• [29]   Higo,  J.,  Kasahara,  K.,  Wada,  M.,  Dasgupta,  B.,  Kamiya,  N.,  Hayami,  T., et al. Free-energy landscape of molecular interactions between endothelin 1 and human endothelin type B receptor: Fly-casting mechanism. Protein Eng. Des. Sel. 32, 297–308 (2019). https://doi.org/10.1093/protein/gzz029
• [30]   Higo,  J.,  Kasahara,  K.,  Bekker,  G.-J.,  Ma,  M.,  Sakuraba,  S.,  Iida,  S., et al. Fly casting with ligand sliding and orientational selection supporting complex formation of a GPCR and a middle sized flexible molecule. Sci. Rep. 12, 13792 (2022). https://doi.org/10.1038/s41598-022-17920-7
• [31]  Xie, Q., Kasahara, K., Higo, J., Takahashi, T. Molecular mechanisms of functional modulation of transcriptional coactivator PC4 via phosphorylation on its intrinsically disordered region. ACS Omaga 8, 14572−14582 (2023). http://doi.org/10.1021/acsomega.3c00364
• [32]   Zgarbová,  M.,  Otyepka,  M.,  Šponer,  J.,  Mládek,  A.,  Banáš,  P.,  Cheatham  T. E. III, et al. Nucleic acids force field based on reference quantum chemical calculations of glycosidic torsion profiles. J. Chem. Theory Comput. 7, 2886–2902 (2011). https://doi.org/10.1021/ct200162x
• [33]   Izadi,  S.,  Onufriev,  A. V. Accuracy limit of rigid 3-point water models. J. Chem. Phys. 145, 074501 (2016). https://doi.org/10.1063/1.4960175
• [34]   Li,  Z.,  Song,  L. F.,  Li  P.,  Merz,  K. M. Jr. Systematic parametrization of divalent metal ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB water models. J. Chem. Theory Comput. 16, 4429–4442 (2020). https://doi.org/10.1021/acs.jctc.0c00194
• [35]   Sengupta,  A.,  Li,  Z.,  Song,  L. F.,  Li,  P.,  Merz,  K. M. Jr. Parameterization of monovalent ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB water models. J. Chem. Inf. Model. 61, 869–880 (2021). https://doi.org/10.1021/acs.jcim.0c01390
• [36]  Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., et al. Gaussian 09, Revision D.01. Gaussian, Inc., Wallingford CT. (2009).
• [37]   Bayly,  C. I.,  Cieplak,  P.,  Cornell,  W.,  Kollman,  P. A. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: The RESP model. J. Phys. Chem. 97, 10269–10280 (1993). https://doi.org/10.1021/j100142a004
• [38]   Wang,  J. M.,  Wolf,  R. M.,  Caldwell,  J. W.,  Kollman,  P. A.,  Case,  D. A. Development and testing of a general amber force field. J. Comput. Chem. 25, 1157–1174 (2004). https://doi.org/10.1002/jcc.20035
• [39]   Higo,  J.,  Kamiya,  N.,  Sugihara,  T.,  Yonezawa,  Y.,  Nakamura,  H. Verifying trivial parallelization of multicanonical molecular dynamics for conformational sampling of a polypeptide in explicit water. Chem. Phys. Lett. 473, 326–329 (2009). https://doi.org/10.1016/j.cplett.2009.03.077
• [40]   Ikebe,  J.,  Umezawa,  K.,  Kamiya,  N.,  Sugihara,  T.,  Yonezawa,  Y.,  Takano,  Y., et al. Theory for trivial trajectory parallelization of multicanonical molecular dynamics and application to a polypeptide in water. J. Comput. Chem. 32, 1286–1297 (2011). https://doi.org/10.1002/jcc.21710
• [41]   Kutzner,  C.,  Páll,  S.,  Fechner,  M.,  Esztermann,  A.,  de Groot  B. L.,  Grubmüller,  H. More bang for your buck: Improved use of GPU nodes for GROMACS 2018. J. Comput. Chem. 40, 2418–2431 (2019). https://doi.org/10.1002/jcc.26011
• [42]   Bonomi,  M.,  Branduardi,  D.,  Bussi,  G.,  Camilloni,  C.,  Provasi,  D.,  Raiteri,  P., et al. PLUMED: A portable plugin for free-energy calculations with molecular dynamics. Comput. Phys. Commu. 180, 1961–1972 (2009). https://doi.org/10.1016/j.cpc.2009.05.011
• [43]   Hess,  B.,  Bekker,  H.,  Berendsen,  H. J. C.,  Fraaije,  J. G. E. M. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 18, 1463–1472 (1997). https://doi.org/10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0.CO;2-H
• [44]   Nosé,  S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 81, 511–519 (1984). https://doi.org/10.1063/1.447334
• [45]   Hoover,  W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 31, 1695–1697 (1985). https://doi.org/10.1103/PhysRevA.31.1695
• [46]   Kamiya,  N.,  Fukuda,  I.,  Nakamura,  H. Application of zero-dipole summation method to molecular dynamics simulations of a membrane protein system. Chem. Phys. Lett. 568–569, 26–32 (2013). https://doi.org/10.1016/j.cplett.2013.03.014
• [47]   Fukuda,  I.,  Kamiya,  N.,  Yonezawa,  Y.,  Nakamura,  H. Simple and accurate scheme to compute electrostatic interaction: Zero-dipole summation technique for molecular system and application to bulk water. J. Chem. Phys. 137, 054314 (2012). https://doi.org/10.1063/1.4739789
• [48]   Fukuda,  I.,  Yonezawa,  Y.,  Nakamura,  H. Molecular dynamics scheme for precise estimation of electrostatic interaction via zero-dipole summation principle. J. Chem. Phys. 134, 164107 (2011). https://doi.org/10.1063/1.3582791
• [49]   Kasahara,  K.,  Sakuraba,  S.,  Fukuda,  I. Enhanced sampling of molecular dynamics simulations of a polyalanine octapeptide: Effects of the periodic boundary conditions on peptide conformation. J. Phys. Chem. B 122, 2495–2503 (2018). https://doi.org/10.1021/acs.jpcb.7b10830
• [50]   Sakuraba,  S.,  Fukuda,  I. Performance evaluation of the zero-multipole summation method in modern molecular dynamics software. J. Comput. Chem. 39, 1551–1560 (2018). https://doi.org/10.1002/jcc.25228
• [51]   Crowley,  P. B.,  Rabe,  K. S.,  Worrall,  J. A. R.,  Canters,  G. W.,  Ubbink,  M. The ternary complex of cytochrome f and cytochrome c: identification of a second binding site and competition for plastocyanin binding. ChemBioChem 3, 526–533 (2002). https://doi.org/10.1002/1439-7633(20020603)3:6<526::AID-CBIC526>3.0.CO;2-N
• [52]   Xu,  X.,  Reinle,  W. G.,  Hannemann,  F.,  Konarev,  P. V.,  Svergun,  D. I.,  Bernhardt,  R., et al. Dynamics in a pure encounter complex of two proteins studied by solution scattering and paramagnetic NMR spectroscopy. J. Am. Chem. Soc. 130, 6395–6403 (2008). https://doi.org/10.1021/ja7101357
• [53]   Kozakov,  D.,  Li,  K.,  Hall,  D. R.,  Beglov,  D.,  Zheng,  J.,  Vakili,  P., et al. Encounter complexes and dimensionality reduction in protein–protein association. eLife 3, e01370 (2014). https://doi.org/10.7554/elife.01370
• [54]   Schilder,  J.,  Ubbink,  M. Formation of transient protein complexes. Curr. Opin. Struct. Biol. 23, 911–918 (2013). https://doi.org/10.1016/j.sbi.2013.07.009
• [55]   Šponer,  J.,  Bussi,  G.,  Krepl,  M.,  Banáš,  P.,  Bottaro,  S.,  Cunha,  R. A., et al. RNA structural dynamics as captured by molecular simulations: A comprehensive overview. Chem. Rev. 118, 4177–4338 (2018). https://doi.org/10.1021/acs.chemrev.7b00427
• [56]   van Mourik,  T.,  Hogan,  S. W. L. DNA base stacking involving adenine and 2-aminopurine. Struct. Chem. 27, 145–158 (2016). https://doi.org/10.1007/s11224-015-0708-3
• [57]  Spolar, R. S., Record, M. T. Jr. Coupling of local folding to site-specific binding of proteins to DNA. Science 263, 777–784 (1994). http://www.jstor.org/stable/2882917
• [58]   Fred,  G. Future paths for integer programming and links to artificial intelligence. Comput. Oper. Res. 13, 533–549 (1986). https://doi.org/10.1016/0305-0548(86)90048-1
• [59]   Beveridge,  D. L.,  DiCapua,  F. M. Free energy via molecular simulation: Applications to chemical and biomolecular systems. Annu. Rev. Biophys. Biophys. Chem. 18, 431–492 (1989). https://doi.org/10.1146/annurev.bb.18.060189.002243
• [60]   Huber,  T.,  Torda,  A. E.,  van Gunsteren,  W. F. Local elevation: A method for improving the searching properties of molecular dynamics simulation. J Comput. Aided Mol. Des. 8, 695–708 (1994). https://doi.org/10.1007/BF00124016
• [61]   Grubmüller,  H. Predicting slow structural transitions in macromolecular systems: Conformational flooding. Phys. Rev. E 52, 2893–2906 (1995). https://doi.org/10.1103/PhysRevE.52.2893
• [62]   Wang,  F.,  Landau,  D. P. Efficient, multiple-range random walk algorithm to calculate the density of states. Phys. Rev. Lett. 86, 2050–2053 (2001). https://doi.org/10.1103/PhysRevLett.86.2050
• [63]   Laio,  A.,  Parrinello,  M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 99, 12562–12566 (2002). https://doi.org/10.1073/pnas.202427399
• [64]   Hamelberg,  D.,  Mongan,  J.,  McCammon,  J. A. Accelerated molecular dynamics: a promising and efficient simulation method for biomolecules. J. Chem. Phys. 120, 11919–11929 (2004). https://doi.org/10.1063/1.1755656
• [65]   Kumar,  S.,  Rosenberg,  J. M.,  Bouzida,  D.,  Swendsen,  R. H.,  Kollman  P. A. THE weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. J. Comput. Chem. 13, 1011–1021 (1992). https://doi.org/10.1002/jcc.540130812
• [66]   Gelman,  A.,  Meng,  X.-L. Simulating normalizing constants: From importance sampling to bridge sampling to path sampling. Statist. Sci. 13, 163–185 (1998). https://doi.org/10.1214/ss/1028905934
• [67]   Fukunishi,  Y.,  Mikami,  Y.,  Nakamura,  H. The filling potential method: A method for estimating the free energy surface for protein-ligand docking. J. Phys. Chem. B 107, 13201–13210 (2003). https://doi.org/10.1021/jp035478e
• [68]   Fukunishi,  Y. Structure-based drug screening and ligand-based drug screening with machine learning. Comb. Chem. High Throughput Screen. 12, 397–408 (2009). https://dx.doi.org/10.2174/138620709788167890
• [69]   de Ruiter,  A.,  Oostenbrink,  C. Protein-ligand binding from distancefield distances and hamiltonian replica exchange simulations. J. Chem. Theory Comput. 9, 883–892 (2013). https://doi.org/10.1021/ct300967a
• [70]   Lier,  B.,  Öhlknecht,  C.,  de Ruiter,  A.,  Gebhardt,  J.,  van Gunsteren,  W. F.,  Oostenbrink,  C. A suite of advanced tutorials for the GROMOS biomolecular simulation software [Article v1.0]. Living J. Comp. Mol. Sci. 2, 18552 (2020). https://doi.org/10.33011/livecoms.2.1.18552
• [71]  Bekker, G.-J., Kamiya, N., Araki, M., Fukuda, I., Okuno, Y., Nakamura, H. Accurate prediction of complex structure and affinity for a flexible protein receptor and its inhibitor. J. Chem. Theory. Comput. 13, 2389−2399 (2017). https://doi.org/10.1021/acs.jctc.6b01127
• [72]  Bekker, G.-J., Araki, M., Oshima, K., Okuno, Y., Kamiya, N. Dynamic docking of a medium-sized molecule to its receptor by multicanonical MD simulations. J. Phys. Chem. B 123, 2479−2490 (2019). https://doi.org/10.1021/acs.jpcb.8b12419
• [73]   Bekker,  G.-J.,  Fukuda,  I.,  Higo,  J.,  Kamiya,  N. Mutual population-shift driven antibody-peptide binding elucidated by molecular dynamics simulations. Sci. Rep. 10, 1406 (2020). https://doi.org/10.1038/s41598-020-58320-z
• [74]   Bekker,  G.-J.,  Kawabata,  T.,  Kurisu,  G. The biological structure model archive (BSM-Arc): An archive for in silico models and simulations. Biophys. Rev. 12, 371–375 (2020). https://doi.org/10.1007/s12551-020-00632-5