Flexible binding simulation by a novel and improved version of virtual-system coupled adaptive umbrella sampling

Chemical Physics Letters 662 (2016) 327–332

Contents lists available at ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

Research paper

Flexible binding simulation by a novel and improved version of virtual-system coupled adaptive umbrella sampling Bhaskar Dasgupta a,b,⇑, Haruki Nakamura a, Junichi Higo a a b

Institute for Protein Research, Osaka University, 3-2, Yamadaoka, Suita, Osaka 565-0871, Japan Technology Research Association for Next Generation Natural Products Chemistry, 2-3-26 Aomi, Koto-ku, Tokyo 135-0064, Japan

a r t i c l e

i n f o

Article history: Received 4 July 2016 In final form 20 September 2016 Available online 21 September 2016 Keywords: Adaptive umbrella sampling (AUS) Virtual-system coupled AUS Markov approximated probability distribution Force-bias parameterization Binding free-energy landscape

a b s t r a c t Virtual-system coupled adaptive umbrella sampling (VAUS) enhances sampling along a reaction coordinate by using a virtual degree of freedom. However, VAUS and regular adaptive umbrella sampling (AUS) methods are yet computationally expensive. To decrease the computational burden further, improvements of VAUS for all-atom explicit solvent simulation are presented here. The improvements include probability distribution calculation by a Markov approximation; parameterization of biasing forces by iterative polynomial fitting; and force scaling. These when applied to study Ala-pentapeptide dimerization in explicit solvent showed advantage over regular AUS. By using improved VAUS larger biological systems are amenable. Ó 2016 Elsevier B.V. All rights reserved.

1. Introduction Molecular dynamics simulation (MD) is useful to understand the complexities of a polyatomic system, however inefficient because of the complicated force-field equations and explicit treatment of solvents with the system of interest [1]. This makes conformational sampling computationally expensive and sometimes inadequate to simulate long-timescale events [2]. Therefore, understanding microscopic details of those rare events require special generalized ensemble techniques [3–8], for example adaptive umbrella sampling (AUS) [9–12]. Understanding details of rare events are particularly crucial for biomolecular systems, because many of the molecular functions related to the survival of an organism are slow processes (e.g. folding, binding, conformational change) [13,14]. Yet, those microscopic details are useful in medical research [15]. In this context, binding is one of the central functions of biomolecules, which can be regarded as a slow process of sampling between two extremes – bound and unbound states [2]. For such a slow process thermodynamic details may include large statistical error. Supposing a free-energy difference of 3 kcal/mol between bound (major basin) and unbound (minor basin) states, one samples 1.7
107 snapshots of bound state per snapshot of unbound state [16], which means large statistical error may be included in binding ⇑ Corresponding author at: Institute for Protein Research, Osaka University, 3-2, Yamadaoka, Suita, Osaka 565-0871, Japan. E-mail address: [email protected] (B. Dasgupta). http://dx.doi.org/10.1016/j.cplett.2016.09.059 0009-2614/Ó 2016 Elsevier B.V. All rights reserved.

free-energy landscape. Additionally, when the bound state includes multiple basins the minor basins may also be insignificantly sampled. To overcome this sampling problem AUS or similar techniques are used. AUS is one of the flat histogram techniques, quite akin to sampling by multicanonical MD (McMD) [2,9]. In AUS, energy of the simulated system is biased by potential of mean force (PMF) along a chosen reaction coordinate ðkÞ. This bias is applied at a given temperature ðTÞ, at which the simulation is performed; therefore multiple basins are sampled at that temperature. By iterative MD initially unknown PMF is updated till convergence. The converged PMF corresponds to an uniform probability distribution of the reaction coordinate ðP obs ðk; TÞÞ. However, we observed that many cases of iterations fail to converge because of hysteresis of Pobs ðk; TÞ. One reason for the hysteresis is that there is theoretically infinite number of microstates in a given k-slice ðk; k þ dkÞ, therefore in a given iteration only a tiny subset of the microstates may be sampled, and in the successive iteration a different subset is sampled. This can be improved by sampling a k region exhaustively. This motivated us to design virtual-system coupled AUS or VAUS [16]. The method VAUS has a precursor named virtual-system coupled multicanonical MD (VMcMD) [4,17–20]. In VAUS we simulate the biomolecular system of interest coupled to an arbitrary virtual system composed of virtual states (see Section 2). Moreover, VMcMD can be generalized by VAUS, i.e., the energy is adopted as the reaction coordinate in VMcMD [21].

