A virtual system-coupled molecular dynamics simulation free from experimental knowledge on binding sites: Application to RNA-ligand binding free-energy landscape

Abstract

Ligand–receptor docking simulation is difficult when the biomolecules have high intrinsic flexibility. If some knowledge on the ligand–receptor complex structure or inter-molecular contact sites are presented in advance, the difficulty of docking problem considerably decreases. This paper proposes a generalized-ensemble method “cartesian-space division mD-VcMD” (or CSD-mD-VcMD), which calculates stable complex structures without assist of experimental knowledge on the complex structure. This method is an extension of our previous method that requires the knowledge on the ligand–receptor complex structure in advance. Both the present and previous methods enhance the conformational sampling, and finally produce a binding free-energy landscape starting from a completely dissociated conformation, and provide a free-energy landscape. We applied the present method to same system studied by the previous method: A ligand (ribocil A or ribocil B) binding to an RNA (the aptamer domain of the FMN riboswitch). The two methods produced similar results, which explained experimental data. For instance, ribocil B bound to the aptamer’s deep binding pocket more strongly than ribocil A did. However, this does not mean that two methods have a similar performance. Note that the present method did not use the experimental knowledge of binding sites although the previous method was supported by the knowledge. The RNA-ligand binding site could be a cryptic site because RNA and ligand are highly flexible in general. The current study showed that CSD-mD-VcMD is actually useful to obtain a binding free-energy landscape of a flexible system, i.e., the RNA-ligand interacting system.

Significance

A novel conformation sampling method, CSD-mD-VcMD, was proposed to compute a free-energy landscape of binding of a flexible receptor and ligand. This method is free from experimental knowledge of the system, and predicts the complex structures. This method was applied to systems consisting of an RNA (the aptamer domain of FMN riboswitch) and a ligand (ribocil A or ribocil B). The resultant free-energy landscape was consistent with that obtained from our previous sampling method, which was done using the knowledge of the experimental complex structure.

 Introduction

In this paper, we develop a novel conformational sampling method based on our previously proposed approach. This method was applied to an RNA system. Below, we first discuss the challenges of studying a structurally flexible system, focusing on an RNA–ligand system. Next, we briefly review computational methods that have contributed to conformational sampling of biomolecular systems with atomistic details. We then describe our previous sampling method and its drawbacks, which naturally lead to the motivation for developing our currently proposed method to overcome the drawbacks.

RNA has been an important target for small-molecular drug candidates [1–8], and efforts for predicting the RNA–ligand complex structure have been accumulated [9–12]. On the other hands, those studies are uncovering that the flexibilities of RNA and ligand make the complex–structure prediction highly difficult [13,14].

The difficulty induced by structural flexibility of molecules can be seen also in proteins: “induced fit” or “induced folding” [15,16], where a part of the protein or ligand varies its structure to form a stabler complex. One of extreme cases for this difficulty of molecular binding is a fuzzy complex [17] or a multi-modal interaction ensemble [18], where many binding modes appear at the room temperature by the thermal fluctuations. We note an interesting phenomenon called a “coupled folding and binding mechanism” that is observed in an intrinsically unstructured protein: An unstructured protein or a partly unstructured protein in an isolated state adopts a well-defined tertiary structure upon binding to a partner molecule to exert its biological activity [19–22]. In previous study, we applied a generalized-ensemble method, the multicanonical molecular dynamics simulation [23] to a system of Ref. 20, which consists of an $\alpha$-helix bundle protein (mSin3B) and a 15-residue segment taken from a neural restrictive silencer factor (NRSF). Although this 15-residue segment is disordered in the isolated state, it adopts an $\alpha$-helix when it binds to the cleft of mSin3B [20]. Therefore, in the initial conformation of our simulation, mSin3B and the segment were dissociated and the segment was disordered. Importantly, in the resultant free-energy landscape, the lowest free-energy basin coincided to the experimentally determined complex structure, where the 15-residue NRSF segment folded into an $\alpha$-helix in the cleft of mSin3B [24]. In addition, the obtained free-energy landscape included various semi-stable free-energy basins with distinct tertiary structures of the system. Another interesting point of this landscape is that some basins were connected by relatively low free-energy pathways, although the other basins were spaced by high free-energy barriers. This means that the free-energy landscape provides possible pathways of conformational changes among different complex structures.

The above examples suggest the usefulness of effective sampling methods to model or predict complex structures of flexible biomolecules and to investigate a mechanism of the complex formation. Methods in a group of generalized ensemble method sample the conformational space of biomolecular system by applying a bias force to the system or modulate a thermodynamic parameter of the to cause readily conformational transitions [25–28]. Importantly, a statistical weight (thermodynamic weight) is assigned to any of sampled conformations (snapshots) based on the physico-chemical principle. This enables us to compute the free-energy landscape (or the conformational distribution of the system) by assigning the density of conformation (potential of mean force) to the conformational space.

One may suspect the accuracy of the resultant free-energy landscape because of the inaccuracy of the force-field parameters, and we admit that no perfect force-field parameters exist so far. However, our previous study, which used a generalized ensemble method named mD-VcMD simulation, developed by us [29], showed qualitative agreement with experimental result [30]: Binding of ribocil B to the aptamer domain of the FMN riboswitch is stronger than that of ribocil A. Our subsequent analysis showed that the free-energy landscape for ribocil B has a funnel-like shape, and the landscape for ribocil A is rather ragged. On the other hands, Ref. 14 noted that the binding poses of ligands to RNA are determined poorly although the overall RNA structures are well reproduced by correct prediction of the coaxial stacking and tertiary contacts in RNA. This difficulty of the complex-structure prediction may increase the value of the generalized ensemble method because our simulation has provided some distinct RNA–ligand complex structures (binding poses) characterized by the statistical weight. Because the binding poses from computation are characterized by thermodynamic weight, the binding-poses can be used for candidates for the complex structure with a priority ranking.

Before explaining to computational methods, here we explain some biophysical significance of the simulation system (an RNA-ligand system) studied in the present paper: The aptamer domain of the FMN riboswitch (RNA) binding to ribocil A or ribocil B (ligand), which was studied using our previous method. In the present study, we applied the new method to the same system to compare the results between the two studies. Details of the system has been presented in our previous paper and the referred paper in the paper [31–36]. Thus, we explain it briefly. The FMN riboswitch exists in mRNA of a bacterium, Bacillus subtilis, which consists of the aptamer domain and an expression platform. The FMN riboswitch relates to gene-expression of riboflavin (also known as vitamin B₂), which is essentially important for the life support of the bacterium. Importantly, binding of a ligand to the deep binding pocket of the aptamer domain causes a structural change in the expression platform although the mechanism for this structural change has not yet been understood well, and finally the riboflavin’s biosynthesis is suppressed by the ligand binding. The stronger the ligand binding, the stronger the biosynthesis suppression. Thus, the investigation of the ligand-receptor binding mechanism is useful for drug discovery of antibiotics.

Analyses to understand statistical mechanics of molecular systems are achieved by taking advantages of the generalized ensemble methods, and those methods were also developed to prevent conformational trapping and then raise the sampling efficiency in the research fields of physics and biophysics. Well-known methods are the replica exchange sampling [37–39], the multicanonical sampling [40–43], and adaptive umbrella sampling [44,45]. Benefit of the generalized ensemble methods are not only their sampling efficiency but also an ability to assign a statistical weight (or thermodynamic weight) to sampled conformations. Furthermore, several generalized ensemble methods have been developed [46–65], which have similarities with the replica exchange, multicanonical, or adaptive umbrella sampling although they were developed independent of replica exchange, multicanonical, or adaptive umbrella sampling. These methods had been used mainly for understanding conformational changes of a single biomolecule in an early stage, and gradually expanded to study on ligand–receptor interactions. In these methods as explained above, a bias force is introduced in the system or a thermodynamic parameter of the system is modulated during the simulation.

Here we note that there is another group of computation methods [66–71] that aim to obtain directly kinetic constants, which are required to investigate conformational transitions among regions/states set in the conformational space of a molecular system: Given two states, say $s_A$ and $s_B$, the free energy difference ($\mathrm{\Delta }G(s_A,s_B)$) between the states is defined by $\mathrm{\Delta }G(s_A,s_B)=-RTln[r_{A\to B}/r_{B\to A}]$, where $r_{A\to B}$ and $r_{B\to A}$ are, respectively, rate constants from stare $s_A$ to $s_B$ and from $s_B$ to $s_A$: $R$ and $T$ are the gas constant and temperature, respectively. By setting the states appropriately, one can obtain a map of $\mathrm{\Delta }G$ in a wide conformational space, i.e., a free-energy landscape.

Because the above-mentioned methods search the conformations in a high-dimensional conformational space, the computation is usually expensive. A way to tackle this problem is to assume a one-dimensional (1D) line, which can be a curve line in the high-dimensional space, and sample the conformations only in the 1D line. We note that the free-energy difference (the difference of the potential of mean force) between the two ends of the 1D line can be calculated with relatively low computational cost. We call this method “line sampling” for convenience in this paper. Importantly, the line sampling has been used to calculate a binding/unbinding free energy (free-energy difference between unbound and bound states) [64,65,72–84]. However, the line sampling cannot discover low-energy regions (basins) out of the 1D pathway in the high-dimensional conformational space. This drawback is crucial when we do not know stable complex structures in advance.

To compute the free-energy landscape of a system consisting of two biomolecules (a ligand and a receptor) in a wide conformational space, which contains unbound configurations, various semi-stable complexes, and native-like complexes, we proposed a generalized ensemble method, a “virtual-system coupled molecular dynamics” (or “mD-VcMD”) simulation method [85,86]. Here we refer to this mD-VcMD as “DSD-mD-VcMD” for convenience (reason of this naming is explained later). The DSD-mD-VcMD method controls the relative positioning of the ligand and the receptor by introducing multiple reaction coordinates (multiple-RCs) and by enhancing the conformational motions in the multiple-RC space. Therefore, this method is suitable to investigate ligand–receptor docking (details are explained later). Furthermore, the DSD-mD-VcMD method can assign a statistical weight (thermodynamic weight) to each of the sampled conformations, and consequently provides a free-energy landscape of molecular binding/unbinding. However, the DSD-mD-VcMD method requires knowledge of the ligand–receptor binding sites to define RCs in advance (details are explained later). This point is a drawback to be applied to a highly flexible system such as the RNA–ligand system because different stable complex structures may be hidden in the conformational space, which is far from the ligand–receptor binding sites postulated in the introduced RCs.

In this study we present a novel sampling method, which uses multiple-RCs free from postulating the ligand–receptor binding sites in advance, and test this method to an RNA–ligand system. Therefore, we do not use the knowledge of the ligand–receptor binding sites, which are known experimentally. Here, we repeat that the present and previous methods are referred to as “CSD-mD-VcMD” and “DSD-mD-VcMD” methods, respectively, and note that the collective term for the two method is “mD-VcMD”. If the present method has proposed a procedure useful to understand the complex formation of the highly flexible molecules (i.e., RNA and ligand) and to present the most and semi-stable complex structures, we believe that the present work has a value for the RNA–ligand complex structure prediction.

 Materials and methods

Because the molecular conformational space to be sampled is vast, we have performed the mD-VcMD simulation iteratively [29]: When the $M$-th iterative stage has finished, the data files for the next iterative stage are prepared, and then the $(M+1)$-th iterative stage starts. Through this cycle, the sampled region of the conformational space is expanded. So far, we have proposed three types of data-preparation methods for mD-VcMD: The “original” [85,86], a “subzone-based” [29,86], and a “genetic-algorithm-guided” [29,86] data-preparation methods. Both the present DSD-mD-VcMD and the previous CSD-mD-VcMD simulations used the subzone-based data-preparation method. See Re. 86 for details of the subzone-based method.

Therefore, the difference of DSD- and CSD-mD-VcMD simulations is not the data-preparation method but the definition (or setting) of the multiple-RCs. It is important that the multiple-RCs used in DSD-mD-VcMD are defined based on a priori knowledge of the ligand–receptor binding sites whereas those for CSD-mD-VcMD are free from binding-site information (reason is explained later).

 Systems studied

We first explain the molecular system studied for the present study. In the previous DSD-mD-VcMD simulation [29], we studied an RNA–ligand system consisting of the aptamer domain of the Fusobacterium nucleatum flavin mononucleotide (FMN) riboswitch (receptor) and ribocil A or ribocil B (ligand). Because ribocil A and B are stereoisomers (Figure 1), we examined two similar systems. We studied these two systems again for the present CSD-mD-VcMD simulation.

Figure 1  Two ligands: (A) ribocil A and (B) ribocil B, which are optical isomers mutually. Ligands are flexible in simulation because all dihedral bonds are treated as rotatable. Red hydrogen atom labeled “H” is the atom different between ribocil A and ribocil B. Four rings in each ribocil are named “R1”, “R2”, “R3”, and “R4” for convenience as shown in figure. This ring naming was used uniformly in this paper and in Ref. 29.

We referred to the system consisting of ribicil A and the aptamer domain of the FMN riboswitch as the “Ribo-A” system, and did to the system of ribocil B and the aptamer domain as the “Ribo-B” system which is consistent with the previous study [29]. Supplementary Figure S1 demonstrates conformations of the Ribo-A and Ribo-B systems, both of which are completely dissociated configurations of the receptor and ligand in an explicit solvent. These were used for the initial configurations for the previous DSD-mD-VcMD study: See Ref. 29 for detailed methods to generate the systems. These configurations were also used for generating the initial conformations for the current CSD-mD-VcMD simulations.