328

B. Dasgupta et al. / Chemical Physics Letters 662 (2016) 327–332

Previously VMcMD and VAUS were successfully used for dimerization of a endothelin-1 derivative [19] and Ab-peptide [16], respectively. However, there were a few theoretical limitations on sampling during iterative process. Moreover, we are also interested in qualitative difference between ensembles obtained by VAUS and AUS. In VAUS we have used a reaction coordinate defined by the distance between geometric centers of a pair of clusters of atoms. The clusters of atoms are included in two different polypeptide chains (Fig. 1), thereby k can be regarded as a inter-molecular distance. This k is simple to realize, but sensitive to the variation in biasing potential. This sensitivity requires that one should be careful in estimating the biasing potential. In AUS and VAUS this estimation is crucial to the convergence of simulation and is one of the primary focuses of the current work. To this end, we have improved the calculation of bias by (1) a novel calculation method of probability distribution by using Markov state modeling, (2) refinement of fitting to better parameterize the probability distribution, and (3) de-sensitizing biasing force. The result was compared with that from a regular AUS procedure.

2. Theory 2.1. VAUS First we briefly introduce VAUS method, whose details are explained in supplementary material S1. In VAUS a virtual degree of freedom is introduced to the real biomolecular system ðrÞ [2,16]. We assume that the virtual degree of freedom is a property of a virtual system ðv Þ and takes discrete states v k , where k ¼ 1; . . . ; m and m is the number of virtual states. In VAUS along a reaction-coordinate k, the state of the entire system is designated by k and v k . Virtual state v k is specified by a window ½kk;min ; kk;max  along the k axis (supplementary table ST1), and the configuration is restrained in the window if the system is in v k . VAUS is a sampling method where the system moves along the k axis continuously and hops among the virtual states (see supplementary material S1 and supplementary Fig. SF1a). The virtual and molecular systems interact to each other [19,22] because potential energy function

involves v k and molecular configurations x. The effective potential energy of the biomolecular system is

EVAUS ðx; v k Þ ¼ Eðx; v k Þ þ RT ln ½Pc ðk; v k ; TÞ=gðv k ; kÞ;

ð1Þ

where E is the original potential energy of the molecular system at configuration x using an empirical force-field, R is the gas constant. The term Pc is a canonical probability distribution as a function of k in virtual state v k , and g is a function introduced to restrain the biomolecular system in the window of virtual state v k . If Pc is given in advance, EVAUS is calculated by Eq. (1), and VAUS simulation is performed using EVAUS . However, Pc is unknown a priori and therefore, it is determined by an iterative procedure, where Pc is updated by using the following equation [2]:

ln Pc ðk; v k ; TÞðcþ1Þ ¼ ln Pc ðk; v k ; TÞðcÞ þ ln PVAUS ðk; v k ; TÞðcÞ ;

ð2Þ

where PVAUS ðk; v k ; TÞ is observed probability distribution of k in virtual state v k . The detail of this update procedure is outlined in supplementary material S1. 2.2. Improvements of VAUS Although VAUS is a powerful sampling method [16], some improvements are required. The improvements are aimed at maintenance of detailed balance [2] and better approximation of the distribution functions. In Monte Carlo sampling the detailed balance is automatically satisfied. In MD, however, it is assumed that the detailed balance is satisfied if a trajectory is significantly long. In enhanced sampling schemes, such as VAUS, additional attentions are required. 2.2.1. Markov approximated probability distribution The trajectory (sequence of recorded k) obtained from an iteration c may not be significantly long, and then, the observed distribution may be different from the equilibrium distribution. To obtain the equilibrated distribution, a method (e.g. Markov state model; MSM) that imposes the detailed balance to the trajectory is required. For this purpose, we assume that even for a short trajectory a configuration propagator follows Markov rules [23–25]. In fact MSM is one of the successful methods to understand equilibrium dynamics [26–29]. In short, P VAUS ðk; v k ; TÞ is converted to

PMA VAUS ðk; v k ; TÞ in which detailed balance is imposed (supplementary Fig. SF2). Next we explain details of the technique. In VAUS, the configuration moves along the k axis by keeping its virtual state at v k for f int steps of simulation (supplementary information S1, [22]), and the system can hop to a different virtual state only at the f int -th step. Hence, a complete VAUS trajectory can be broken to a set of subtrajectories of length f int and each subtrajectory can be assigned to a virtual state. Accordingly, MSM consists of two layers; one layer is composed of transitions among bins along the k axis for each virtual state (within a subtrajectory), and the other layer is composed of transitions among virtual states along virtual degree of freedom. A virtual-state transition probability matrix (X) can be calculated from the virtual-state sequence, whose ðl; kÞ-th element is the virtual-state transition probability Xðv l jv k Þ from v k to v l (Fig. SF2b). In the other layer, by collecting subtrajectories belonging to v k , we calculate an intra-state transition probability Fig. 1. Equilibrated two Ace-(Ala)5-Nme peptides system (without solvent) used for initial conformation of simulation. Traces of peptides are shown in magenta, and atoms are in stick representation with CPK color scheme. Cyan and red spheres are Ca- and carbonyl O-atoms, respectively. Geometric center of cyan and red-spheres for each peptide is shown by black sphere, and inter-peptide distance ðkÞ, shown by broken line, is the distance between the geometric centers (in Å). Thick black arrows define inter-Ca vector from the first (residue 2) to fifth (residue 6) Ala residues. Angle between these two vectors is referred by h.

matrix Xkintra whose matrix element Xkintra ðbj jbi Þ provides a transition probability from bin bi to bj (Fig. SF2b), where bi is i-th bin in v k . To compute Xðv l jv k Þ from the virtual-state sequence, a countmatrix ðCÞ was defined, where Cðv k ; v l Þ denotes transition counts from v k to v l . If the simulation is long enough, an equality stands theoretically: Cðv k ; v l Þ ¼ Cðv l ; v k Þ. It is known that the use of multiple parallel runs increases sampling statistics [16,30,31].

B. Dasgupta et al. / Chemical Physics Letters 662 (2016) 327–332

However, if the sampling is inadequate, an inequality Cðv k ; v l Þ – Cðv l ; v k Þ is resulted. Such inequality (observed frequently in MD) means that the detailed balance may not be maintained accurately [24]. To impose the equality on the count matrix, the Markov state model uses the maximum likelihood method [24,32]. Note that imposition of equality to the count matrix is symmetrization of the matrix. We adapted this method from MSMBuilder2 [23,24,32] and converted it to R-script with necessary modifications. Likewise, from the subtrajectories we computed a count-matrix element C kintra ðbi ; bj Þ, which is the number of counts from bi to bj within virtual state v k , and obtained symmetrized count matrix. Finally, the transition probability matrices

X and Xkintra were calculated from the symmetrized count matrices. By using MSM, equilibrium population for virtual states  fw1 k jk ¼ 1; 2; . . . ; mg is estimated by applying X on an initial population ðfw1 k jk ¼ 1; 2; . . . ; mgÞ of virtual states for infinite times.
 Similarly, an equilibrium population c1 ik for the bin bi in v k is cal


culated by applying Xkintra for infinite times on the initial counts ðcik Þ. A probability distribution can be calculated ordinarily by using those initial counts without using Markov state modeling (see supplementary material S1, equation S1.15). Practically, we and obtained the eigenperformed diagonalization of X and X vectors assigned to the largest eigenvalue (i.e., 1). Theoretically, these eigenvectors are equivalent to the equilibrium values for 1 c1 ik and wk [24]. Using these equilibrium values we define probability in bin bi in v k as k intra

, p1 ik

¼

1 w1 k c ik

k

Nb m X X 1 w1 k c ik k

ð3Þ

i

where p1 ik is a Markov approximated probability in bin bi of

v k and