Biophysical properties of the aptamer domain regarding the ligand binding have been studied experimentally [87–90], and an X-ray crystallography reported a complex structure of the aptamer domain with ribocil B (PDB ID: 5c45; resolution 2.93 $\mathring{\mathrm{A}}$) [91]. In this complex, ribocil B binds to the deep pocket of the aptamer domain. Interestingly, an experiment [30] has shown that binding of ribocil B to the aptamer domain is considerably stronger than that of ribocil A, in spite that the two complex structures (the Ribo-A and Ribo-B complexes) are almost identical. Our previous DSD-mD-VcMD simulation provided the following results [29]. i) In the lowest free-energy conformation, both ribocil A and B were in the binding pocket of the aptamer domain, although the ligands shifted by a few angstroms from the position in the X-ray complex structure. ii) Ribocil B bound to the aptamer domain more strongly than ribocil A did, which is consistent with the experimental observation [30] as mentioned above. iii) Both ribocil A and B approached the binding site with passing two tunnels that were completely different to each other. iv) For both ligands, a ligand’s orientational selection mechanism exerted at the entrance of the either tunnel, which is an important mechanism to lead the ligand to the ligand-binding site smoothly. v) There were some ligand-binding modes in the deep binding pocket of the aptamer domain other than the experimentally reported complex structure whereas the overall tertiary structure of the aptamer domain was maintained around the X-ray structure (PDB ID: 5c45).

In the present study, we perform the CSD-mD-VcMD simulation assuming that we do not know either the RNA–ligand X-ray complex structure nor their ligand–receptor binding sites, and compare the results with those from the DSD-mD-VcMD simulation. Remember that DSD-mD-VcMD used the knowledge of the binding sites of the complex [29]. As it will be explained later, both the DSD-mD-VcMD and the CSD-mD-VcMD simulations started from the unbound conformations.

 Multiple-reaction coordinates

Now, we explain multiple reaction coordinates (multiple-RCs), which are set before executing the mD-VcMD simulation [86]. We first explain how to control the RCs for sampling the conformation space widely. After that, we explain a drawback of DSD-mD-VcMD when the ligand-binding site is not known well. Each RC is introduced as a function of atomic coordinates in the system, otherwise we cannot control the ligand–receptor mutual positioning. We refer to a space constructed by the RCs as “RC space”. In the previous studies, we applied distances between two atom groups for the multiple-RCs [29,92–96]. If one of RCs is defined by an inter-molecular distance, then the RC can control association or dissociation of the two molecules. This is why mD-VcMD is suitable for studying the ligand–receptor binding/unbinding. A graphical image of an RC is given in Supplementary Figure S1 of Ref. 86 or Supplementary Figure S4 of Ref. 29. To clarify that the used RC space is a distance space, we have referred to those mD-VcMD simulations as “distance-space-division mD-VcMD” simulations or “DSD-mD-VcMD”.

Whereas we have used multiple RCs usually, Figure 2 presents only a single RC to explain DSD-mD-VcMD. Figure 2 presents schematically a system in that the single RC is represented by the red line between two positions fixed in either the receptor or the ligand. The two molecules go away mutually when the RC increases, and the two molecules approach when the RC decreases.

Figure 2  Two types of distance-based reaction coordinate (RC), which are set to the same ligand–receptor system to perform DSD-mD-VcMD. Although multiple RCs are usually used in DSD-mD-VcMD, only one RC is introduced in each panel for simplicity. (A) RC (thick red line) is defined by distance between ligand and a position near the ligand-binding site (filled red circle) of receptor. (B) RC is distance between the ligand and a position distant from the ligand-binding site of the receptor. We assume that red and black filled circles are the lowest and second-lowest energy sites for ligand binding, respectively, on the receptor surface. In both panels, RC axes are divided into small intervals (called “zones” in this paper), which are numbered from the 1’st to 9’th. Checked shells are mentioned in text.

In mD-VcMD, each RC is divided into small intervals (called “zones”), to which ordinal numbers, such as 1’st, …, 9’th, are assigned as in Figure 2. Actually, we divide the RC axis into more than 9 zones usually. Now, we consider a process that a ligand approaches the ligand-binding site of the receptor from a fully dissociated position. This process is recognized as a series of inter-zone transitions in the scheme of DSD-mD-VcMD. For instance, if the ligand is in the 8’th zone in Figure 2A, the ligand’s approach to its binding site may be a series of inter-zone transitions such as 8’th → 7’th → … → 3’rd → 2’nd. We did not consider the transition from the 2’nd to the 1’st zone in Figure 2A because the 1’st zone is in the interior of the receptor in this figure: An atom collision occurs between the ligand and the receptor in the 1’st zone, which is rejected energetically during a simulation. The inter-zone transitions are controlled by assigning a transition probability between adjacent zones (inter-zone transition probability). Therefore, division of RC into zones is essentially important as discussed in detail in Ref. 86. This transition probability is also explained briefly in this paper later. After finishing the simulation, a statistical weight (thermodynamic weight) is assigned to any sampled conformation by using a re-weighting technique in mD-VcMD [86].

The RC in Figure 2A has a benefit for realizing the molecular binding. The ligand-binding site of the receptor (the red filled circle in the figure) is involved in the 2’nd zone and the volume of the 2’nd zone is smaller than the other zones except for the 1’st zone, which means that the conformation can encounter the ligand-binding site readily when the ligand moves from the 8’th zone to the 2’nd zone. The strategy to search the binding site from a point distant from the binding site is simple: If the zone ordinal No. is set to decrease with reaching the binding site, the ligand reaches the ligand-binding site by moving the ligand so that its zone ordinal No. decreases. An important thing is: The knowledge on the ligand-binding sites is mandatory to set such a convenient RC.

When the ligand-binding site is unknown, one may set an end of the RC randomly at a position of the receptor as exemplified by the “$\times$” mark in Figure 2B, which is unfortunately distant from the ligand-binding site (the red filled circle). In Figure 2B the ligand-binding site is in the 6’th zone, and the volume of the 6’th zone is larger than that of the 2’nd zone. Thus, the ligand in the 6’th zone in Figure 2B has a smaller probability of detecting the ligand-binding site than in the 2’nd zone in Figure 2A. In other words, a longer simulation may be required to search the ligand-binding site when the ligand-binding site is unknown.

As described in the Introduction section, the DSD-mD-VcMD may work inefficiently when applied to a highly flexible system, when the real complex structure may be far from the postulated complex structure, or when the complex structure may have a variety by the conformational fluctuations. Therefore, a robust method free from the postulation on the multiple-RCs is required, which is the “cartesian-space division mD-VcMD” or “CSD-mD-VcMD” method.

 Setting zones directly to the cartesian space

Figure 3 shows the difference of the multiple-RC setting between DSD-mD-VcMD and CSD-mD-VcMD. Figure 3A is a scheme for DSD-mD-VcMD, in which three distance-based RCs (RC1, RC2, and RC3) are introduced, all of which are bridging the receptor and the ligand. Therefore, this 3D-RC space is a distance space. Those multiple-RCs are divided into small zones and the molecular motions are controlled by the inter-zone transition probability as explained above. Because Figure 3A involves three RCs (remember that Figure 2 involved a single RC), we can control elaborately the mutual positioning of the two molecules.

Figure 3  (A) 3D-RCs, RC1, RC2, and RC3, introduced for DSD-mD-VcMD. Although all RCs are defined by distances between receptor and ligand in this figure, some of RCs may be set by an intra-receptor or intra-ligand distance, which specifies the distance between two portions within the receptor. (B) 3 RCs directly assigned to periodic boundary box (i.e., 3D cartesian coordinate space). Terms and symbols appearing in panel are mentioned in text. See figure caption for Figure 2 for red and black filled circles.

An RC can be an arbitrary function expressed by the system’s atomic coordinates [86]. Figure 3B shows a scheme for CSD-mD-VcMD. First, we put the centroid of the aptamer domain of the FMN riboswitch at the body center of the periodic boundary box, as explained in see the Systems–Studied subsection of the Materials–and–Methods section and Supplementally Figure S1. The orientation of the aptamer domain regarding to the periodic box was randomly set. The x-, y-, and z-axes of the cartesian coordinate space are parallel to sides of the periodic boundary box. Therefore, the three RCs $[RC1,RC2,RC3]$ are defined as the 3D cartesian coordinates set to the ligand (marked as “$\times$”). In the present study, the ligand’s position was represented by the centroid of the ring R3 of the ligand (see Figure 1 for the position of ring R3). Therefore, $[RC1,RC2,RC3]=[x_{R_3},y_{R_3},z_{R_3}]$, where $x_{R_3}$, $y_{R_3}$, and $z_{R_3}$ are, respectively, the x, y, and z coordinates of the centroid of ring R3 in the 3D cartesian coordinate space. We did not used the entire-ribocil centroid for the RCs because the ligand rotation around the entire-ribocil centroid does not vary the values of RCs. In other words, the entire-ribocil centroid does not have an enough resolution to discriminate the ligand’s orientational variety. Later in the Result–section, we exemplify a case in that the entire-ligand centroid specifies insufficiently the orientational variety.

Because the shape of the aptamer domain was relatively spherical and compact, the orientation of the aptamer domain was set randomly in the periodic box. On the other hands, one may set the orientation artificially when the receptor’s shape is deviated considerably from the spherical shape (for instance, when the receptor has a rod shape).

 3D-RC space for CSD-mD-VcMD

Here, we explain details of the RC setting for CSD-mD-VcMD. We first assume that the x, y, and z coordinate axes of the 3D cartesian space are parallel to the three sides of the cubic boundary box (Figure 3B). Next, we set the $RC1$, $RC2$, and $RC3$ axes to be parallel to the x, y, and z coordinate axes, respectively: $RC1//x$, $RC2//y$, and $RC3//z$. Then, we divide the 3D cartesian space (i.e., periodic boundary box) into small 3D zones: The yellow region labeled “a zone” in Figure 3B is an example of the zone for CSD-mD-VcMD. We set a periphery region (denoted as “margin” in Figure 3B) in the periodic boundary box, and did not allow [$RC1$, $RC2$, $RC3$] (i.e., the centroid of the ring R3 of the ligand) to enter into the margin during the simulation. In contrast, the ligand moves freely in the region called “the ligand-movable box” in Figure 3B. Therefore, only the ligand-movable box is divided into small 3D zones. A ligand’s density (i.e., the probability of detecting the ring R3 of the ligand) in a zone is calculated only in the ligand-movable box during the mD-VcMD simulation. In both CSD- and DSD-mD-VcMD methods, the inter-zone transition procedure is implemented in the PLUMED plug-in [97] combined with GROMACS [98]. Therefore, the computational cost for CSD-mD-VcMD and DSD-mD-VcMD simulations are almost the same as that for a conventional MD simulation, because the assignment of a conformation to a virtual zone $[RC1,RC2,RC3]$ is done simply by PLUMED. Thus, the GROMACS program runs without modification. We explain later how to control the inter-zone transitional motions of the ligand.

 Zones in Multi-RC space and Inter-zone transitions

Details for division of the 3D-RC space into small 3D zones are explained elsewhere [86]. We note that the space-division procedure is independent of the type of RC whatever the 3D-RCs are distances (Figure 3A) or 3D-cartesan coordinates (Figure 3B). In the mD-VcMD simulation, the conformation is confined in a 3D-RC zone (“current zone”) for a while and then transitions from the current zone to one of neighboring zones. Supplementary Section 1 explains details for the division of the 3D-RC space (Supplementary Figures S2A and B) and procedure for the inter-zone transition (Supplementary Figures S2C). The mD-VcMD simulation produces a distribution function (density function) of the system assigned to each zone, which is sampled during the simulation. A zone is denoted by $\mathit{L}$, which is a 3D vector when multiple RCs are introduced for sampling.

The density assigned to a zone $\mathit{L}$ is expressed as $Q_{cano}(\mathit{L})$. When we specify the iterative stage $M$ to the density, $Q_{cano}(\mathit{L})$ is denoted as $Q_{cano}^{[M]}(\mathit{L})$. The method to calculate $Q_{cano}^{[M]}(\mathit{L})$ from the simulation is explained in Supplementary Equations S1 and 2 of Supplementary Section 2. The inter-zone transition probability is explained in Supplementary Section 3, where Supplementary Equations S3–7 provide rigorous derivation of the inter-zone transition probability, and Supplementary Figure S3 shows schematically the transitions. With traveling among various 3D-RC zones via the inter-zone transitions in the ligand-movable box (Figure 3B), we can obtain the density of the system in each 3D zone, which is explained below.

 CSD-mD-VcMD simulation

As mentioned in the beginning of the Materials–and–Methods section, the mD-VcMD simulation is done iteratively. This is because we set the ligand-movable box (Figure 3B) largely and accordingly the majority of the biophysically important conformations are not sampled in a single simulation run. Because iterative simulations have been performed in DSD-mD-VcMD, we also did the CSD-mD-VcMD in the present study. We numbered the iterations as “the first iterative stage”, “the secondt iterative stage”, and so on. When the $M$-th iterative stage has finished, we calculate a density function $Q_{cano}^{[M]}(\mathit{L})$, which is a probability of detecting the system in zone $\mathit{L}$ computed from the trajectories from the 1-st to $M$-th iterative stages. See Supplementary Section 2 or Ref. 86 for details for calculation of $Q_{cano}^{[M]}(\mathit{L})$.

In DSD-mD-VcMD, to increase the sampling efficiency further, we have performed $N_{run}$ runs in parallel in each iterative stage, where the runs start from different conformations. This parallel-job procedure is also used in CSD-mD-VcMD, and set $N_{run}=1,000$ in this study. In the previous DSD-mD-VcMD simulation we set $N_{run}=1,024$. Those MD trajectories starting from different initial conformations converges to a canonical ensemble in theory [99,100].