the quantity N kb is the total number of bins in v k . It is expected that k 1 the distribution PMA VAUS ðk; v k ; TÞ ¼ fpik ji ¼ 1; 2; . . . ; N b g is more reasonable (supplementary Fig. SF3) than P VAUS ðk; v k ; TÞ because detailed balance was imposed on the count matrixes. Different virtual states may have different k widths. Then, we 1 converted w1 k to wk =ðkk;max  kk;min Þ for normalization, where the denominator is the window size for v k .

329

differentiating Eq. (2) with respect to k (equation S1.14 in supplementary material S1). We note that such force update might lead to an abrupt change of biasing force between successive iterations and that this type of force change cannot be reduced by any of the aforementioned methods. Therefore, to reduce the effect of the abrupt change of biasing force we introduced a scaling factor ðrÞ empirically and write update equation for c > 1 as, ðcþ1Þ Dk ln P c ðk; v k Þðcþ1Þ ¼ Dk ln Pc ðk; v k ÞðcÞ þ rDk ln PMA ; VAUS ðk; v k Þ

ð4Þ where Dk ¼ d=dk and ‘T’ is dropped for simplicity. We set r by manual inspection of biasing forces (usually r ¼ 0:5). This empirical scaling factor can be interpreted as indirect pseudo-counts introduced in the probability distribution with assumption that if a longer simulation is performed, undersampled k bins are counted more than the observed counts, and over-sampled k bins are counted less than the observed ones. 3. Materials and methods We used two Ala-pentapeptides (Ace-(Ala)5-Nme) to study their dimerization using the VAUS protocol and the AUS protocol (see below). The reaction-coordinate was based on the intermolecular distance between the two peptides (Fig. 1). The k range to be sampled was set from 5 to 20 Å. Details for the peptide system and the MD procedure are given in supplementary material S2 [33–42]. In VAUS, we performed iterative refinement of bias using the improvements of VAUS described above with seven virtual states ðm ¼ 7Þ (supplementary Table ST1, Fig. SF1b). However, it was not clear how sampling influenced the distribution edges of each virtual state. To avoid such unclear artifact (if present), we fused those seven virtual states to a single virtual state ðm ¼ 1Þ after iteration 9 ðc ¼ 9Þ. For iterations c > 9, we refer to observed probabil-

ity distribution as P MA VAUS ðk; v 1 ; TÞ. Note that the single virtual-state VAUS is equivalent to conventional AUS. We referred this protocol as VAUS protocol, although the later part of sampling was done with a single virtual state.

2.2.2. Improved parameterization of the observed probability distribution During virtual-state transition ðv l ! v k Þ, P c ðk; v l ; TÞ changes to P c ðk; v k ; TÞ with fixed k value [22]. Here, we note that such switching may cause an abrupt change of bias if P c ðk; v l ; TÞ – P c ðk; v k ; TÞ in the overlapping region of v k and v l . Such a mismatch does not strictly maintain detailed balance between virtual states. Although ðcÞ P c ðk; v k ; TÞðcþ1Þ is computed from PMA (see below), a small VAUS ðk; v k ; TÞ but abrupt change frequently observed at the virtual-state transitions. We introduce the procedure below to decrease this artifact.

The distributions lnP VAUS ðk; v k ; TÞðcÞ for different v k calculated for iteration c are fitted by polynomials that may not match smoothly in the overlapping region. Therefore, the polynomials are refined by an iterative method. In the first iteration of polynomial fitting, MA

an initial polynomial is fitted to lnP VAUS ðk; v k ; TÞðcÞ by taking into account the distributions of surrounding virtual states (v k1 , for instance). The obtained polynomials are refit again using the surrounding distributions, and new polynomials are obtained (second iteration). This procedure is repeated till the obtained polynomials overlap smoothly (supplementary Fig. SF4). MA

2.2.3. De-sensitizing biasing forces MD requires forces on each atom. In the current method, biasing forces include dPc ðk; v k ; TÞ=dk, which is obtained by