We set the simulation length (denoted as $\mathrm{\Delta }{t}_{run}$) of each of $N_{run}$ runs to $2.0ns$ ($1\times {10}^6$ steps with time-step of $2$ fs) for both of the Ribo-A and Ribo-B systems. Therefore, the total length for an iterative stage was $2.0\mu s$ (=$2.0ns\times 1,000runs$). As reported later, we performed 85 iterations for each system. Thus, the total simulation length for each system is $170.0\mu s$ ($=2.0ns\times 1,000\times 85iterations$). We show in the Results–and–Discussion section that the density $Q_{cano}^{[M]}(\mathit{L})$ converges before the 85-th iterative stage. We saved a snapshot at every $1\times {10}^5$ steps ($0.2ns$) in each run of each iterative stage. Thus, the number of saved snapshots from the 85 iterations is $8.5\times {10}^5$ ($=170\mu s/0.2ns$) for each system. Because we set $\mathrm{\Delta }{\tau }_{trns}$ to 200 fs ($100$ steps) in the present study, inter-zone transition events occur ${10}^4$ times ($=1\times {10}^6/100$) in a single run including self-transitions.

The mD-VcMD algorithm was first implemented to an MD simulation program myPresto/omegagene [27] and next to GROMACS [98] with the PLUMED plug-in [97]. Both the CSD- and DSD-mD-VcMD simulations can be done using GROMACS with PLUMED. The simulation conditions were the following: LINCS [101] to constrain the covalent-bond lengths relating to hydrogen atoms, the Nosé–Hoover thermostat [102,103] to control the simulation temperature, and the zero-dipole summation method [104–106] to compute the long-range electrostatic interactions. No structural restraint was applied to the system except for weak position restraints on some receptor atoms to suppress the translational and rotational motions of the receptor as explained later. As mentioned above, the time-step of simulation was set to 2 fs ($\mathrm{\Delta }t=2$ fs) and the simulation temperature to $300K$. Runs were performed on the TSUBAME3.0 (Tokyo Institute of Technology) and TSUBAME4.0 (Institute of Science Tokyo) supercomputers, both of which use GPU. The force field parameters are explained in the next subsection.

The TSUBAME3.0 and 4.0 computer systems are constructed by many cpus and GPUs, and each user can use a part of them depending on how the computer system is crowded. Usually, we have used 30–100 GPUs at the same time and executed 1 run by a combination of 1 cpu and 1 GPU. In the present study a run ($1\times {10}^6$ steps) required about 16 minutes using TSUBAME4.0 system. Then, 85 iterations, each of which consisted of 1,000 runs, took about 9 days for each of Ribo-A and Ribo-B systems executing 30 runs in the same time. Although this computer resource is large, we note that the computation speed and the number of cpu/GPU are increasing rapidly. We are planning to examine larger biomolecules in future.

One might infer that the resultant conformational ensemble is affected by the confinement of the system’s conformation in 3D zones for intervals of $\mathrm{\Delta }{\tau }_{trns}$. This criticism is valid for a conventional MD simulation that uses an artificial potential. In mD-VcMD, however, the bias effect by the zone confinement is removed from the resultant ensemble by a reweighting technique [29,86]. Thus, the thermodynamic weight assigned to a snapshot is a quantity without the confinement. On the other hands, the position restraints introduced in CSD-mD-VcMD may affect the resultant ensemble. Therefore, we set the restraints considerably weak as explained later and checked that the weakness of the restraints. Ref. 29 explains that mD-VcMD is different from umbrella sampling.

 Force fields, ligand-movable box, and weak position restraints on RNA

The systems studied are displayed in Supplementary Figure S1 and the force fields used were explained in Ref. 29. Here, we explain briefly the force fields: Those for the aptamer domain are from the ff99bsc0χOL3 (or χ_OL3) [107] and that for water molecule is from the 3-point optimal point charge (OPC3) water model [108]. The force field for ${Mg}^{2+}$ ion is from a recent study on divalent ions [109] and those for $K^{+}$ and ${Cl}^{-}$ ions are from another recent study on monovalent ions [110]. We have modeled the force field parameters for ribocil A and B by us: The atomic partial charges of ribocil A or B are calculated quantum-chemically using Gaussian09 [111] at the HF/6-31G* level, followed by RESP fitting [112]. Then, the obtained atomic partial charges are incorporated with the general AMBER force field 2 (GAFF2) [113], which was designed to be compatible with conventional AMBER force-fields.

Figure S1 demonstrates that the aptamer domain of the FMN riboswitch was put at the center of the periodic box (cubic shape) filled by a solvent and that the ligand was put near a corner of the periodic boundary box for both systems. Therefore, the ligand and the receptor were completely dissociated. For multiple-run execution in each iterative stage, we generated $N_{run}$ ($=1,000$) conformations from this the dissociated conformation for the initial conformations of the CSD-mD-VcMD simulation. See the legend of Figure S1 for the number of constituent atoms and the size of the periodic boundary box.

Figure 3B shows schematically the periodic boundary box and the ligand-movable box, which is put inside the periodic boundary box. See subsection of “Setting Zones Directly to the Cartesian Space” above for the method to put the aptamer domain in the ligand-movable box. The periodic boundary box was set so that their sides were parallel to the x, y, and z coordinate axes with dimension from $0.0\mathring{\mathrm{A}}$ to $97.3511\mathring{\mathrm{A}}$ for the Ribo-A system and from $0.0\mathring{\mathrm{A}}$ to $97.3304\mathring{\mathrm{A}}$ for the Ribo-B system. The ligand-movable box was set as $[{RC}_{\min }^{(h)},{RC}_{\max }^{(h)}]=[5.0\mathring{\mathrm{A}},92.0\mathring{\mathrm{A}}]$ ($\forall h$) for both systems in the present study. Note that solvent molecules can move in the margin region, although the centroid of the ring R3 of the ligands cannot enter into the margin. Each range of $[{RC}_{\min }^{(h)},{RC}_{\max }^{(h)}]$ was divided into 144 zones ($n_{vs}^{(h)}=144$) for all of the three RC axes (Supplementary Table S1). The lower and upper boundaries of each zone, $[{\zeta }_i^{(h)}{]}_{low}$ and $[{\zeta }_i^{(h)}{]}_{up}$, are listed in Supplementary Table S2. The zone width was set to $1.2\mathring{\mathrm{A}}$ for all zones: $[{\zeta }_i^{(h)}{]}_{up}-[{\zeta }_i^{(h)}{]}_{low}=1.2\mathring{\mathrm{A}}$ ($\forall h$ and $\forall i$). This value of $1.2\mathring{\mathrm{A}}$ was determined from a preliminary study (data not shown). The RC-axis division was the same for all RC axes (Supplementary Tables S1 and S2).

The CSD-mD-VcMD simulation produces the density of the ligand position in each zone at the simulation temperature ($300$ K) in the ligand-movable box. Note that the zones are fixed to the 3D cartesian space (not to the receptor). Therefore, if the receptor can translate and rotate freely in the cartesian space, the ligand’s density may converge to a flat distribution function in the cartesian space. However, in that case, a considerably long simulation is required to obtain the flat density (i.e., convergence). Thus, to avoid the long simulation, we apply weak position restraints on some atoms of the aptamer domain: Eight atoms (P, OP1, OP2, O5', C5', C4', C3', and O3) of 15 nucleotides presented by the stick model in Supplementary Figure S4 are the position-restrained atoms: 120 ($=8\times 15$) atoms. Those atoms were selected not to suppress the motions of the binding pocket of the aptamer domain because those atoms are relatively distant from the ligand–binding site. Furthermore, the position-restraint was not applied to P1 and P4 helices because an experiment has reported that these two helices move largely in thermal fluctuations [91].

We used the following position-restraint function, which has been implemented in GROMACS [98]:

  

$$E_{pr}=c_{pr}\sum _i^{}{\left|{\mathit{r}}_i–{\mathit{r}}_i^{(ref)}\right|}^2,$$(1)

where ${\mathit{r}}_i$ is the position of the $i$-th position-restrained atom and ${\mathit{r}}_i^{(ref)}$ is the reference position (Supplementary Figure S4), to which the restraint is applied. To realize a weak position restraint, we set the coefficient $c_{pr}$ in Equation 1 to a small value: $c_{pr}=10kJ/mol\times n{m}^2=0.024kcal/mol\times {\mathring{\mathrm{A}}}^2$. We show in the Results–and–Discussion section that the aptamer domain fluctuates largely and that both ribocil A and B could bind to the deep binding pocket of the aptamer domain during the simulation.

 Introduction of sub-box: Hierarchical space division

As described above, we set $144$ zones for each RC axis and prepared 1,000 initial conformations for each system. So, we may start the CSD-mD-VcMD simulation. However, the ${144}^3$ zones existing in the 3D-RC space require a large amount of memory in computation. Therefore, to decrease the memory space, we introduce smaller boxes (denoted as “sub-boxes”) in the ligand-movable box and confine runs to move only in a sub-box. The sub-box is a cubic consisting of ${n_{SB}}^3$ zones, which is currently set as $n_{SB}=61$. Thus, the number of 3D zones in the sub-box is considerably smaller than that in the ligand-movable box: $(n_{SB}/n_{vs}^{(h)}{)}^3=(61/144{)}^3=0.076$. Figure 4A illustrates two sub-boxes, ${sb}_1$ and $s{b}_2$, put in the ligand-movable box, and Figure 4B is a two-dimensional version of Figure 4A. The sub-box $s{b}_1$ is fully involved in the ligand-movable box, and $s{b}_2$ is done partially. Therefore, a run confined in $s{b}_1$ can visit all zones (i.e., ${61}^3$ zones) in the sub-box, whereas a run confined in $s{b}_2$ does not sample zones outside the ligand-movable box. The red-colored regions of the zones in Figure 4B are those sampled by the runs. Figure 4C illustrates a situation where many sub-boxes are put in the ligand-movable box.

Figure 4  (A) Two sub-boxes, denoted as $s{b}_1$ and ${sb}_2$, put in ligand-movable box (blue region). Boxes shown by gray lines represent the periodic boundary box. Although the number of sub-boxes is 20 in the actual CSD-mD-VcMD, it is smaller than 20 in this figure. Small cube in each sub-box is the central cube of the sub-box. (B) Two-dimensional representation of panel A. Red area for each sub-box is mentioned in main text. (C) More sub-boxes than 2 are put in the sampling region, where the sub-boxed overlap mutually, and some regions are not occupied by sub-boxes in this panel. Because the sub-boxes are put randomly at the beginning of each iterative stage, the unoccupied regions are occupied with increasing the iterative stage.

Position of each sub-box is determined as follows: First, we pick up a single 3D zone randomly from the ligand-movable box. Next, a sub-box is put so that the picked zone is at the center of the sub-box. The small cubes in Figure 4A and the small squares in Figure 4B are the picked zones, which are located at the centers of sub-boxes. We put 20 sub-boxes at the beginning of each iterative stage of CSD-mD-VcMD in the present study. The sub-boxes may overlap mutually as illustrated in Figure 4C. In an iterative stage, we confined 50 runs in each sub-box. The total number of runs $N_{run}$ done in an iterative stage was 1,000 ($N_{run}=20\times 50$) for the present study. The number of sub-boxes and the number of runs in each sub-box are adjustable parameters.

Supplementary Section 4 explains how to prepare the 50 initial conformations in a sub-box for the $(M+1)$-th iterative stage after the $M$-th iteration has finished. We note that the initial conformations for the first iterative stage should be prepare carefully, which is also explained in Supplementary Section 4 and Supplementary Figure S5 demonstrates actual positioning of ligands for the first iterative stage.

 Ratios of sampled zones

Before starting the simulation, all the zones are not sampled (i.e., all zones are empty). With proceeding the iterative stages, the number of sampled zones increases. After finishing the $M$-th iterative stage, we classify each 3D zone into three types depending on the situation of the zone: “Perfectly-sampled”, “imperfectly-sampled”, and “empty”. In a perfectly-sampled zone, all of the eight subzones in the 3D zone have been sampled during the first to the $M$-th iterative stages. In an imperfectly sampled zone, on the other hand, one or more subzones have not yet sampled during the iterative stages. In an empty zone, no subzone has been sampled yet.

Here, we define two quantities: $R_p^{[M]}$ and $R_{p+ip}^{[M]}$. $R_p^{[M]}$ is the ratio of the perfectly-sampled zones to all the 3D zones (i.e., $[n_{vs}^{(h)}{]}^3={144}^3$) up to the $M$-th iterative stage, and $R_{p+ip}^{[M]}$ is the ratio of perfectly-sampled and imperfectly-sampled zones (or ${144}^3-$ the number of empty zones). To analyze the convergence of $Q_{cano}^{[M]}(\mathit{L})$ with proceeding the iterative stage, we use $R_p^{[M]}$ and $R_{p+ip}^{[M]}$.

 Correlation of spatial pattern between ${\mathit{Q}}_{\mathit{c}\mathit{a}\mathit{n}\mathit{o}}^{\left(\mathit{I}\right)}\mathbf{(}\mathit{L}\mathbf{)}$ and ${\mathit{Q}}_{\mathit{c}\mathit{a}\mathit{n}\mathit{o}}^{\left(\mathit{J}\right)}\mathbf{(}\mathit{L}\mathbf{)}$

Above, we introduced ratios $R_p^{[M]}$ and $R_{p+ip}^{[M]}$ to analyze the sampled zones. Now, we introduce another quantity, a correlation coefficient, $C_{corr}(I;J)$, to analyze a similarity of the 3D-RC spatial patterns between $Q_{cano}^{[I]}(\mathit{L})$ and $Q_{cano}^{[J]}(\mathit{L})$. The value of $C_{corr}(I;J)$ ranges from $-1$ to $+1$. The larger the $C_{corr}(I;J)$, the more similar the $Q_{cano}^{[I]}(\mathit{L})$ to $Q_{cano}^{[J]}(\mathit{L})$. Remember that $Q_{cano}^{[M]}(\mathit{L})$ is the system’s spatial density calculated using from the first to the $M$-th iterative stages. Therefore, $Q_{cano}^{[J]}(\mathit{L})$ expresses a future form of $Q_{cano}^{[I]}(\mathit{L})$. Because $Q_{cano}^{[M]}(\mathit{L})$ converges to a canonical distribution (thermally equilibrated distribution assigned to $\mathit{L}$ at the simulation temperature) with increasing $M$ [99,100], the correlation coefficient can be a good quantity to assess the convergence of $Q_{cano}^{[M]}(\mathit{L})$. The actual form of $C_{corr}(I;J)$ is given in Supplementary Equation S8, and derivation of $C_{corr}(I;J)$ is explained by Supplementary Equations S9–14 of Supplementary Section 5.

 Cumulative spatial density of ligand around receptor

The function $Q_{cano}\left(\mathit{L}\right)$ estimated by the CSD-mD-VcMD simulation expresses the spatial density of the ligand (strictly speaking, the centroid of the R3 ring of ribocil) because $\mathit{L}$ is the position of a 3D zone fixed to the 3D cartesian space. Therefore, the larger the $Q_{cano}\left(\mathit{L}\right)$ in a zone $\mathit{L}$, the larger the ligand’s density in the zone. This gives us an image of the property of the free energy landscape (funnel-like or ragged-like). In this subsection, we analyze quantitatively the degree of the funnel-like property of the free-energy landscape.

For this purpose, we introduce a “cumulative spatial density” of the ligand around the aptamer domain. First, all zones existing in the ligand-movable box are expressed as $\{{\mathit{L}}_{\mathbf{1}},{\mathit{L}}_{\mathbf{2}},{\mathit{L}}_{\mathbf{3}},\cdots ,{\mathit{L}}_{{\mathit{N}}_{\mathit{e}}}\}$ where $N_e$ is the number of 3D zones in the ligand-movable box: $N_e=n_{vs}^{(1)}{n}_{vs}^{(2)}{n}_{vs}^{(3)}={144}^3$. Then, we rearrange the zones so that their densities $\{Q_{cano}\left({\mathit{L}}_{\mathbf{1}}\right),\cdots ,Q_{cano}\left({\mathit{L}}_{\mathit{N}\mathit{e}}\right)\}$ are arranged in a descending order, and the 3D zones are renamed as: ${\mathbf{\Lambda }}_{\mathbf{1}}\mathbf{,}{\mathbf{\Lambda }}_{\mathbf{2}}\mathbf{,}{\mathbf{\Lambda }}_{\mathbf{3}}\mathbf{,}\dots \mathbf{,}{\mathbf{\Lambda }}_{{\mathit{N}}_{\mathit{e}}}$ according to the descending order of densities. By this definition, the following inequality is satisfied:

  

$$Q_{cano}\left({\mathbf{\Lambda }}_{\mathbf{1}}\right)\ge Q_{cano}\left({\mathbf{\Lambda }}_{\mathbf{2}}\right)\ge \dots \ge Q_{cano}\left({\mathbf{\Lambda }}_{{\mathit{N}}_{\mathit{e}}}\right).$$(2)

We call these 3D zones, ${\mathbf{\Lambda }}_{\mathbf{1}}\mathbf{,}{\mathbf{\Lambda }}_{\mathbf{2}}{\mathbf{,}\mathbf{\Lambda }}_{\mathbf{3}}\mathbf{,}\dots \mathbf{,}{\mathbf{\Lambda }}_{{\mathit{N}}_{\mathit{e}}}$, as “density-rearranged zones”.

Next, we accumulate $Q_{cano}\left({\mathbf{\Lambda }}_{\mathit{i}}\right)$ in a descending order of ${\mathbf{\Lambda }}_{\mathit{i}}$ as:

  

$${\mathrm{\Delta }}_n=\sum _{i=1}^n{Q}_{cano}\left({\mathbf{\Lambda }}_{\mathit{i}}\right).$$(3)

Because ${\mathrm{\Delta }}_n$ is probability, it satisfies an inequality: ${0\le \mathrm{\Delta }}_n\le 1$ (remember that $Q_{cano}(\mathit{L})$ is originally normalized by Equation 1). If we set ${\mathrm{\Delta }}_n=0.5$ for instance, then we can assemble zones that contribute to the system’s probability by 50% from the highest-density zones.

Finaly, by setting ${\mathrm{\Delta }}_n$ to a value, we can define a cumulative density $Q_{cum}$ as:

  