Fig. 2. Plot of standard deviation (sd) of the observed mean logarithmic-probability ðcÞ (ln P MA or ln P AUS ðk; TÞðcÞ ) for each iteration ðcÞ from VAUS (empty VAUS ðk; v k ; TÞ circles) and AUS (filled circles) (supplementary material S2). Solid and dotted lines are linear fit to sd plots for VAUS and AUS, respectively (supplementary Table ST3). Dashed gray horizontal line indicates sd = 0.1. Gray vertical dashed line indicates the last multi-virtual state iteration ðc ¼ 9Þ for VAUS.

330

B. Dasgupta et al. / Chemical Physics Letters 662 (2016) 327–332

To analyze the efficiency of VAUS we performed AUS simulation (single virtual-state VAUS) as a control, which is referred here as AUS protocol. In AUS, we did not use the methods proposed under the Improvements of VAUS in Theory to analyze the efficiency of VAUS. Here we refer to observed probability distribution as PAUS ðk; v 1 ; TÞ or simply P AUS ðk; TÞ. For both VAUS and AUS, the identical number of iterations was performed with identical simulation lengths (supplementary Table ST2). The probability distributions for iteration 1 were set as Pc ðk; v k ; TÞð1Þ ¼ 1, for all k for VAUS and Pc ðk; v 1 ; TÞð1Þ ¼ 1 for AUS (single virtual state). The initial conformations were those obtained from iteration 0 (see supplementary material S2). Total cumulative sampling was 3.4 ls to refine the biasing force, and the production run was 320 ns. Furthermore, we performed conventional canonical simulations of the system (supplementary material S2).

4. Results 4.1. Convergence of simulation We sampled the Ala-pentapeptide system by accurately biasing k, and obtained nearly uniform P MA VAUS ðk; v k ; TÞ from iteration 26 for VAUS and P AUS ðk; TÞ for AUS (supplementary Fig. SF5). The performance of Markovian method was verified by comparing observed probability distributions estimated by using and not using MSM for c ¼ 1. We observed that in c ¼ 1 by using Markov approximation observed probability distribution can be extrapolated by at least 30% along time than probability distribution without using MSM (see supplementary Fig. SF3 and corresponding legend for calculation details). Therefore the distribution obtained under Markov approximation was more acceptable. An estimate of the uniformity of the distributions was calculated from standard deviation (sd) of the mean logarithmicprobability. From visual inspection we observed that an sd < 0.1 resulted in approximately uniform distribution (supplementary Fig. SF5). The plot of sd with iteration step revealed how the convergence was achieved (Fig. 2). VAUS protocol gave sd < 0.1 for the last four iterations unlike AUS protocol. A linear fit to the sd revealed that VAUS protocol led to a faster convergence than AUS protocol. Recall that in VAUS, multiple virtual states were used until 9th iteration. From c ¼ 10 onwards single virtual state was used in the VAUS protocol. It is likely that these first iterative refinements helped for better convergence rate. 4.2. PMF Before the simulation converge, one can calculate PMF (defined as ‘RT ln P c ðk; TÞ’) [19] for each iterative simulation ðPMFðcÞÞ from

ln Pc ðk; v k ; TÞðcÞ . To analyze how PMF changed over iterations we compared PMF from the production run ðc ¼ 26Þ to other PMFsðc < 26Þ (i.e., root-mean-square deviation (RMSD) between

Fig. 3. Root-mean-square deviation (RMSD) between ln P c ðk; TÞð26Þ and ln P c ðk; TÞðcÞ as a function of c, where ln P c ðk; TÞð26Þ is canonical distribution from the final production runs of VAUS and AUS (supplementary material S2).

ln Pc ðk; v k ; TÞð26Þ and ln P c ðk; v k ; TÞðcÞ ) in Fig. 3. We did not include the term ‘RT’ in the RMSD calculation to highlight the difference of probability distribution. Overall, RMSDs from VAUS were lower than those from AUS, indicating that VAUS provided a better convergence to the final distribution than AUS did. Moreover, in VAUS,

Fig. 4. Comparison of ensemble features between VAUS and AUS protocols. (a) Frequency distribution of number of inter-peptide chain hydrogen bonds [43]. (b) Cumulative proportion of variances (see supplementary Table ST4 legend) from the principal component analysis (PCA) of subsets of the complete ensemble obtained under VAUS and AUS (Table ST4). Black lines indicate PCA of a subset defined by 5.0 Å 6 k 6 5.5 Å, while gray lines indicate PCA of complete ensemble (5.0 Å 6 k 6 20 Å).

B. Dasgupta et al. / Chemical Physics Letters 662 (2016) 327–332

PMFðc ¼ 8—10Þ were also similar to PMFðc ¼ 26Þ, indicating final PMF was approximated well in the early iterations. The final canonical distributions ln P c ðk; TÞ were compared between the two protocols (supplementary Fig. SF6), and we found that in both protocols the free-energy minimum (the most stable

331

state) was near 8 Å. Moreover, the curves agreed well between 8 and 20 Å. However, they did not match well below 8 Å; conformations from VAUS were more stable than those from AUS. We could not judge which distribution is qualitatively more reasonable for k < 8 Å only from Fig. SF6. We discuss this point later. We observed several conformational differences in the final conformational ensemble between the two protocols. (1) We observed more inter-chain hydrogen bonds [43] in VAUS than in AUS (Fig. 4a), most remarkably at the four inter-chain hydrogen bonds. Doubtlessly, those inter-chain hydrogen bonds helped to stabilize bound structures and to decrease PMF in k ¼ 5—8 Å. (2) We generated different subsets of configurations existing in different k-ranges (e.g. a subset of configurations in k range of 5.0–5.5 Å) and performed principle component analysis (PCA) to these subsets. PCA provided a way to assess conformational variability. Fig. 4b and supplementary Table ST4 show that VAUS had a greater conformational variability than AUS for the ensemble with k-range 5.0–5.5 Å. For an ensemble with wider k range, the difference between the two protocols became insignificant. Moreover, we monitored intra-chain hydrogen bonds leading to secondary structures and turn formation [44–46] as a function of k in supplementary Fig. SF7. Previously 310-helix and a-helix had been observed experimentally in Ala-rich peptides [47]. Fig. SF7 shows that inter-strand b-sheets were formed at lower k-region, whereas intra-chain secondary structures were formed as the two peptides were separated more. 4.3. 2D free-energy landscape identifying major conformational states To understand detailed characteristics of conformational ensembles, we computed two-dimensional (2D) free-energy landscapes at 300 K (Fig. 5a and b). Here, the two dimensional space were constructed by the reaction-coordinate k and an angle ðhÞ representing molecular orientation (Fig. 1). A comparison of landscapes between VAUS and AUS revealed that while a stable conformational state was found in the small k region (5.0–5.5 Å) for VAUS (Fig. 5a), almost no conformational states stood out for AUS

Fig. 5. 2D free-energy landscape PMFðk; hÞ using final production run of VAUS and AUS ((a and b), respectively) and (c) conventional canonical MD. Color scale is shown in the top right panel. In each panel, the lowest-PMF site is set as PMF = 0 (supplementary material S2). The black boxes in (a) highlight two major basins obtained from VAUS. In (a), structures picked from the marked cells, which are represented by white (relatively stable) and black (relatively unstable) numbers, are used to initiate conventional canonical simulations. The landscape from AUS (b) is shallower than that from (a). (c) The landscape from the conventional canonical MD is computed from the resultant trajectories.

Fig. 6. Three representative structures from basins in Fig. 5a are shown. Hydrogen bonds are shown by dotted lines. Residue numbers are indicated to show direction of chains. Traces of chains are shown in cartoon representation.

332

B. Dasgupta et al. / Chemical Physics Letters 662 (2016) 327–332