$$Q_{cum}({\mathbf{\Lambda }}_{\mathit{i}};{\mathrm{\Delta }}_n)=\{\begin{array}{c}\hfill 1\left(if{\mathbf{\Lambda }}_{\mathit{i}}isinvolvedin{\mathrm{\Delta }}_n\right)\hfill \\ \hfill 0(else)\hfill \end{array}.$$(4)

Because we know the correspondence between ${\mathit{L}}_{\mathit{i}}$ and ${\mathbf{\Lambda }}_{\mathit{k}}$, we can re-write this equation as:

  

$$Q_{cum}({\mathit{L}}_{\mathit{i}};{\mathrm{\Delta }}_n)=\{\begin{array}{c}\hfill 1\left(if{\mathit{L}}_{\mathit{i}}isinvolvedin{\mathrm{\Delta }}_n\right)\hfill \\ \hfill 0(else)\hfill \end{array},$$(5)

The spatial pattern of $Q_{cum}({\mathit{L}}_{\mathit{i}};{\mathrm{\Delta }}_n)$ teaches us the shape of the free-energy landscape: Funnel-like or ragged shape.

 Ligand binding pathways to the binding site of the aptamer domain

In the crystallographic complex structure (PDB ID: 5c45; resolution 2.93 $\mathring{\mathrm{A}}$), ribocil B is buried deeply in the pocket of the aptamer domain [91]. Supplementary Figure S6 shows that the binding ligand is visible from completely opposite directions (front and rear views) of the surface of the aptamer domain. The naming of “front” and “rear” in Figure S6 is not universal but assigned by us arbitrarily in our previous work [29] for convenience. This figure clarifies that the ligand-binding site is in a deep pocket and that the pocket forms a tunnel connecting the front and rear surfaces of the aptamer domain. We cannot judge which is the plausible ligand-binding pathway only from the crystallographic structure.

Our computation [29] showed that the ligand can reach the binding site from both surfaces for both the Ribo-A and Ribo-B systems because the ligand’s density $Q_{cano}(\mathbf{L})$ did not vanish in the entire tunnel region connecting the front and rear surfaces (see Figure 4b and d of Ref. 29). Thus, we named this tunnel as “density tunnel”. On the other hands, we presumed that entering from the front surface leads the ligand to the binding site more smoothly than entering from the rear surface does. Our analyses showed that ligand’s orientation should rotate largely in the density tunnel before reaching the binding sites when the ligand enters from the rear surface [29]. In other words, entering from the rear surface requires climbing a free-energy barrier resulted from the molecular rotation. This point is checked also in the present study.

 Molecular orientation of ribocil around the aptamer domain

Our previous study [29] quantified the ligand’s molecular orientation by two vectors ${\mathit{e}}_{\uparrow }$ and ${\mathit{e}}_{\leftarrow }$, which are fixed in the ligand. Note that the vector ${\mathit{e}}_{\alpha }$ ($\alpha =\uparrow or\leftarrow$) moves according to the ligand’s motion. Thus, the vectors can be good quantities to represent the molecular orientation. We also analyze the ordering of ${\mathit{e}}_{\alpha }$ in the present study. Supplementary Equations S15–17 of Supplementary Section 6 and Supplementary Figure S7 present the method to calculate ${\mathit{e}}_{\alpha }$.

Because the vector ${\mathit{e}}_{\alpha }$ can be assigned to a spatial site in the 3D cartesian space, we can investigate the spatial pattern of ${\mathit{e}}_{\alpha }$ around the aptamer domain. Then, we divided the cartesian space into small cubes (size is ${2\times 2\times 2\mathring{\mathrm{A}}}^3$), whose positions are represented by the body centers of the cubes ${\mathit{r}}_{cube}$. By using the cubes, we can calculate the thermal average of ${\mathit{e}}_{\alpha }$ at $300K$ in each cube: $<{\mathit{e}}_{\alpha }\left({\mathit{r}}_{\mathit{c}\mathit{u}\mathit{b}\mathit{e}}\right)>$ (Equation S14). The larger the $|<{\mathit{e}}_{\alpha }\left({\mathit{r}}_{\mathit{c}\mathit{u}\mathit{b}\mathit{e}}\right)>|$ (i.e., the larger the norm of $<{\mathit{e}}_{\alpha }\left({\mathit{r}}_{\mathit{c}\mathit{u}\mathit{b}\mathit{e}}\right)>$) at site ${\mathit{r}}_{\mathit{c}\mathit{u}\mathit{b}\mathit{e}}$, the more the ordering of ${\mathit{e}}_{\alpha }$ at the site.

We next defined a scalar product, ${SP}_{\alpha }\left({\mathit{r}}_{\mathit{c}\mathit{u}\mathit{b}\mathit{e}}\right)$, between ${\mathit{e}}_{\alpha }^{(X-ray)}$ and $<{\mathit{e}}_{\alpha }\left({\mathit{r}}_{\mathit{c}\mathit{u}\mathit{b}\mathit{e}}\right)>$ (Equation S16), where ${\mathit{e}}_{\alpha }^{(X-ray)}$ is ${\mathit{e}}_{\alpha }$ of the ligand in the X-ray complex structure (PDB ID: 5c45). The larger the ${SP}_{\alpha }\left({\mathit{r}}_{\mathit{c}\mathit{u}\mathit{b}\mathit{e}}\right)$ at site ${\mathit{r}}_{\mathit{c}\mathit{u}\mathit{b}\mathit{e}}$, the similar the $<{\mathit{e}}_{\alpha }\left({\mathit{r}}_{\mathit{c}\mathit{u}\mathit{b}\mathit{e}}\right)>$ to ${\mathit{e}}_{\alpha }^{(X-ray)}$ at this site.

 Ligand’s RMSD and radial distribution-like function

We quantify the ligand’s structural difference between a snapshot and the X-ray complex structure by a well-known quantity “a root-mean-square-deviation” (“$rms{d}_{lig}$”):

  

$$rms{d}_{lig}=\sqrt{\frac{{\mathrm{\Sigma }}_k{|{\mathit{r}}_k-{\mathit{r}}_{\mathit{k}}^{\left(\mathit{X}\right)}|}^2}{N_{lig}}},$$(6)

where ${\mathit{r}}_k$ and ${\mathit{r}}_{\mathit{k}}^{\mathbf{(}\mathit{X}\mathbf{)}}$ are, respectively, the position of the $k$-th heavy atom of the ligand in the snapshot and that of the corresponding atom in the X-ray complex structure, and $N_{lig}$ is the number of heavy atoms in ribocil ($N_{lig}=27$ for both ribocil A and B).

Remember that we applied a weak position restraint on some atoms in the aptamer domain (Equation 1) and that the restrained atoms were set to be distant from the ligand-binding pocket (Supplementary Figure S4). Because of the small force constant $c_{pr}$ in Equation 1, the aptamer domain exhibited whole-molecular motions (we show large fluctuations later), and then the resultant $rms{d}_{lig}$ was affected by the whole molecular motions of the aptamer domain. This means that the whole-molecular motions should be removed from ligand’s fluctuations to analyze $rms{d}_{lig}$, and consequently the removal of the whole-molecular motions is equivalent to superposition of some atoms of the aptamer domain of the snapshots to those in the X-ray complex structure. We selected the superimpose atoms to surround the ligand-binding pocket (Supplementary Figure S8). Note that this superposition was applied to the system after the CSD-mD-VcMD simulation. All the atoms in the snapshots move according to the superposition of the selected atoms, and Equation 6 is calculated. If a snapshot has a large $rms{d}_{lig}$, it means that the ligand structure in the snapshot and that in the X-ray structure are different largely to each other translationally, orientationally, and structurally. Supplementary Section 7 provides some notes for calculating $rms{d}_{lig}$ and the radial distribution-like function ${\rho }_{rd}(r)$.

 Degeneracy problem in DSD-mD-VcMD and alleviation of the problem in CSD-mD-VcMD

Here we note a problem, which we call “degeneracy problem” in this paper, which may induce misleading in analyzing the resultant density $Q_{cano}(\mathit{L})$. In Figures 2 and 3, the lowest and second-lowest binding-energy sites for receptor–ligand binding are shown by the red- and black-filed circles. See the figure caption of Figure 2. Here we note that a serious difficulty, degeneracy problem, may occur, depending on the mutual positioning of low binding-energy sites (the sites of the red- and black-filed circles in this case). The word “degeneracy” is usually used to discuss the energy of the system in physics. However, we use this word to discuss the structure of the system.

In short, when we set the reaction coordinate according to the Figure 2B in the DSD-mD-VcMD sampling, the conformations with the lowest and second-lowest binding-energy sites cannot be discriminated in the density function $Q_{cano}$ (i.e., two conformations are degenerated at the same position in the density function). Then, this may cause misinterpretation of the resultant free-energy landscape. In contrast, in the CSD-mD-VcMD sampling, the lowest and second-lowest binding-energy sites are distinguished (the degeneracy vanished) in the resultant density function. We explained this problem in Supplementary Section 8.

 Results

As mentioned, the present CSD-mD-VcMD and previous DSD-mD-VcMD simulations [29] used the same molecular systems, which consist of the aptamer domain of the FMN riboswitch (receptor) and ribocil A or B (ligand). In DSD-mD-VcMD, we introduced three distance-based RCs, where two of the three RCs were set using the knowledge of the native-complex structure (PDB ID: 5c45; resolution 2.93 $\mathring{\mathrm{A}}$) [91]. Therefore, the ligand approached the position in the native complex position (see Figure 1 of Ref. 29 for actual setting of RCs) when the two RCs became small. The remaining one RC was set to control the cleft width of the aptamer domain, which may also help the ligand to enter into the binding pocket. Consequently, we performed 45-th iterative stages till the density function converged in the 3D distance space. In contrast, in CSD-mD-VcMD, the sub-boxes were set randomly in the ligand-movable box, and performed 85 iterative stages were performed for both systems as explained later.

The extra iterative stages in CSD-mD-VcMD compering to DSD-mD-VcMD may be because we did not use the knowledge of the binding sites of the complex structure. However, we note that if the knowledge is inaccurate or the molecular system is highly flexible, the DSD-mD-vcMD simulation may require more extra iterative stages than the CSD-mD-VcMD simulation. We discuss this point later.

In CSD-mD-VcMD, we set the number of runs, $N_{run}$, performed in each iterative stage as $N_{run}=1,000$, which is slightly smaller than that taken for our previous DSD-mD-VcMD study: $N_{run}=1,024$. The simulation length of each run was exactly the same ($1\times {10}^6$ steps with time step of $2$ fs) with that of the previous study.

 Convergence process

In each iterative stage, 20 sub-boxes were put randomly in the ligand-movable box, and 50 runs were confined in each of the 20 sub-boxes as explained in the Materials–and–Methods section. Sub-panels “Ite $[1:M$]” of Figures 5A and B demonstrate the sampled regions during the first to the $M$-th iterative stages for Ribo-A and the Ribo-B, respectively. Red-regions in the figures represent sampled zones during the iterative stages, and no color region is assigned to zones that were not sampled during the iterative stages (empty-zone region).

Figure 5  Pannels in (A) and (B) are calculated for Ribo-A and Ribo-B systems, respectively. If a zone is colored in red in a sub-panel labeled “Ite $[1:M]$”, then the zone is sampled during the first to the $M$-th iterative stages. Other zones, to which no color is assigned, are non-sampled zones (empty zones) during the iterative stages. The aptamer domain of the FMN riboswitch is put at the center of periodic boundary box, and the P1 and P4 helices are indicated by broken-line circles in panels of “Ite $[1:M]$”.

In the first sub-panel “Ite $[1:1$]” of each of Figures 5A and B, the positions of the individual sub-boxes are seen relatively clearly, although sub-boxes are fused one another with increasing $M$: See sub-panels “Ite $[1:5$]” for instance. With proceeding $M$, sampled regions increased rapidly for both systems. However, the volume occupied by the aptamer domain was not involved in the red contours in both systems because the ligand could not overlap the aptamer domain. On the other hands, the positions of the P1 and P4 helices of the aptamer domain were involved in the red contours when $M$ increased. This is because the two helices fluctuated largely during the simulations. Remember that we did not apply the position restraints to the P1 and P4 helices (Supplementary Figure S4) because an experiment [91] reports that the P1 and P4 helices fluctuate largely.

To analyze the process of the system reaching an equilibrium, we plot the ratios $R_p^{[M]}$ and $R_{p+ip}^{[M]}$ as a function of iterative stage $M$. As shown in Figure 6A, the dependency of the ratios $R_p^{[M]}$ and $R_{p+ip}^{[M]}$ on $M$ was similar between the two systems: Both $R_p^{[M]}$ and $R_{p+ip}^{[M]}$ increased rapidly up to $M=30$ and slowly up to $M=40$. After that, the ratios reached a plateau more slowly in range of $40<M\le 85$. Therefore, one may consider that $Q_{cano}^{[M]}(\mathit{L})$ almost reached an equilibrium for $M\ge 40$ by seeing the similar manner of $R_p^{[M]}$ and $R_{p+ip}^{[M]}$ between the two systems. However, the correlation coefficient $C_{corr}(M;85)$ converged differently between the two systems as shown in the next paragraph. We note that both $R_p^{[M]}$ and $R_{p+ip}^{[M]}$ did not reach 100% with slowly increasing to a value smaller than 100%. This behavior is natural because the ligand cannot crash into the interior of the receptor as shown in Figure 5.

Figure 6  (A) Ratios $R_p^{[M]}$ and $R_{p+ip}^{[M]}$ as a function of iterative stage $M$. Procedure to calculate the ratios is given in Materials–and–Methods section. (B) Correlation coefficient $C_{corr}(M;85)$ (Supplementary Equation S7) as a function of $M$.

Figure 6B plots the correlation coefficient $C_{corr}(M;85)$ (Supplementary Equation S8) as a function of $M$. The coefficient $C_{corr}$ for the Ribo-B system increased quickly in range of $1\le M\le 40$ and reached almost plateau in range of $M>40$. This behavior of $C_{corr}$ is similar to those of the ratios $R_p^{[M]}$ and $R_{p+ip}^{[M]}$ in Figure 6A, although a peak at $M=9$ in Figure 6B for Ribo-B is not recognized in Figure 6A. In contrast, the behavior of $C_{corr}$ for Ribo-A was different from that for Ribo-B. The $C_{corr}$ for Ribo-A increased quickly in $33<M<37$, did slowly in $38\le M,\le 79$, and finally reached plateau in $79\le M\le 84$ (see inset of Figure 6B).

The bumpy behavior of $C_{corr}$ in early iterative stages for both systems (see $C_{corr}$ for $5\le M\le 20$) is because the number of sampled zones used to calculate $C_{corr}$ in the early iterative stages was smaller than that used in latter iterative stages. Remember that an empty zone was eliminated in Supplementary Equation S8. Therefore, in the early iterative stages, $C_{corr}$ involved relatively large statistical errors.

 Spatial density of ligand around the aptamer domain

As explained in the above subsection, we performed the simulation up to the 85–th iterative stage although $Q_{cano}^{[M]}(\mathit{L})$ has been converged up to the 79-th iterative stage for both systems. We rewrite $Q_{cano}^{[M]}(\mathit{L})$ at M=85 as $Q_{cano}(\mathit{L})$ with omitting term “$[M]$” and consider that $Q_{cano}(\mathit{L})$ is an equilibrated special density. Figure 7 illustrates $Q_{cano}(\mathit{L})$, which demonstrates the ligand distribution around the aptamer domain. Apparently, the highest-density spots (i.e., the red-contour regions of Figure 7) were seen in the deep ligand-binding pocket in both the Ribo-A and Ribo-B systems. However, the two systems show different spatial patterns of $Q_{cano}(\mathit{L})$: The red-contour region of Ribo-B was found only in the binding-pocket as a single spot (Figure 7B), although the red-contour region of Ribo-A was split into two spots (Figure 7A) in the binding-pocket. Furthermore, in Ribo-A, the red-contours were found in two regions other than in the binding-pocket: See black broken-line circles in Figure 7A. Therefore, it is likely that the density for ribocil B was localized well in the small volume in the binding-pocket whereas the density for ribocil A dispersed. This tendency of density localization is found for the second-highest iso-density regions (blue-colored contours) and third-highest iso-density regions (yellow-colored contours). These results gave us a strong impression that the free-energy landscape for Ribo-B is funnel like and that for Ribo-A is rugged. In the next subsection, we analyze this impression more quantitatively.

Figure 7  Spatial density $Q_{cano}\left(\mathit{L}\right)$ of ribocil A (A) and ribocil B (B) in the ligand-movable box. Note that $Q_{cano}\left(\mathit{L}\right)$ is normalized by Supplementary Equation S1. Shown density is from the iterative stage 85, in which $Q_{cano}\left(\mathit{L}\right)$ converged (Figure 6B). Differently colored iso-density levels are presented by differently colored contours as indicated in inset. ${\rho }_{\mathit{max}}$ is the maximum density assigned to a zone for each system: ${\rho }_{\mathit{max}}={10}^{\mathrm{-}2.45}$ for Ribo-A and ${\rho }_{\mathit{max}}={10}^{\mathrm{-}2.03}$ for Ribo-B. Positions of ligand-binding pocket, P1 helix, and P4 helix are indicated by words “binding site”, “P1 helix” and “P4 helix”, respectively. Ligand-binding pocket is also indicated by brown broken-line circle. Black broken-line circles are mentioned in main text.

Although Figure 7 demonstrated the difference of the free-energy landscape around the ligand-binding pocket between the two systems, in a broader perspective of the free-energy landscape the two systems have a similarity: The ligand density decreased rapidly when the ligand departed from the surface of the receptor. Therefore, our study shows computationally that the complex formation tends to have a superiority against being in unbound conformations, whereas the system’s state appearing in a real solution depends on the concentration of the ligand and the receptor.

The P1 and P4 helices of the riboswitch fluctuated largely during the simulation, which was also observed in our previous DSD-mD-VcMD simulation [29]. As explained in the Introduction–section, ribocil B and ribocil A are a strong and weak suppressor of the riboflavin biosynthesis. Thus, we expected that the thermal fluctuations of the P1 and P4 helices may be different between the Ribo-A anhd Ribo-B systems. However, unfortunately, the fluctuation patterns of the P1 and P4 helices did not show substantial difference between the two systems (data not shown). This result is also found in our previous DSD-mD-VcMD simulation. The expression platform was removed in both simulations and the lengths of the P1 and P4 helix may be too short to produce the motional difference between the two systems.

 Cumulative density

The above subsection gave us an impression that the free-energy landscape of Ribo-B may be more funnel-like than that of Ribo-A. To analyze more the degree of the funnel feature of the landscape, we introduced the cumulative spatial density $Q_{cum}({\mathit{L}}_{\mathit{i}};{\mathrm{\Delta }}_n)$ (Equation 5). Figure 8 demonstrates $Q_{cum}({\mathit{L}}_{\mathit{i}};{\mathrm{\Delta }}_n)$ at ${\mathrm{\Delta }}_n=0.6,0.8$ and $0.9$. All of the zones in the ligand-movable box are colored in red when we set to ${\mathrm{\Delta }}_n=1.0$, and in contrast, all of them are not colored when we set to ${\mathrm{\Delta }}_n=0.0$. Figure 8A shows that 60% (${\mathrm{\Delta }}_n=0.6$) of the density of ribocil A existed not only in the binding-pocket but also in two other regions on the surface of the aptamer domain in the Ribo-A system. Note that these two regions are also found in Figure 7A. In contrast, Figure 8B demonstrates that the 60% of the ligand’s density localized only in the binding pocket in the Ribo-B system. With increasing ${\mathrm{\Delta }}_n$ to 0.8 and 0.9, the volume of the occupied zones increased for both systems. However, the tendency of density localization in the deep binding pocket was more remarkable in the Ribo-B system than in the Ribo-A system. Therefore, Figure 8 indicates clearly that the free-energy landscape for Ribo-B is funnel-like and that the landscape for Ribo-A was rather rugged. Supplementary Figure S9 presents schematically the free-energy landscapes for the two systems: The Ribo-A system has a rugged free-energy landscape, whereas the Ribo-B does a funnel-like landscape.

Figure 8  Cumulative spatial density $Q_{cum}({\mathit{L}}_{\mathit{i}};{\Delta }_n)$ (Equation 5) at three values of ${\Delta }_n$: ${\Delta }_n=0.6,0.8$ and $0.9$. (A) $Q_{cum}$ for Ribo-A. (B) $Q_{cum}$ for Ribo-B.

Not only the current simulation (CSD-mD-vcMD) but also the previous simulation (DSD-mD-VcMD) showed that ribocil B binds more strongly to the aptamer domain than ribocil A does. In other words, the lowest free-energy basin in the deep binding pocket in the aptamer domain for the Ribo-B system is deeper than that for the Ribo-A system in the free-energy landscape. Importantly, these computational results agree qualitatively with an experiment [30]. As investigated in detail in the previous simulation, this difference of the free-energy landscape between the two systems caused by the difference of inter-molecular stacking: The inter-molecular stacking was formed more firmly in the Ribo-B system than in Ribo-A system. However, the simulation showed that the force-fields involve some inaccuracy. Therefore, we consider that more examination of the force field should be done carefully to discuss quantitatively the free-energy landscape.

 Radial distribution-like function

Figures 7 and 8 indicated that the lowest free energy was assigned to the ligand binding in the ligand-binding pocket. Now, we show concretely the conformations in the ligand-binding pocket. To pick up those conformations we used the radial distribution-like function ${\rho }_{rd}\left(rms{d}_{lig}\right)$ for each system. Figure 9 indicates that ${\rho }_{rd}$ for the Rino-A system has double peaks (i.e., two clusters) at $rms{d}_{lig}\approx 3.75\mathring{\mathrm{A}}$ (denoted as $A_1$ in the figure) and $rms{d}_{lig}\approx 5.25\mathring{\mathrm{A}}$ ($A_2$) and that ${\rho }_{rd}$ for the Ribo-B does a single peak (i.e., single cluster) at $rms{d}_{lig}\approx 6.25\mathring{\mathrm{A}}$ ($B_1$). This figure also shows existence of a small peak (denoted as “native-like”) in the range of $rms{d}_{lig}\le 2.0\mathring{\mathrm{A}}$ for both systems. These clusters were named by us for convenience. Supplementary Figure S10 indicates apparently that the ribocil conformations from those clusters were in the binding-pocket. Supplementary Figure S10A displays the positions of the entire-ligand centroid of snapshots taken from the $A_1$, $A_2$ and native-like clusters from the Rino-A system, and Supplementary Figure S10B does those from the $B_1$ and native-like clusters.

Figure 9  Radial distribution-like function ${\rho }_{\mathrm{r}\mathrm{d}}\left(rms{d}_{lig}\right)$ as a function of a root mean square deviation $rms{d}_{lig}$ (Equation 6) for Ribo-A (blue line) and for Ribo-B (red line). Names of peak of the functions, “peak $A_1$”, “peak $A_2$”, “peak $B_1$”, and “native-like”, are mentioned in main text.

Remember that the high-density region (red-colored contours) of the Ribo-A system were split into two spots (Figure 7A), although the two spots fused into one in Supplementary Figure S10A. This inconsistency is because Figure 7A was generated by the centroid of the R3 ring (i.e., the 3D-RCs were set to the centroid of the R3 ring in CSD-mD-VcMD), whereas Figure S10A was calculated using the position of the entire-ligand centroid. As mentioned in the Materials–and–Methods–section, the ligand’s rotation around the entire-ligand centroid does not move the entire-ligand centroid position largely. In contrast, the R3 ring quantifies the ruggedness of the binding pocket, which results in the distinct high-density spots in the binding pocket in Figure 7A. We presume that it may be better to select a small portion of ligand for the multiple-RCs in the CSD-mD-VcMD.

The entire-centroids of ribicol B in the cluster $B_1$ in Figure 9 forms a single cluster in Figure S10B, and remember that ribocil B exhibited a single density peak in the ligand-binding pocket. This may mean that the binding-pocket of the aptamer domain behaved as a less bumpy tunnel when ribocil B is in the pocket.

 Complex structures in the ligand-binding pocket

We show the complex structures of clusters identified in Figure 9. However, before it, we show the X-ray complex structure (PDB ID: 5c45) in Supplementary Figure S11, which indicates existence of three inter-molecular stackings (A48–R4, G62–R1, and A85–R3). The stackings formed in the X-ray complex structure [91] are referred to as “native stackings” in this paper. We also found two inter-molecular hydrogen bonds (H-bonds): One is from the atom O07 of ribocil B to the O2’ atom of A48 and the other is from the atom O07 of ribocil B to the N7 atom of A99, as shown in Supplementary Figure S11. We show that H-bonds regarding the atom O07 of ribocil are formed frequently to surrounding nucleotides in the clusters.

Now, we show the complex structures of the five clusters in the main text and in the supplementary file: Native-like cluster of Ribo-A is Figure 10, $A_1$ cluster from Ribo-A is in Supplementary Figure S12, and $A_2$ cluster from Ribo-A is in Supplementary Figure S13. The structural comparison of $A_1$ and $A_2$ clusters is shown in Supplementary Figure S14. The native-like cluster of Ribo-B is in Figure 11, and $B_1$ cluster from Ribo-B is in Supplementary Figure S15.

Figure 10  (A) Ten snapshots of ribocil A, picked randomly from the native-like cluster of Ribo-A (Figure 9), which are presented on the frame of the aptamer domain of Figure S1. (B) Ten ribocil-A structures with removing the aptamer domain from panel (A) in green-colored stick model. Magenta-colored stick model is ribocil B in the X-ray complex structure (PDB ID: 5c45). See Figure 1 for four rings R1, R2, R3, and R4 in ribocil. (C) A conformation randomly picked from ten snapshots in panel (A), where the aptamer domain was replaced by one pairing with the picked ribocil A. Three nucleotides, which form native stackings (A48, B62, and A85) with ribocil are shown. Panel (C) provides a native-like H-bond formed in the shown complex structure. The number of H-bonds depends on the conformation in the native-like cluster.

Figure 11  (A) Ten snapshots of ribocil B picked randomly from the native-like cluster of Ribo-B (Figure 9), which are presented on the frame of the aptamer domain of Figure S1. (B) Ten ribocil-A structures with removing the aptamer domain from panel (A) in green-colored stick model. Magenta-colored stick model is ribocil B in the X-ray complex structure (PDB ID: 5c45). See Figure 1 for four rings R1, R2, R3, and R4 in ribocil. (C) A conformation randomly picked from ten snapshots in panel (A), where the aptamer domain was replaced by one pairing with the picked ribocil A. Three nucleotides, which form native stackings (A48, B62, and A85) with ribocil are shown. Panel (C) provides two native-like H-bonds formed in the shown complex structure.

Figure 10 demonstrates complex structures from the native-like cluster of the Ribo-A system. Figure 10A displays 10 snapshots of ribocil A picked randomly from the native-like cluster peak of Ribo-A. These ribocil-A structures overlap well one another and also to ribocil B in the X-ray structure (Figure 10B). Most of these complex structures form three native stackings, and an H-bond regarding the O07 atom of ribocil A are formed (Figure 10C). This figure was drawn by picking up a single snapshot randomly from the ten conformations in Figure 10A.

Figure 11 is complex structures assigned to the native-like cluster from the Ribo-B system. Figure 11A displays 10 snapshots of ribocil B picked randomly from the native-like cluster of Ribo-B. These ribocil-A structures also overlap well one another and also to the X-ray position. Most of these complex structures form three native stackings between ribocil and the aptamer domain, and the two H-bonds regarding the O07 atom of ribocil A are formed (Figure 10C). This figure was drawn by picking up a single snapshot randomly from Figure 11A. Of course, the number of native-like H bonds depends on the snapshots although the majority of the complex structures involves two native-like H-bonds (data not shown).

The native-like clusters (Figures 10 and 11) were involved the native-stackings and the native-like H-bonds, which are presumably stabilizing factors of the native-like complex structure. In the other clusters, we also found various stackings (native and non-native stackings). We also found that the clusters $A_1$, $A_2$, and $B_1$ are frequently characterized by H-bonds regarding the O07 atom of ribocil, although the binding atoms for the H-bonds were not those in the native-like structures. We repeat here that these stackings and H-bonds stabilized the complexes of the clusters $A_1$, $A_2$, and $B_1$.

The present study showed that CSD-mD-VcMD simulation sampled various complex structures, and the native-like cluster for Ribo-B was not the largest. The native-like complex structures were characterized by the native-like stackings and the native-like H-bonds. This result may suggest that the force field is inaccurate: We argued the same problem in our previous DSD-mD-VcMD study. If the force field is more accurate, the native-like complex structure is more stabilized.

Existence of the non-native clusters ($A_1$, $A_2$, and $B_1$) suggests a possibility of the force-field inaccuracy. On the other hands, Supplementary Section 9 discussed that those non-native clusters are stabilized naturally by inter-molecular stackings and inter-molecular H-bonds. As mentioned in the Introduction–section, the flexibilities of RNA and ligand make the complex–structure prediction highly difficult [3,13,14]. Even keeping the coaxial stacking and tertiary contacts in RNA, some RNA–ligand complex structures are possible [14]. Furthermore, the system in our simulation was in solution, whereas the experimental structure of this system is from an X-ray crystallography [91]. Change of the environmental condition may vary the system’s flexibility (this point is discussed later). Based on these arguments, we presume that the complex structure is fluctuating among different conformations in solution. Therefore, a sampling method that can propose a variety of complex structures is important when the system is flexible.

 Ligand’s density in the density tunnel and two binding pathways

In the above subsection, we analyzed ligand’s structures in the ligand-binding pocket. In this subsection, we analyze the ligand’s density in the ligand-binding pocket. As shown in Supplementary Figure S6, there is a tunnel that connects the front and rear surfaces of the aptamer domain. Our previous study showed that non-zero ligand’s density occupies the tunnel, and we named the tunnel “density tunnel” [29]. In this subsection, we analyze again the ligand’s density in the density tunnel in detail.

Figures 12A and B focus on the density tunnel, where two and single high-density spots exist in the Ribo-A and Ribo-A systems, respectively. This figure demonstrates apparently that the front surface and the rear surface of the aptamer domain was connected by non-zero density region (the bright-yellow contours): $Q_{cano}\left(\mathit{L}\right)={\rho }_{\mathit{max}}\times {10}^{-2.55}$ for Ribo-A (Figure 12A) and $Q_{cano}\left(\mathit{L}\right)={\rho }_{\mathit{max}}\times {10}^{-3.47}$ for Ribo-B (Figure 12B). Thus, the ligand can go through the density tunnel from either surface. This result agrees with that from our previous DSD-mD-VcMD study [29]. On the other hands, the density of the cyan contours ($Q_{cano}\left(\mathit{L}\right)={\rho }_{\mathit{max}}\times {10}^{-1.85}$ for Ribo-A and $Q_{cano}\left(\mathit{L}\right)={\rho }_{\mathit{max}}\times {10}^{-2.77}$ for Ribo-B) was disconnected by the red dotted line in each panel of Figure 12. This means that a free-energy barrier exists long the red dotted line in both systems. The red dotted line is at a position slightly shifted from the highest-density spots (the red-contours) toward the rear surface. We note that the free-energy barrier was also found in the DSD-mD-VcMD study [29] and the barrier position was also similar with the present study.

Figure 12  Spatial density of ligand in the ligand-binding pocket. Density is shown at three counter levels indicated by insets. Value for ${\rho }_{\mathit{max}}$ is given in main text. (A) Ribo-A and (B) Ribo-B. Terms appearing in the two panels are mentioned in main text. Position of P1 and P4 helices are indicated by “P1 helix” and “P4 helix”, respectively.

Of course, the connection of the density recovers when we decrease the density level of the cyan contours. We presume from Figure 12 that free-energy barrier height from the lowest free-energy spot for Ribo-B is higher than that for Ribo-A because the ratio of density between the largest density spot ($Q_{cano}\left(\mathit{L}\right)={\rho }_{\mathit{max}}$) and cyan-colored contours is $\frac{{\rho }_{\mathit{max}}}{{\rho }_{\mathit{max}}\times {10}^{-1.85}}\approx 70.8$ for Ribo-A and $\frac{{\rho }_{\mathit{max}}}{{\rho }_{\mathit{max}}\times {10}^{-3.02}}\approx 1047$.3 for Ribo-B.

In the present CSD-mD-VcMD study, $Q_{cano}\left(\mathit{L}\right)$ was computed using all conformations appearing in the entire simulation (i.e., the number of all simulation steps): $85\times {10}^9$ conformations ($={10}^6$ steps $\times {10}^3$ runs $\times 85$ iterations) for each system. In contrast, in DSD-mD-VcMD, $Q_{cano}\left(\mathit{L}\right)$ was computed from snapshots output during the simulation: $4.608\times {10}^6$ snapshots ($=10$ snapshots $\times 1024$ runs $\times 45$ iterations) for each system [29] because a snapshot was output at every ${10}^5$ steps of simulation (i.e., a single run produced only $10$ snapshots). In this sense, the accuracy of $Q_{cano}\left(\mathit{L}\right)$ from CSD-mD-VcMD is considerably higher than that from DSD-mD-VcMD: $(85\times {10}^9)/(4.608\times {10}^6)\approx 2\times {10}^4$. On the other hands, in DSD-mD-VcMD, a small-RC region was sampled selectively: If the small-RC region has not been sampled well up to an iterative stage (say, the $M$-th iterative state), the small-RC region can be sampled selectively in the $(M+1)$-th iterative stage [29,86]. Consequently, DSD-mD-VcMD produced many conformations in the small-RC range than CSD-mD-VcMD did.

On the other hand, we should note that the only quantity usable to analyze the molecular structures in detail is snapshots, which were output taking an interval of ${10}^5$ steps in both CSD-mD-VcMD and DSD-mD-VcMD. Because we performed presently the CSD-mD-VcMD simulation assuming that the native complex structure is unknown, the number of snapshots in the density tunnel is small. This point is discussed in the next subsection.

Last in this subsection, we note that a poorly sampled RC region can be sampled selectively even in the CSD-mD-VcMD simulation: We can set a sub-box so that the poorly sampled region is involved in the sub-box in the next iterative stage. If the native-complex structure is known in advance, the sub-boxes can be set so that the sub-box are near the native-complex structure, although the present CSD-mD-VcMD simulation did not use the knowledge of the native-complex structure.

 Ordering of ligand’s molecular orientations in the density tunnel

The quantities regarding the ligand’s molecular orientation around the aptamer domain were defined in Supplementary Section 6: ${\mathit{e}}_{\alpha }$, $<{\mathit{e}}_{\alpha }>$, and ${SP}_{\alpha }$ ($\alpha =\uparrow$ or $\leftarrow$) at a fixed site ${\mathit{r}}_{\mathit{c}\mathit{u}\mathit{b}\mathit{e}}$ in the periodic boundary box (or ligand-movable box). Remember that the weak position restraints were applied to 15 nucleotides of the aptamer domain (Supplementary Figure S4). The spatial pattern of vector $<{\mathit{e}}_{\uparrow }>$ (Supplementary Equation S15) is provided in Figures 13A and B of the main text. We did not provide the vector $<{\mathit{e}}_{\leftarrow }>$ because it was not informative, which means that $<{\mathit{e}}_{\leftarrow }>$ was randomized in the density tunnel.

Figure 13  Vector $\mathbf{<}{\mathit{e}}_{\mathbf{\uparrow }}\mathbf{>}$ (Supplementary Equation S14) at each site for Ribo-A (A) and Ribo-B (B). Norms of $\mathbf{<}{\mathit{e}}_{\mathbf{\uparrow }}\mathbf{>}$ are presented in inset. Words appearing in the figures are explained in main test. Scalar product $S{P}_{\mathbf{\uparrow }}$ (Supplementary Equation S16) at each site for Ribo-A (C) and for Ribo-B (D). Only two types of $S{P}_{\mathbf{\uparrow }}$ are displayed by different colors depending on their norm (see inset).

Now, we divided the density tunnel into “front tunnel” and “rear tunnel”. Similarly, each of the front and rear tunnels was divided into “entrance” (a cyan broken-line rectangle) and “middle” (a red broken-line rectangle) as shown in Figures 13A and B. Apparently, the ordering of ${\mathit{e}}_{\uparrow }$ was strong in the entrance region for both the front and rear tunnels because relatively large $|<{\mathit{e}}_{\uparrow }>|$ vectors were assigned to the entrances. This strong ordering of ${\mathit{e}}_{\uparrow }$ in the entrances is similar with that found in our previous DSD-mD-VcMD study [29]. It is likely that a mechanism called “ligand’s orientation selection” may work in the current ligand–RNA systems, which was observed computationally in the binding of a medicine and a GPCR [92]. Figures 13C and D demonstrate that the scalar product ${SP}_{\uparrow }$ is relatively large (${|SP}_{\uparrow }|\ge 0.5$) in both entrances for both systems although the ordering in the front entrance is stronger than that in the rear entrance. This result suggests that when the ligand enters into the tunnel from the rear entrance, the ligand should change its molecular orientation largely before the ligand reaches the binding position with the experimentally observed binding pose. It is likely that the orientational rearrangement requires overcoming the free energy barrier (red broken line in Figures 12A and B). We note that this mechanism was also found in our previous DSD-mD-VcMD study [29].

Figures 13A and B showed that the ligand’s orientational ordering was weak in the middle regions (red broken line rectangles). This weak ordering in the middle regions disagrees with our previous DSD-mD-VcMD study: We could detect many $<{\mathit{e}}_{\uparrow }>$ vectors in the middle regions in the study by the effect of selective sampling explained in the above subsection. We have analyzed the snapshots from the current CSD-mD-VcMD simulation and found a simple reason: The number of snapshots output in the middle regions of the density tunnel was small in the present CSD-mD-VcMD simulation. We note that even if a small number of snapshots are output to a zone, it does not necessarily mean that the density assigned to the zones is small [86]. We also noticed that the front middle region had more $<{\mathit{e}}_{\uparrow }>$ vectors than the rear middle region did in Figures 13A and B. This is simply because more snapshots were output to the front middle region than to the rear middle region.

 Discussion

First, we describe the similarities of results between the present and previous [29] methods in many aspects. In both studies, ribocil B bound to the binding pocket of the aptamer domain more strongly than ribocil A did. The current study has shown that the free-energy landscape for the Ribo-B system is funnel-like whereas the landscape for the Ribo-A system is rather rugged. The spatial density $Q_{cano}\left(\mathit{L}\right)$ showed the existence of a tunnel (named density tunnel) connecting the front and rear entrances to the ligand-binding pocket for both systems: It is likely that the ligand can reach the ligand-binding site of the aptamer domain by entering either the front or rear entrance. However, ligand’s entering from the rear entrance requires overcoming a free-energy barrier to reach the ligand-binding position, which agreed with the result from the previous paper.

Furthermore, for both systems, the ligand’s molecular orientation $<{\mathit{e}}_{\uparrow }>$ was ordered in both the front and rear entrances (Figures 13A and B), and $<{\mathit{e}}_{\uparrow }>$ in the front entrance of the density tunnel tended to be parallel to that in the X-ray complex structure (Figures 13C and D). Therefore, the ligand can reach the ligand-binding site smoothly without a large orientational rearrangement when the ligand enters from the front entrance. We presume that this ligand’s orientational ordering suggests a ligand-binding mechanism (ligand’s orientation selection) that was proposed in our work of a ligand (bosentan) binding to the deep binding site of a GPCR protein [92]. In contrast, the ligand entering from the rear entrance tended to have an anti-parallel molecular orientation to that in the X-ray complex structure. Thus, the ligand’s orientation is required to rotate largely to reach the ligand-binding site in that case. In other words, the approach from the rear entrance to the binding-site accompanies a free-energy barrier. These results also agree with our previous study [29].

To investigate $<{\mathit{e}}_{\uparrow }>$ in the middle regions of the density tunnel (see Figures 13A and B for definition of the middle regions), we tried to pick up snapshots from the middle regions. However, the number of snapshots from the middle regions was considerably small. Remember that the sub-boxes were set randomly in the ligand-movable box in the present study and that $<{\mathit{e}}_{\uparrow }>$ was calculated using only the output snapshots (a snapshot outputs every ${10}^5$ simulation steps). Contrarily, in DSD-mD-VcMD, we can sample selectively small-RC regions, which involve the native-like complex conformations, with proceeding the iterative stage. This difference of data collection between CSD- and DSD-mD-VcMD procedures is resulted from whether the native-complex structure is known or unknown. However, we note that selective sampling of poorly-sampled RC zones can be used even when the native complex structure is unknown in CSD-mD-VcMD.

Next, we argue the ligand’s poses detected in the binding pocket. We obtained five distinct clusters: The native-like, $A_1$ and $A_2$ clusters from the Ribo-A system, and the native-like and $B_1$ clusters from the Ribo-B system. As mentioned in Results section, the native-like cluster was less stable than the other clusters in Figure 9. However, existence of the native-like complex structure suggests that the current sampling method has an ability to sample a variety of conformations involving the native-like complex. Thus, the stability of the native complex structure may be stabilized by modifying the force fields.

One important issue of this paper is the effect of the structural flexibility (conformational fluctuations) of RNA and ligand on the complex structure. Recent studies on the RNA–ligand complex have uncovered the importance of the structural flexibility. In other words, the flexibilities of RNA and ligand make the complex–structure prediction highly difficult [3,13,14]. Interestingly, even though the overall shape of RNA is well conserved, different binding poses of the ligand to RNA are possible [14]. Our computation results may relate to these studies. The $A_1$, $A_2$, and $B_1$ clusters were different from the native-like cluster, and each cluster was stabilized by the inter-molecular stackings and inter-molecular hydrogen bonds (Supplementary Figures S12–15). Our results may be example that some stable clusters other than the experimentally observed complex can exist in the binding-pocket of RNA with satisfying a statistically acceptable criterion (these clusters are basins in the free-energy landscape at the room temperature in solution). We calculated the structural fluctuations of the aptamer domain (see Supplementary Equations S18–21 for the calculation method of the fluctuations) and found that the fluctuations were larger than those in the crystal (Supplementary Section 10, Supplementary Table S3, and Supplementary Figure S16), whereas the patterns of the structural fluctuations were very similar between in solution and in crystal. This result also supports that multiple complex forms are possible in the binding pocket in solution.

As discussed above, the prediction of the RNA–ligand complex structure is difficult when the high structural flexibility of the molecules or the force-field inaccuracy is involved in the computation. We admit that we cannot solve this difficulty only from the current results of the DSD-mD-VcMD method. However, the current method can propose possible complex structures (candidates) based on the statistical mechanics (or thermodynamics) with assigning a probability to each candidate. Second, our method can use various force fields by switching an input file of GROMACS (i.e., topology file), which may useful to evaluate accuracy and inaccuracy of the examined force fields. Third, the free-energy landscape can provide not only the complex structures but also the pathways along that various semi-stable complex structures appear. Then, we may obtain a deeper understanding on the RNA-ligand binding process and their structure flexibility.

Last in this paper, we discuss the extensibility and reliability of CSD-mD-VcMD to a system other than the RNA system. As explained, we developed the CSD-mD-VcMD method to be applied to a flexible system and/or a system in that the receptor–ligand interaction sites are not known well. Therefore, one suitable target for our method is an intrinsically disordered protein (or an intrinsically disordered region of an ordered protein) interacting to a flexible ligand. Besides, both the CSD- and DSD-mD-VcMD methods are achieved with using PLUMED. The difference of CSD and DSD methods is only the input file for PLUMED. Thus, we can examine various reaction coordinates other than CSD and DSD by preparing a different PLUMED input file.

 Conflict of interest

The authors declare no conflicts of interest.

 Author contributions

The fundamental idea of CSD-mD-VcMD was generated by J. H., K. K. and S. S. Then, Y. F. proposed to apply the CSD-mD-VcMD method to investigate the RNA-ligand binding mechanism and helped some data analyses to make clearer the meaning of the simulation results. K. K. and S. S. implemented the cartesian-space-division (CSD) procedure to GROMACS and PLUMED. G.-J. B. and N. K. determined the RNA–ligand simulation system to the aptamer domain of the FNM riboswitch and robocil A (or B), and generated the input files for the computation. I. F. and T. T. discussed some methodologically important points and helped some data analyses. Simulation and paper-writing were mainly done by J. H.

 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/77 We are planning to make the code available in near future after taking the necessary steps to address the regulatory issues due to compliance with the Foreign Exchange and Foreign Trade Act.

 Acknowledgements

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, CR-23-02 and CR24-02. J.H. and K.N. acknowledge support by the HPCI System Research Project (Project IDs: hp170020, hp180050, hp190018, hp200063, hp210005, hp220015, hp230003, and hp230011, hp240005, hp240032, hp240167). 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 24ae0121028h0004: 2021–) from AMED and the Japan Biological Informatics Consortium (JBiC). GA-mD-VcMD was performed on TSUBAME3.0 and TSUBAME4.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]   Jones,  D.,  Metzger,  H. J.,  Schatz,  A.,  Waksman,  S. A. Control of Gram-negative bacteria in experimental animals by streptomycin. Science 100, 103–105 (1944). https://doi.org/10.1126/science.100.2588.103
• [2]   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
• [3]   Menichelli,  E.,  Lam,  B. J.,  Wang,  Y.,  Wang,  V. S.,  Shaffer,  J.,  Tjhung,  K. F., et al. Discovery of small molecules that target a tertiary-structured RNA. Proc. Natl. Acad. Sci. U.S.A. 119, e2213117119 (2022). https://doi.org/10.1073/pnas.2213117119
• [4]   Seth,  P. P.,  Miyaji,  A.,  Jefferson,  E. A.,  Sannes-Lowery,  K. A.,  Stephen,  A.,  Osgood,  S. A., et al. Discovery of a new class of RNA-binding small molecules for the hepatitis C virus: Internal ribosome entry site IIA subdomain. J. Med. Chem. 48, 7099–7102 (2005). https://doi.org/10.1021/jm050815o
• [5]   Parsons,  J.,  Castaldi,  M. P.,  Dutta,  S.,  Dibrov,  S. M.,  Wyles,  D. L.,  Hermann,  T. Conformational inhibition of the HCV IRES RNA. Nat. Chem. Biol. 5, 823–825 (2009). https://doi.org/10.1038/nchembio.217
• [6]   Dibrov,  S. M.,  Ding,  K.,  Brunn,  N. D.,  Parker,  M. A.,  Bergdahl,  B. M.,  Wyles,  D. L., et al. Structure of a hepatitis C virus RNA domain in complex with a translation inhibitor reveals a binding mode reminiscent of riboswitches, Proc. Natl. Acad. Sci. U.S.A. 109, 5223–5228 (2012). https://doi.org/10.1073/pnas.1118699109
• [7]   Koehn,  J. T.,  Felder,  S.,  Weeks,  K. W. Innovations in targeting RNA by fragment-based ligand discovery. Curr. Opin. Struc. Biol. 79, 102550 (2023). https://doi.org/10.1016/j.sbi.2023.102550
• [8]   Sivaramakrishnan,  M.,  McCarthy,  K. D.,  Campagne,  S.,  Huber,  S.,  Meier,  S.,  Augustin,  A., et al. Binding to SMN2 pre-mRNA-protein complex elicits specificity for small molecule splicing modifiers. Nat. Commun. 8, 1476, 2017. https://doi.org/10.1038/s41467-017-01559-4
• [9]   Cruz,  J. A.,  Blanchet,  M.-F.,  Boniecki,  M.,  Bujnicki,  J. M.,  Chen,  S.-J.,  Cao,  S. RNA-Puzzles: A CASP-like evaluation of RNA three-dimensional structure prediction. RNA 18, 610–625 (2012). https://doi.org/10.1261/rna.031054.111
• [10]   Miao,  Z.,  Adamiak,  R. W.,  Antczak,  M.,  Batey,  R. T.,  Becka,  A. J.,  Biesiada,  M., et al. RNA-Puzzles Round II: assessment of RNA structure prediction programs applied to three large RNA structures. RNA 21, 1066–1084 (2015). https://doi.org/10.1261/rna.049502.114
• [11]   Miao,  Z.,  Adamiak,  R. W.,  Antczak,  M.,  Batey,  R. T.,  Becka,  A. J.,  Biesiada,  M., et al. RNA-Puzzles Round III: 3D RNA structure prediction of five riboswitches and one ribozyme. RNA 23, 655–672 (2017). https://doi.org/10.1261/rna.060368.116
• [12]   Miao,  Z.,  Adamiak,  R. W.,  Antczak,  M.,  Boniecki,  M. J.,  Bujnicki,  J.,  Chen,  S.-J., et al. E. RNA-Puzzles Round IV: 3D structure predictions of four. ribozymes and two aptamers. RNA 26, 982–995 (2020). https://doi.org/10.1261/rna.075341.120
• [13]    Fulle,  S.,  Holger Gohlke,  H. Analyzing the flexibility of RNA structures by constraint counting. Biophys. J. 94, 4202–4219 (2008). https://doi.org/10.1529/biophysj.107.113415
• [14]   Zhou,  Y.,  Jiang,  Y.,  Chen,  S.-J. RNA–ligand molecular docking: Advances and challenges. WIREs Comput. Mol. Sci. 12, e1571 (2022). https://doi.org/10.1002/wcms.1571
• [15]   Monod,  J.,  Wyman,  J.,  Changeux,  J. P. On the nature of allosteric transitions: A plausible model. J. Mol. Biol. 12, 88–118 (1965). https://doi.org/10.1016/s0022-2836(65)80285-6
• [16]   Spolar,  R. S.,  Record,  M. T. Jr. Coupling of local folding to site-specific binding of proteins to DNA. Science 263, 777–784 (1994). https://doi.org/10.1126/science.8303294
• [17]   Tompa,  P.,  Fuxreiter,  M. Fuzzy complexes: Polymorphism and structural disorder in protein–protein interactions. Trends Biochem. Sci. 33, 2–8 (2008). https://doi.org/10.1016/j.tibs.2007.10.003
• [18]   Kasahara,  K.,  Shiina,  M.,  Higo,  J.,  Ogata,  K.,  Nakamura,  H. Phosphorylation of an intrinsically disordered region of Ets1 shifts a multi-modal interaction ensemble to an auto-inhibitory state. Nucleic Acids Res. 46, 2243–2251 (2018). https://doi.org/10.1093/nar/gkx1297
• [19]   Wright,  P. E.,  Dyson,  H. J. Intrinsically unstructured proteins: Re-assessing the protein structure-function paradigm. J. Mol. Biol. 293, 321–331 (1999). https://doi.org/10.1006/jmbi.1999.3110
• [20]   Nomura,  M.,  Uda-Tochio,  H.,  Murai,  K.,  Mori,  N.,  Nishimura,  Y. The neural repressor NRSF/REST binds the PAH1 domain of the Sin3 corepressor by using its distinct short hydrophobic helix. J. Mol. Biol. 354, 903–915 (2005). https://doi.org/10.1016/j.jmb.2005.10.008
• [21]   Dyson,  H. J.,  Wright,  P. E. Intrinsically unstructured proteins and their functions. Nat. Rev. Mol. Cell Biol. 6, 197–208 (2005). https://doi.org/10.1038/nrm1589
• [22]   Sugase,  K.,  Dyson,  H.,  Wright,  P. Mechanism of coupled folding and binding of an intrinsically disordered protein. Nature 447, 1021–1025 (2007). https://doi.org/10.1038/nature05858
• [23]   Nakajima,  N.,  Nakamura,  H.,  Kidera,  A. Multicanonical ensemble generated by molecular dynamics simulation for enhanced conformational sampling of peptides. J. Phys. Chem. B 101, 817–824 (1997). https://doi.org/10.1021/jp962142e
• [24]   Higo,  J.,  Nishimura,  Y.,  Nakamura,  H. A free-energy landscape for coupled folding and binding of an intrinsically disordered protein in explicit solvent from detailed all-atom computations. J. Am. Chem. Soc. 133, 10448–10458 (2011). https://doi.org/10.1021/ja110338e
• [25]   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
• [26]   Ikebe,  J.,  Umezawa,  K.,  Higo,  J. Enhanced sampling simulations to construct free energy landscape of protein-partner substrate interaction. Biophys. Rev. 8, 45–62 (2016). https://doi.org/10.1007/s12551-015-0189-z
• [27]   Kasahara,  K.,  Terazawa,  H.,  Takahashi,  Y.,  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
• [28]   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
• [29]   Higo,  J.,  Bekker,  G.-J.,  Kamiya,  N.,  Fukuda,  I.,  Fukunishi,  Y. Binding free-energy landscapes of small molecule binder and non-binder to FMN riboswitch: All-atom molecular dynamics. Biophys. Physicobiol. 20, e200047 (2023). https://doi.org/10.2142/biophysico.bppb-v20.0047
• [30]   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
• [31]   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
• [32]   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
• [33]   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
• [34]   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
• [35]   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
• [36]   Serganov,  A.,  Nudler,  E. A decade of riboswitches. Cell 152, 17–24 (2013). https://doi.org/10.1016/j.cell.2012.12.024
• [37]   Swendsen,  R. H.,  Wang,  J. S. Replica Monte Carlo simulation of spin-glasses. Phys. Rev. Lett. 57, 2607–2609 (1986). https://doi.org/10.1103/PhysRevLett.57.2607
• [38]   Hukushima,  K.,  Nemoto,  K. Exchange Monte Carlo method and application to spin glass simulations. J. Phys. Soc. Jpn. 65, 1604–1608 (1996). https://doi.org/10.1143/JPSJ.65.1604
• [39]   Sugita,  Y.,  Okamoto,  Y. Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 314, 141–151 (1999). https://doi.org/10.1016/S0009-2614(99)01123-9
• [40]   Berg,  B.,  Neuhaus,  T. Multicanonical ensemble: A new approach to simulate first-order phase transitions. Phys. Rev. Lett. 68, 9–12 (1992). https://doi.org/10.1103/PhysRevLett.68.9
• [41]   Hansmann,  U. H. E.,  Okamoto,  Y. Prediction of peptide conformation by multicanonical algorithm: New approach to the multiple-minima problem. J. Comp. Chem. 14, 1333–1338 (1993). https://doi.org/10.1002/jcc.540141110
• [42]   Kidera,  A. Enhanced conformational sampling in Monte Carlo simulations of proteins: Application to a constrained peptide. Proc. Natl. Acad. Sci. U.S.A. 92, 9886–9889 (1995). https://doi.org/10.1073/pnas.92.21.9886
• [43]   Nakajima,  N.,  Nakamura,  H.,  Kidera,  A. Multicanonical ensemble generated by molecular dynamics simulation for enhanced conformational sampling of peptides. J. Phys. Chem. 101, 817–824 (1997). https://doi.org/10.1021/jp962142e
• [44]   Paine,  G. H.,  Scheraga,  H. A. Prediction of the native conformation of a polypeptide by a statistical–mechanical procedure. I. Backbone structure of enkephalin. Biopolymers 24, 1391–1436 (1985). https://doi.org/10.1002/bip.360240802
• [45]   Mezei,  M. Adaptive umbrella sampling: Self-consistent determination of the non-Boltzmann bias. J. Comput. Phys. 68, 237–248 (1987). https://doi.org/10.1016/0021-9991(87)90054-4
• [46]   Torrie,  G. M.,  Valleau,  J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 23, 187–199 (1977). https://doi.org/10.1016/0021-9991(77)90121-8
• [47]   Lee,  J. New Monte Carlo algorithm: Entropic sampling. Phys. Rev. Lett. 71, 211–214 (1993). https://doi.org/10.1103/PhysRevLett.71.211
• [48]   Fukunishi,  Y.,  Tateishi,  T.,  Suzuki,  M. Octane/water interfacial tension calculation by molecular dynamics simulation. J. Colloid Interface Sci. 180, 188–192 (1996). https://doi.org/10.1006/jcis.1996.0288
• [49]   Iba,  Y.,  Chikenji,  G.,  Kikuchi,  M. Simulation of lattice polymers with multi-self-overlap ensemble. J. Phys. Soc. Jpn 67, 3327–3330 (1998). https://doi.org/10.1143/jpsj.67.3327
• [50]   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
• [51]   Darve,  E.,  Pohorille,  A. Calculating free energies using average force. J. Chem. Phys. 115, 9169–9183 (2001). https://doi.org/10.1063/1.1410978
• [52]   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
• [53]   Fukunishi,  H.,  Watanabe,  O.,  Takada,  S. On the Hamiltonian replica exchange method for efficient sampling of biomolecular systems: Application to protein structure prediction. J. Chem. Phys. 116, 9058–9067 (2002). https://doi.org/10.1063/1.1472510
• [54]   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
• [55]   Deng,  Y.,  Roux,  B. Computations of standard binding free energies with molecular dynamics simulations. Phys. Chem. B 113, 2234–2246 (2009). https://doi.org/10.1021/jp807701h
• [56]   Moritsugu,  K.,  Terada,  T.,  Kidera,  A. Scalable free energy calculation of proteins via multiscale essential sampling. J. Chem. Phys. 133, 224105 (2010). https://doi.org/10.1063/1.3510519
• [57]   Itoh,  S. G.,  Okumura,  H. Replica-permutation method with the Suwa–Todo algorithm beyond the replica-exchange method. J. Chem. Theory Comput. 9, 570–581 (2013). https://doi.org/10.1021/ct3007919
• [58]   Ikebe,  J.,  Sakuraba,  S.,  Kono,  H. Adaptive lambda square dynamics simulation: An efficient conformational sampling method for biomolecules. J. Comput. Chem. 35, 39–50 (2014). https://doi.org/10.1002/jcc.23462
• [59]   Ikebe,  J.,  Sakuraba,  S.,  Kono,  H. H3 histone tail conformation within the nucleosome and the impact of K14 acetylation studied using enhanced sampling simulation. PLoS Comput. Biol. 12, e1004788 (2016). https://doi.org/10.1371/journal.pcbi.1004788
• [60]   Peter,  E. K.,  Shea,  J.-E. A hybrid MD-kMC algorithm for folding proteins in explicit solvent. Phys. Chem. Chem. Phys. 16, 6430–6440 (2014). https://doi.org/10.1039/C3CP55251A
• [61]   Dasgupta,  B.,  Nakamura,  H.,  Higo,  J. Flexible binding simulation by a novel and improved version of virtual-system coupled adaptive umbrella sampling. Chem. Phys. Lett. 662, 327–332 (2016). https://doi.org/10.1016/j.cplett.2016.09.059
• [62]   Ekimoto,  T.,  Ikeguchi,  M. Multiscale molecular dynamics simulations of rotary motor proteins. Biophys. Rev. 10, 605–615 (2018). https://doi.org/10.1007/s12551-017-0373-4
• [63]   Knight,  J.L.,  Brooks,  C. L., III. \lambda$\lambda$-dynamics free energy simulation methods. J. Comput. Chem. 30, 1692–1700 (2009). https://doi.org/10.1002/jcc.21295
• [64]   Kirkwood,  J. G. Statistical mechanics of fluid mixtures. J. Chem. Phys. 3, 300–313 (1935). https://doi.org/10.1063/1.1749657
• [65]   Gelman,  A.,  Meng,  X.-L. Simulating normalizing constants: From importance sampling to bridge sampling to path sampling. Stat. Sci. 13, 163–185 (1998). https://doi.org/10.1214/ss/1028905934
• [66]   Elber,  R.,  Fathizadeh,  A.,  Ma,  P.,  Wang,  H. Modeling molecular kinetics with milestoning. Wiley Interdiscip. Rev. Comput. Mol. Sci. 11, e1512 (2021). https://doi.org/10.1002/wcms.1512
• [67]   Ojha,  A. A.,  Votapka,  L. W.,  Amaro,  R. E. Advances and challenges in milestoning simulations for drug−target kinetics. J. Chem. Theory Comput. 20. 9759–9769 (2024). https://doi.org/10.1021/acs.jctc.4c01108
• [68]   Zuckerman,  D. M.,  Chong,  L. T. Weighted ensemble simulation: Review of methodology, applications, and software. Annu. Rev. Biophys. 46, 43−57 (2017). https://doi.org/10.1146/annurev-biophys-070816-033834
• [69]   Noé,  F.,  Feliks Nüske,  F. A variational approach to modeling slow processes in stochastic dynamical systems. Multiscale Model. Simul. 635–655 (2013). https://doi.org/10.1137/110858616
• [70]   Dellago,  C.,  Bolhuis,  P. G.,  Csajka,  F. S.,  Chandler,  D. Transition path sampling and the calculation of rate constants. J. Chem. Phys. 108, 1964–1977 (1998). https://doi.org/10.1063/1.475562
• [71]   Harada,  R.,  Kitao,  A. Parallel cascade selection molecular dynamics (PaCS-MD) to generate conformational transition pathway. J. Chem. Phys. 139, 035103 (2013). https://doi.org/10.1063/1.4813023
• [72]   Zwanzig,  R. W. High‐temperature equation of state by a perturbation method. I. Nonpolar gases. J. Chem. Phys. 22, 1420–1426 (1954). https://doi.org/10.1063/1.1740409
• [73]   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
• [74]   Merz,  K. M. Jr.,  Kollman,  P. A. Free energy perturbation simulations of the inhibition of thermolysin: prediction of the free energy of binding of a new inhibitor. J. Am. Chem. Soc. 111, 5649–5658 (1989). https://doi.org/10.1021/ja00197a022
• [75]   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
• [76]   Bartels,  C. Analyzing biased Monte Carlo and molecular dynamics simulations. Chem. Phys. Lett. 331, 446–454 (2000). https://doi.org/10.1016/S0009-2614(00)01215-X
• [77]   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
• [78]   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
• [79]   Bekker,  G.-J.,  Araki,  M.,  Oshima,  K.,  Okuno,  Y.,  Kamiya,  N. Exhaustive search of the configurational space of heat-shock protein 90 with its inhibitor by multicanonical molecular dynamics based dynamic docking. J. Comput. Chem. 41, 1606–1615 (2020). https://doi.org/10.1002/jcc.26203
• [80]   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
• [81]   Bekker,  G.-J.,  Araki,  M.,  Oshima,  K.,  Okuno,  Y.,  Narutoshi,  N. Accurate binding configuration prediction of a G‑protein-coupled receptor to its antagonist using multicanonical molecular dynamics-based dynamic docking. J. Chem. Inf. Model. 61, 5161−5171 (2021). https://doi.org/10.1021/acs.jcim.1c00712
• [82]   Bekker,  G.-J.,  Fukuda,  I.,  Higo,  J.,  Fukunishi,  Y.,  Kamiya,  N. Cryptic‑site binding mechanism of medium‑sized Bcl‑xL inhibiting compounds elucidated by McMD‑based dynamic docking simulations. Sci. Rep. 11, 5046 (2021). https://doi.org/10.1038/s41598-021-84488-z
• [83]   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
• [84]   Bekker,  G.-J.,  Kamiya,  N. N‑terminal-driven binding mechanism of an antigen peptide to human leukocyte antigen-A*2402 elucidated by multicanonical molecular dynamic-based dynamic docking and path sampling simulations. J. Phys. Chem. B 125, 13376−13384 (2021). https://doi.org/10.1021/acs.jpcb.1c07230
• [85]   Hayami,  T.,  Higo,  J.,  Nakamura,  H.,  Kasahara,  K. Multidimensional virtual-system coupled canonical molecular dynamics to compute free-energy landscapes of peptide multimer assembly. J. Comput. Chem. 40, 2453–2463 (2019). https://doi.org/10.1002/jcc.26020
• [86]   Higo,  J.,  Kusaka,  A.,  Kasahara,  K.,  Kamiya,  N.,  Hayato,  H.,  Xie,  Q., 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
• [87]   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
• [88]   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
• [89]   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
• [90]   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
• [91]   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
• [92]   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
• [93]   Higo,  J.,  Kawabata,  Y.,  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
• [94]   Higo,  J.,  Kasahara,  K.,  Bekker,  G.-J.,  Ma,  B.,  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
• [95]   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 Omega 8, 14572−14582 (2023). http://doi.org/10.1021/acsomega.3c00364
• [96]   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
• [97]   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
• [98]   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
• [99]   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
• [100]   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
• [101]   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
• [102]   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
• [103]   Hoover,  W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 31, 1695–1697 (1985). https://doi.org/10.1103/PhysRevA.31.1695
• [104]   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
• [105]   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
• [106]   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
• [107]   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
• [108]   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
• [109]   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
• [110]   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
• [111]   Frisch,  M. J.,  Trucks,  G. W.,  Schlegel,  H. B.,  Scuseria,  G. E.,  Robb,  M. A.,  Cheeseman,  J. R., et al. (2009). Gaussian 09, Revision D.01. Gaussian, Inc., Wallingford CT.
• [112]   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
• [113]   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