(Fig. 5b). Thus VAUS and AUS provided different distributions in this small k region (Fig. SF6). Some representative structures taken from the basin in the small k region of Fig. 5a were illustrated in Fig. 6. To further clarify the 2D-energy landscape of VAUS, we picked 11 relatively stable conformations from basins of the landscape given in Fig. 5a and 19 random conformations from other regions. These conformations were used for the initial conformations of unbiased conventional canonical simulations (2 ns long for each run). Fig. 5c is the resultant 2D free-energy landscape, which showed better correlation to the landscape from VAUS (correlation coefficient = 0.75) than that from AUS (0.61). This result clarified qualitative accuracy of the landscape given in Fig. 5a; the stable conformations sampled from VAUS were also stable in the conventional MD simulations. A few canonical trajectories were shown in the supplementary media (SM). 5. Conclusion We have introduced a set of novel techniques for VAUS to simulate molecular binding: Markov-approximated probability distribution calculation, iterative polynomial fitting, and desensitizing biasing forces by a pseudo-counting approach. The VAUS method was used to adaptively update canonical distribution for some initial iterations followed by the regular AUS. We compared the VAUS protocol (multi-virtual state VAUS followed by single-virtual state VAUS) to the AUS protocol (single-virtual state VAUS without improvements of VAUS) (Figs. 2 and 3), which revealed some distinctive characteristic features of canonical ensemble (Fig. 5). The difference in the canonical ensemble between the two protocols may be due to the initial VAUS iterations, which used multiplevirtual states. This helps to cover wider conformational space around bound states (Fig. 4b), strongly bound states (Figs. 4a, 6 and SF7), and unbound states. The validity of the virtual-system coupling was conceptual in a previous study [2] because a toy model was used and sampling was based on a Monte-Carlo scheme. Assessment of VAUS in the present study was more strict and practical than the previous study. We believe that the current VAUS protocol is effective for larger biomolecular system to understand binding phenomena or dynamics of ligand around the active site of a receptor [48]. Acknowledgement The authors thank Dr. Ikuo Fukuda, Dr. Akira Kinjo and Gert-Jan Bekker for helpful discussions. B. D. was supported by JBIC, Japan. J. H. was supported by JSPS KAKENHI Grant No. 16K05517, the Development of core technologies for innovative drug development based upon IT from Japan Agency for Medical Research and Development, AMED, and by the HPCI System Research Project (Project IDs: hp150015). Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.cplett.2016.09. 059.

References [1] R. Anandakrishnan, A. Drozdetski, R.C. Walker, A.V. Onufriev, Biophys. J. 108 (2015) 1153. [2] J. Higo, J. Ikebe, N. Kamiya, H. Nakamura, Biophys. Rev. 4 (2012) 27. [3] A. Mitsutake, Y. Sugita, Y. Okamoto, Biopolymers 60 (2001) 96. [4] N. Nakajima, H. Nakamura, A. Kidera, J. Phys. Chem. B 101 (1997) 817. [5] N. Nakajima, J. Higo, A. Kidera, H. Nakamura, J. Mol. Biol. 296 (2000) 197. [6] N. Nakajima, Chem. Phys. Lett. 288 (1998) 319. [7] J. Higo, N. Nakajima, H. Shirai, A. Kidera, H. Nakamura, J. Comput. Chem. 18 (1997) 2086. [8] J. Ikebe, S. Sakuraba, H. Kono, J. Comput. Chem. 35 (2014) 39. [9] M. Mezei, J. Comput. Phys. 68 (1987) 237. [10] G.H. Paine, H.A. Scheraga, Biopolymers 24 (1985) 1391. [11] G.M. Torrie, J.P. Valleau, J. Comput. Phys. 23 (1977) 187. [12] Y. Deng, B. Roux, J. Phys. Chem. B 113 (2009) 2234. [13] J.L. Adelman, M. Grabe, J. Chem. Phys. 138 (2013) 44105. [14] R.J. Allen, D. Frenkel, P.R. ten Wolde, J. Chem. Phys. 124 (2006) 24102. [15] H. Alonso, A.A. Bliznyuk, J.E. Gready, Med. Res. Rev. 26 (2006) 531. [16] J. Higo, B. Dasgupta, T. Mashimo, K. Kasahara, Y. Fukunishi, H. Nakamura, J. Comput. Chem. 36 (2015) 1489. [17] B.A. Berg, T. Neuhaus, Phys. Rev. Lett. 68 (1992) 9. [18] U.H.E. Hansmann, Y. Okamoto, J. Comput. Chem. 14 (1993) 13. [19] J. Higo, K. Umezawa, H. Nakamura, J. Chem. Phys. 138 (2013) 184106. [20] J. Higo, O.V. Galzitskaya, S. Ono, H. Nakamura, Chem. Phys. Lett. 337 (2001) 169. [21] C. Bartels, M. Karplus, J. Phys. Chem. B 102 (1998) 865. [22] S. Iida, H. Nakamura, J. Higo, Biochem. J. 473 (2016) 1651. [23] F. Noé, S. Fischer, Curr. Opin. Struct. Biol. 18 (2008) 154. [24] J.-H. Prinz, H. Wu, M. Sarich, B. Keller, M. Senne, M. Held, J.D. Chodera, C. Schütte, F. Noé, J. Chem. Phys. 134 (2011) 174105. [25] J.D. Chodera, F. Noé, Curr. Opin. Struct. Biol. 25 (2014) 135. [26] D. Shukla, C.X. Hernández, J.K. Weber, V.S. Pande, Acc. Chem. Res. 48 (2015) 414. [27] T.J. Lane, G.R. Bowman, K. Beauchamp, V.A. Voelz, V.S. Pande, J. Am. Chem. Soc. 133 (2011) 18413. [28] G.R. Bowman, D.L. Ensign, V.S. Pande, J. Chem. Theory Comput. 6 (2010) 787. [29] J.D. Chodera, W.C. Swope, J.W. Pitera, K.A. Dill, Multiscale Model. Simul. 5 (2006) 1214. [30] J. Ikebe, K. Umezawa, N. Kamiya, T. Sugihara, Y. Yonezawa, Y. Takano, H. Nakamura, J. Higo, J. Comput. Chem. 32 (2011) 1286. [31] H. Junichi, K. Narutoshi, S. Takanori, Y. Yasushige, N. Haruki, Chem. Phys. Lett. 473 (2009) 326. [32] K.A. Beauchamp, G.R. Bowman, T.J. Lane, L. Maibaum, I.S. Haque, V.S. Pande, J. Chem. Theory Comput. 7 (2011) 3412. [33] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. [34] D.J. Evans, G.P. Morriss, Phys. Lett. A 98 (1983) 433. [35] N. Kamiya, Y.S. Watanabe, S. Ono, J. Higo, Chem. Phys. Lett. 401 (2005) 312. [36] W.D. Cornell, P. Cieplak, C.I. Bayly, I.R. Gould, K.M. Merz, D.M. Ferguson, D.C. Spellmeyer, T. Fox, J.W. Caldwell, P.A. Kollman, J. Am. Chem. Soc. 117 (1995) 5179. [37] P. Kollman, R. Dixon, W. Cornell, T. Fox, C. Chipot, A. Pohorille, in: W.F. van Gunsteren, P.K. Weiner, A.J. Wilkinson (Eds.), Computer Simulation of Biomolecular Systems: Theoretical and Experimental Applications, vol. 27, Springer Science & Business Media, 2013. [38] J. Ryckaert, G. Ciccotti, H. Berendsen, J. Comput. Phys. 23 (1977) 327. [39] I. Fukuda, N. Kamiya, Y. Yonezawa, H. Nakamura, J. Chem. Phys. 137 (2012) 54314. [40] T. Mashimo, Y. Fukunishi, N. Kamiya, Y. Takano, I. Fukuda, H. Nakamura, J. Chem. Theory Comput. 9 (2013) 5599. [41] W. Humphrey, A. Dalke, K. Schulten, J. Mol. Graph. 14 (1996) 33. [42] K. Morikami, T. Nakai, A. Kidera, M. Saito, H. Nakamura, Comput. Chem. 16 (1992) 243. [43] I.K. McDonald, J.M. Thornton, J. Mol. Biol. 238 (1994) 777. [44] B. Dasgupta, L. Pal, G. Basu, P. Chakrabarti, Proteins 55 (2004) 305. [45] W. Kabsch, C. Sander, Biopolymers 22 (1983) 2577. [46] L. Pal, P. Chakrabarti, G. Basu, J. Mol. Biol. 326 (2003) 273. [47] G.L. Millhauser, C.J. Stenland, P. Hanson, K.A. Bolin, F.J. van de Ven, J. Mol. Biol. 267 (1997) 963. [48] W. Wang, O. Donini, C.M. Reyes, P.A. Kollman, Annu. Rev. Biophys. Biomol. Struct. 30 (2001) 211.