The article information
 Pengfei Li, Xiangyu Jia, Meiting Wang, Ye Mei
 李鹏飞, 贾相瑜, 王美婷, 梅晔
 Comparison of Accuracy and Convergence Rate between Equilibrium and Nonequilibrium Alchemical Transformations for Calculation of Relative Binding Free Energy^{†}
 平衡和非平衡方法计算相对结合自由能的精度和效率比较
 Chinese Journal of Chemical Physics, 2017, 30(6): 789799
 化学物理学报, 2017, 30(6): 789799
 http://dx.doi.org/10.1063/16740068/30/cjcp1711204

Article history
 Received on: November 7, 2017
 Accepted on: December 27, 2017
b. NYUECNU Center for Computational Chemistry at NYU Shanghai, Shanghai 200062, China;
c. Department of Chemistry and Biochemistry, University of Oklahoma, Norman OK 73019, USA
The decrease of free energy is the driving force of any chemical or biological processes [17]. Estimation of free energy within the chemical accuracy is critical for quantitative understanding of biomolecular functions. To this end, much effort has been devoted to the development of new methods for the predictions of proteinligand binding affinity over the past several decades. However, there are still many grand challenges in structurebased rational drug design [812]. Much of the difficulty lies in the accuracy of the employed force fields, the attainable spatiotemporal scale of computer simulations, and the reliability of the postprocessing methods. There are an arsenal of methods with different features for the prediction of proteinligand binding affinity, including docking and scoring [13, 14], the linear interaction energy (LIE) [1517] and linear response approximation (LRA) [18] methods, molecular mechanics/PoissonBoltzmann (generalized Born) surface area (MM/PB(GB)SA) simulations [1921], and thermodynamic perturbation (TP) [22] based free energy methods. The endpoint approximation methods (e.g., LIE, LRA, MMPBSA/GBSA) compute ''bound'' and the ''free'' states. These methods often overestimate the enthalpic benefit upon ligand binding (by inadequate accounting for water movement) and underestimate the entropic penalty (due to undersampling of the relevant conformers) [23]. Recently, alchemical free energy simulations have been widely used to predict the relative binding free energy [2429], in which a nonphysical thermodynamic path is constructed to connect the two thermodynamic states in question. Due to that sufficient phase space overlap between the two states is the precondition for an accurate free energy estimation, a series of immediate states joining the two states with a slow varying coupling parameter
All the equilibrium (EQ) free energy estimators, such as TI, TP, Bennett acceptance ratio (BAR) [32, 33] and its multiple state variant (Multistate BAR or MBAR) [34, 35], assume that equilibrium has been reached in each molecular dynamics simulation and the sampled conformations are uncorrelated. Rather than spending large computational costs on long equilibrium samplings, in 1997, Jarzynski proved that the equilibrium free energy difference can be extracted from a series of nonequilibrium (NE) trajectories, which has now become the wellknown ''Jarzynski equality'' [36]. In our previous study [37], this extraordinary equality has been used to ''correct'' the solvation free energies of thirteen amino acid side chains analogs from a molecular mechanical (MM) level to a quantum mechanical/molecular mechanical (QM/MM) level. It has been found that when the phase space overlap between the two Hamiltonians is insufficient, the free energy difference between the MM Hamiltonian and the QM/MM Hamiltonian estimated by TP method becomes obscure. In contrast, averaging over many short nonequilibrium pulling simulations using the Jarzynski equality can yield significantly improved results. The performance of nonequilibrium MD simulations is determined by the transformation velocity between the thermodynamic states and the number of nonequilibrium simulation trajectories. With fixed total computational time, these two factors conflict with each other, and there must be a systemdependent balance between them that leads to the highest efficiency. However, in contrast to the free energy calculations based on equilibrium trajectories, nonequilibrium free energy methods are still in their infancy. Their performance in all aspects, for instance their performance in the calculations of relative binding free energy studied in this work, needs thorough investigations.
Calculating the binding free energy difference as the difference of absolute binding free energies (the vertical legs shown in FIG. 1) is not costeffective, and quite often it may encounter numerical difficulties [3844]. Instead, it can be calculated as the difference of two alchemical mutation free energies in the binding site and in the aqueous solution for the ligands [4547]. In this work, the equilibrium and nonequilibrium simulations for the relative binding free energy calculations are compared in terms of their accuracy and convergence rate. Two wellstudied systems were chosen. One is T4lysozyme mutant L99A (L99A) as shown in FIG. 2. It is a good model system for the comparison of free energy methods [42, 43, 48, 49], due to its rigidity upon binding of benzene and phenol. Another receptor is cyclooxygenase2 (COX2), of which the structure is also shown in FIG. 2. This system is relatively more challenging [5053]. Two ligands bound to COX2 were studied. One ligand has a hydroxyl group, and the other has a methyl substituent.
Ⅱ. METHODThe relative binding free energy
$ \begin{equation} \Delta \Delta G_{\rm{bind}} = \Delta G_{\rm{bind}}({\rm{B}})\Delta G_{\rm{bind}}(\rm{A}) \end{equation} $  (1) 
where
$ \begin{equation} \Delta \Delta G_{\rm{bind}} = \Delta G_{\rm{pert}}({\rm{bound}})\Delta G_{\rm{pert}}(\rm{free}) \end{equation} $  (2) 
using the thermodynamic cycle depicted in FIG. 1. The process in the upper rung with a free energy change
Within these stages, a series of intermediate states characterized by
$ \begin{equation} \Delta G=\frac{1}{\beta}\ln \frac{\left< f(U_{i}U_{i+1}+C)\right>_{i+1}}{\left< f(U_{i+1}U_{i}C)\right>_{i}}+C \label{Eq:BAR1} \end{equation} $  (3) 
and
$ \begin{equation} C=\Delta G+\frac{1}{\beta}\ln \frac{N_{i+1}}{N_{i}} \end{equation} $  (4) 
where
$ \begin{eqnarray} \delta^{2}\hspace{0.15cm}&=&\hspace{0.15cm}\frac{1}{\beta^{2}N_{i}}\left[\frac{\langle f^{2}(U_{i+1}U_{i}C)\rangle_{i}}{\langle f(U_{i+1}U_{i}C) \rangle^{2}_{i}}1\right]+ \nonumber \\ &&\frac{1}{\beta^{2}N_{i+1}}\left[\frac{\langle f^{2}(U_{i}U_{i+1}+C)\rangle_{i+1}}{\langle f(U_{i}U_{i+1}+C)\rangle_{i+1}^{2}}1\right] \end{eqnarray} $  (5) 
In 1997, Jarzynski proved for the first time that the equilibrium free energy difference
$ \begin{equation} \exp(\beta \Delta G)=\left<\exp(\beta W)\right> \end{equation} $  (6) 
Here, the angled bracket
$ \begin{equation} \delta W(t\rightarrow t+\delta t)=U({{\bf{r}}}(t), \lambda(t+ \delta t))U({{\bf{r}}}(t), \lambda(t)) \end{equation} $  (7) 
where
$ \begin{equation} W=\sum\limits_{t=0}^{\tau  \delta t}U\left({\bf{r}}(t), \lambda(t+ \delta t)\right)U\left({\bf{r}}(t), \lambda(t)\right) \end{equation} $  (8) 
The variance of
$ \begin{equation} \sigma^{2}=\frac{1}{\beta^2N}\left(\frac{\sigma_{x}}{\left<x\right>}\right)^2, \, \, \, {\rm{with}} \, \, x \equiv \exp\left(\beta W\right) \end{equation} $  (9) 
where
$ \begin{eqnarray} \sum\limits_{i=1}^{n_{\rm{F}}}&&\hspace{0.5cm}\frac{1}{1+\exp \left[\beta(M+W_{i}\Delta G)\right]} \nonumber\\ &=&\hspace{0.15cm}\sum\limits_{j=1}^{n_{\rm{R}}}\frac{1}{1+\exp \left[\beta(M+W_{j}\Delta G)\right]} \end{eqnarray} $  (10) 
where
For T4lysozyme mutant L99A (L99A), two ligands (benzene and phenol) are studied. The initial structure of L99A is obtained from protein data bank (PDB ID: 181L). For cyclooxygenase2 (COX2), two celecoxib analogues are studied as the ligands. The initial structure of COX2 is also obtained from protein data bank (PDB ID: 1CX2). The initial structures of the ligands were optimized at B3LYP/631G
$ \begin{eqnarray} V_{\rm{SCvdW}}&=&4\epsilon_{ij}(1\lambda)\cdot\nonumber\\ &&\left\{\frac{1}{\left[\alpha\lambda+\left(\displaystyle{\frac{r_{ij}}{\sigma_{ij}}}\right)^6\right]^2} \frac{1}{\alpha\lambda+\left(\displaystyle{\frac{r_{ij}}{\sigma_{ij}}}\right)^6}\right\} \end{eqnarray} $  (11) 
where
In equilibrium simulations, at each stage, eleven alchemical windows (
In nonequilibrium simulations, at each stage, one hundred forward transformation measurements and one hundred backward transformation measurements were combined to calculate the free energy changes. In each nonequilibrium trajectory,
All the MD simulations were carried out with the pmemd program in AMBER 14 package [88], and all the QM calculations were carried out with Gaussian 09 package [89].
Ⅲ. RESULTS AND DISCUSSION A. Relative binding free energies for L99A lysozyme systemAs shown in FIG. 4, convergence of free energy calculations by the equilibrium and nonequilibrium methods is satisfactory. In FIG. 4(a), the calculated free energy difference at each stage can reach converged results in 5 ns for each window. Moreover, the standard deviations at each stage are in a range of 0.010.05 kcal/mol, which indicates that the precision of the calculated results is quite satisfactory. In fact, the results calculated in a 2 ns simulation for each window have already converged as can be seen in this figure. The calculated free energy difference at each stage, depicted as the black lines in FIG. 4(b), can also reach convergence with 100 nonequilibrium simulation trajectories. The standard deviations at each stage range in 0.010.07 kcal/mol, which have a similar precision to the equilibrium simulations. The free energy changes calculated by Jarzynski equality from the forward and backward nonequilibrium transformations are plotted in red and blue, respectively. It is well known to the community that unidirectional methods such as thermodynamic perturbation and the Jarzynski equality overestimate the free energy difference between the target Hamiltonian and the sampled Hamiltonian, mainly due to that in most cases the important region of the target Hamiltonian has not been wellsampled [5760, 77]. Therefore, usually the free energy difference (
Therefore, all the data listed in Table Ⅰ are converged and reliable. As listed in Table Ⅰ, the free energy changes associated with the three transformation stages of the ''bound'' state calculated by the equilibrium simulations are 7.30, 4.57, and 34.34 kcal/mol, respectively. And those from the nonequilibrium simulations are 7.33, 4.53, and 34.33 kcal/mol, respectively. In the ''free'' state, the free energy changes are slightly weaker. The free energy changes in the neutralization, alchemy and charging stages are 6.45, 1.39, 33.85 kcal/mol from the equilibrium simulations, and those are 6.47, 1.36, and 33.92 kcal/mol from the nonequilibrium simulations. At all three stages, the free energy differences calculated by the nonequilibrium simulations are very close to those from the equilibrium simulations, indicating that the nonequilibrium simulations can yield statistically the same results as the equilibrium simulations.
The relative binding free energy for benzene and phenol binding into L99A system are (1.84
Listed in Table Ⅱ are the cost of wall clock time at each stage. All the calculations were done on servers containing two Intel Xeon CPU E52660 2.20 GHz CPU with 16 cores in total. For the equilibrium simulations, the wall clock time for the neutralization, alchemy and charging stages of the ''bound'' state are 97, 97, and 96 h, respectively. Before the nonequilibrium simulations, an equilibrium simulation is required to generate an ensemble of structures of the initial state as the starting structures for the nonequilibrium simulations. Therefore, the wall clock time used for the equilibrium simulation should also be taken into account as the cost of the nonequilibrium simulations. In the ''bound'' state, the initial equilibrium simulations took 3 h each. The total wall clock time of the neutralization, alchemy and charging stages of the ''bound'' state are 10 (3+7), 7 (3+4), and 7 (3+4) h, which are far less than that of the equilibrium simulation (shown in the parentheses are the wall clock time for the initial equilibrium simulation and the actual nonequilibrium simulation, respectively). In addition, the wall clock time of the neutralization, alchemy and charging stages of the ''free'' state respectively calculated by the equilibrium simulations are 5 h each, and those for the nonequilibrium simulations are 1.3 (0.4+0.9), 1.2 (0.4+0.8) and 0.9 (0.4+0.5) h, respectively. It shows again the nonequilibrium simulation methods are far more efficient than the equilibrium methods. The total time for the equilibrium and the nonequilibrium simulations are 305 and 27 h, respectively. The nonequilibrium simulation method yielded consistent results with comparable variances as the equilibrium simulation method using only about 9 percent of the wall clock time. In addition, while the explicit TIP3P water model was used in this work, implicit solvent models can be used to further enhance the efficiency of nonequilibrium sampling due to the reduction of dissipation as shown in a recent work of Liu et al. [90].
B. Relative binding free energies for COX2 systemAs shown in FIG. 5, convergence of the free energy calculations by the equilibrium and nonequilibrium method is also satisfactory. In FIG. 5(a), the calculated free energy difference at each stage can reach converged results in 5 ns for each window. Moreover, the standard deviations at each stage range in 0.010.04 kcal/mol, which also indicates that the precision of free energy calculations is quite satisfactory. The calculated free energy difference at each stage, depicted as the black line in FIG. 5(b), can also reach convergence with 100 nonequilibrium simulation trajectories. The standard deviations at each stage range in 0.010.06 kcal/mol, which also has the similar precision as done by equilibrium simulations. The free energy changes calculated by Jarzynski equality from the forward and backward nonequilibrium transformations are plotted in red and blue, respectively. It can also been seen that all the results done by unidirectional nonequilibrium simulations can converge to the results from the analysis using the bidirectional BAR method, which combines forward and backward directions. The bidirectional BAR results indeed converge faster. Therefore, all the data listed in Table Ⅲ are reliable. The free energy changes associated with the three transformation stages of the ''bound'' state calculated by equilibrium simulations are 40.57, 11.49, and 27.75 kcal/mol, respectively. And those from the nonequilibrium simulations are 40.08, 11.09, and 27.80 kcal/mol. In the ''free'' state, the free energy changes in the neutralization, alchemy and charging stages are 39.61, 9.59, and 26.69 kcal/mol from the equilibrium simulations, and are 39.58, 9.59, and 26.72 kcal/mol from the nonequilibrium simulations. At each stage, the free energy difference calculated by the nonequilibrium simulation is also very close to that from the equilibrium simulation, which also indicates that the nonequilibrium simulations can yield statistically the same results as the equilibrium simulations. The relative binding free energy for ligand ''OH'' and ligand ''CH
Listed in Table Ⅳ are the wall clock time for each stage of the transformation using the same computers as above. For the equilibrium simulations, the wall clock time for the neutralization, alchemy and charging stages of the ''bound'' state respectively are all 235 h. For the nonequilibrium simulations, the wall clock time are 100 (45+55), 83 (37+46), and 17 (8+9) h, which are also far cheaper than the equilibrium simulations. For the equilibrium simulations of the ''free'' state, the wall clock time are around 17 h to 18 h for the neutralization, alchemy and charging transformations. Using the nonequilibrium simulations, the computer time reduces to 1.5 (0.5+1), 1.5 (0.5+1) and 2.3 (1.7+0.6) h. The total time for equilibrium and nonequilibrium simulations are 757 and 205 h, respectively, and the nonequilibrium simulation method saved about 73 percent of the wall clock time as compared to the equilibrium simulation method to reach the same uncertainties for the results.
Ⅳ. CONCLUSIONComputeraided rational drug design is still facing the gap in spatiotemporal resolution between the experimental measurement and that can be afforded by computer simulations. With decades of improvement, equilibrium simulationsbased methods have become the standard technique for the calculations of free energies. However, these methods are still too expensive, for instance when used in highthroughput screening. With the discovery of the Jarzynski equality, it becomes more and more clearly that equilibrium simulations methods are nothing but special cases of nonequilibirum methods, and the latter gradually attracts wide attention from the community. In this work, we applied both the equilibrium and the nonequilibrium simulation methods to the calculations of the relative binding free energies of two ligands to L99A Lysozyme and two ligands to COX2. The results show that the nonequilibrium molecular dynamics simulations can yield statistically the same results as their equilibrium counterparts in the calculations of the relative binding free energies, but with higher efficiency in term of the wall clock time. The nonequilibrium simulation method saved about 91 percent and 73 percent of wall clock time as compared to the equilibrium method for L99A Lysozyme and COX2, respectively. Moreover, the bidirectional analysis of the NE trajectories with BAR converges much faster than the unidirectional Jarzynski equality.
Ⅴ. ACKNOWLEDGMENTSThis work is supported by the National Natural Science Foundation of China (No.21773066). CPU time was provided by the Supercomputer Center of East China Normal University.
[1]  P. A. Bash, U. C. Singh, R. Langridge, and P. A. Kollman, Science 236 , 564 (1987). DOI:10.1126/science.3576184 
[2]  W. L. Jorgensen, Acc. Chem. Res. 22 , 184 (1989). DOI:10.1021/ar00161a004 
[3]  P. Kollman, Chem. Rev. 93 , 2395 (1993). DOI:10.1021/cr00023a004 
[4]  A. Pohorille, C. Jarzynski, and C. Chipot, J. Phys. Chem. B 114 , 10235 (2010). DOI:10.1021/jp102971x 
[5]  C. D. Christ, and A. E. Mark, W. F. van Gunsteren, J. Comput. Chem. 31 , 1569 (2010). 
[6]  N. Hansen, W. F. van Gunsteren, J. Chem. Theory Comput. 10 , 2632 (2014). DOI:10.1021/ct500161f 
[7]  C. Chipot, A. Pohorille, and E ds, Free Energy Calculations:Theory and Applications in Chemistry and Biology, Springer Series in Chemical Physics, Vol. 86, Verlag Berlin Heidelberg:Springer (2007). 
[8]  W. L. Jorgensen, Acc. Chem. Res. 42 , 724 (2009). DOI:10.1021/ar800236t 
[9]  W. L. Jorgensen, Science 303 , 1813 (2004). DOI:10.1126/science.1096361 
[10]  L. Wang, Y. Wu, Y. Deng, B. Kim, L. Pierce, G. Krilov, D. Lupyan, S. Robinson, M. K. Dahlgren, J. Greenwood, D. L. Romero, C. Masse, J. L. Knight, T. Steinbrecher, T. Beuming, W. Damm, E. Harder, W. Sherman, M. Brewer, R. Wester, M. Murcko, L. Frye, R. Farid, T. Lin, D. L. Mobley, W. L. Jorgensen, B. J. Berne, R. A. Friesner, and R. Abel, J. Am. Chem. Soc. 137 , 2695 (2015). DOI:10.1021/ja512751q 
[11]  D. J. Huggins, W. Sherman, and B. Tidor, J. Med. Chem. 55 , 1424 (2012). DOI:10.1021/jm2010332 
[12]  C. Bissantz, B. Kuhn, and M. Stahl, J. Med. Chem. 53 , 5061 (2010). DOI:10.1021/jm100112j 
[13]  E. Yuriev, M. Agostino, and P. A. Ramsland, J. Mol. Recognit. 24 , 149 (2011). DOI:10.1002/jmr.v24.2 
[14]  E. Yuriev, and P. A. Ramsland, J. Mol. Recognit. 26 , 215 (2013). DOI:10.1002/jmr.v26.5 
[15]  J. Aqvist, C. Medina, and J. Samuelsson, Protein Eng. 7 , 385 (1994). DOI:10.1093/protein/7.3.385 
[16]  J. Aqvist, and J. Marelius, Comb. Chem. High Throughput Screening 4 , 613 (2001). DOI:10.2174/1386207013330661 
[17]  X. Jia, J. Zeng, J. Z. H. Zhang, and Y. Mei, J. Comput. Chem. 35 , 737 (2014). DOI:10.1002/jcc.v35.9 
[18]  N. Makri, J. Phys. Chem. B 103 , 2823 (1999). DOI:10.1021/jp9847540 
[19]  J. Srinivasan, T. E. Cheatham, P. Cieplak, P. A. Kollman, and D. A. Case, J. Am. Chem. Soc. 120 , 9401 (1998). DOI:10.1021/ja981844+ 
[20]  N. Homeyer, and H. Gohlke, Mol. Inf. 31 , 114 (2012). DOI:10.1002/minf.v31.2 
[21]  T. Hou, J. Wang, Y. Li, and W. Wang, J. Chem. Inf. Model. 51 , 69 (2011). DOI:10.1021/ci100275a 
[22]  R. W. Zwanzig, J. Chem. Phys. 22 , 1420 (1954). DOI:10.1063/1.1740409 
[23]  B. J. WilliamsNoonan, E. Yuriev, and D. K. Chalmers, J. Med. Chem (2017). DOI:10.1021/acs.jmedchem.7b00681(2017) 
[24]  S. Bruckner, and S. Boresch, J. Comput. Chem. 32 , 1303 (2011). DOI:10.1002/jcc.v32.7 
[25]  P. V. Klimovich, M. R. Shirts, and D. L. Mobley, J. Comput. Aided Mol. Des. 29 , 397 (2015). DOI:10.1007/s1082201598409 
[26]  M. R. Shirts, and V. S. Pande, J. Chem. Phys. 122 , 114107 (2005). DOI:10.1063/1.1861872 
[27]  F. M. Ytreberg, R. H. Swendsen, and D. M. Zuckerman, J. Chem. Phys. 125 , 184114 (2006). DOI:10.1063/1.2378907 
[28]  D. Rodriguez, E. GomezDarve, and A. Pohorille, J. Chem. Phys. 120 , 3563 (2004). DOI:10.1063/1.1642607 
[29]  H. Paliwal, and M. R. Shirts, J. Chem. Theory Comput. 7 , 4115 (2011). DOI:10.1021/ct2003995 
[30]  J. G. Kirkwood, J. Chem. Phys. 3 , 300 (1935). DOI:10.1063/1.1749657 
[31]  N. Homeyer, F. Stoll, A. Hillisch, and H. Gohlke, J. Chem. Theory Comput. 10 , 3331 (2014). DOI:10.1021/ct5000296 
[32]  C. H. Bennett, J. Comput. Phys. 22 , 245 (1976). DOI:10.1016/00219991(76)900784 
[33]  A. Hahn, and H. Then, Phys. Rev. E 80 , 031111 (2009). DOI:10.1103/PhysRevE.80.031111 
[34]  M. R. Shirts, and J. D. Chodera, J. Chem. Phys. 129 , 124105 (2008). DOI:10.1063/1.2978177 
[35]  M. R. Shirts, Reweighting from the Mixture Distribution as a Better Way to Describe the Multistate Bennett Acceptance Ratio, arXiv.org (2017). 
[36]  C. Jarzynski, Phys. Rev. Lett. 78 , 2690 (1997). DOI:10.1103/PhysRevLett.78.2690 
[37]  M. Wang, P. Li, X. Jia, W. Liu, Y. Shao, W. Hu, J. Zheng, B. R. Brooks, and Y. Mei, J. Chem. Inf. Model. 57 , 2476 (2017). DOI:10.1021/acs.jcim.7b00001 
[38]  D. Hamelberg, and J. A. McCammon, J. Am. Chem. Soc. 126 , 7683 (2004). DOI:10.1021/ja0377908 
[39]  S. Boresch, F. Tettinger, M. Leitgeb, and M. Karplus, J. Phys. Chem. B 107 , 9535 (2003). DOI:10.1021/jp0217839 
[40]  J. Wang, Y. Deng, and B. Roux, Biophys. J. 91 , 2798 (2006). DOI:10.1529/biophysj.106.084301 
[41]  M. Aldeghi, A. Heifetz, M. J. Bodkin, S. Knapp, and P. C. Biggin, Chem. Sci. 7 , 207 (2016). DOI:10.1039/C5SC02678D 
[42]  D. L. Mobley, A. P. Graves, J. D. Chodera, A. C. McReynolds, B. K. Shoichet, and K. A. Dill, J. Mol. Biol. 371 , 1118 (2007). DOI:10.1016/j.jmb.2007.06.002 
[43]  Y. Deng, and B. Roux, J. Phys. Chem. B. 113 , 2234 (2009). DOI:10.1021/jp807701h 
[44]  W. Liu, X. Jia, M. Wang, P. Li, X. Wang, W. Hu, J. Zheng, and Y. Mei, RSC Adv. 7 , 38570 (2017). DOI:10.1039/C7RA06215J 
[45]  M. R. Reddy, M. D. Varney, V. Kalish, V. N. Viswanadhan, and K. Appelt, J. Med. Chem. 37 , 1145 (1994). DOI:10.1021/jm00034a012 
[46]  M. R. Reddy, and M. D. Erion, J. Am. Chem. Soc. 123 , 6246 (2001). DOI:10.1021/ja0103288 
[47]  C. Oostenbrink, W. F. van Gunsteren, Proc. Natl. Acad. Sci. USA 102 , 6750 (2005). DOI:10.1073/pnas.0407404102 
[48]  A. Morton, W. A. Baase, and B. W. Matthews, Biochemistry 34 , 8564 (1995). DOI:10.1021/bi00027a006 
[49]  Y. Deng, and B. Roux, J. Chem. Theory Comput. 2 , 1255 (2006). DOI:10.1021/ct060037v 
[50]  T. D. Penning, J. J. Talley, S. R. Bertenshaw, J. S. Carter, P. W. Collins, S. Docter, M. J. Graneto, L. F. Lee, J. W. Malecha, J. M. Miyashiro, R. S. Rogers, D. J. Rogier, S. S. Yu, G. D. Anderson, E. G. Burton, J. N. Cogburn, S. A. Gregory, C. M. Koboldt, W. E. Perkins, K. Seibert, A. W. Veenhuizen, Y. Y. Zhang, and P. C. Isakson, J. Med. Chem. 40 , 1347 (1997). DOI:10.1021/jm960803q 
[51]  J. Michel, M. L. Verdonk, and J. W. Essex, J. Chem. Theory Comput. 3 , 1645 (2007). DOI:10.1021/ct700081t 
[52]  J. Michel, M. L. Verdonk, and J. W. Essex, J. Med. Chem. 49 , 7427 (2006). DOI:10.1021/jm061021s 
[53]  M. L. Plount Price, and W. L. Jorgensen, J. Am. Chem. Soc. 122 , 9455 (2000). DOI:10.1021/ja001018c 
[54]  X. Jia, M. Wang, Y. Shao, G. Konig, B. R. Brooks, J. Z. H. Zhang, and Y. Mei, J. Chem. Theory Comput. 12 , 499 (2016). DOI:10.1021/acs.jctc.5b00920 
[55]  C. Jarzynski, Phys. Rev. E 56 , 5018 (1997). DOI:10.1103/PhysRevE.56.5018 
[56]  G. Crooks, J. Statis. Phys. 90 , 1481 (1998). DOI:10.1023/A:1023208217925 
[57]  B. P. Cossins, S. Foucher, and C. M. Edge, J. W. Essex, J. Phys. Chem. B 113 , 5508 (2009). DOI:10.1021/jp803532z 
[58]  M. Goette, and H. Grubmuller, J. Comput. Chem. 30 , 447 (2009). DOI:10.1002/jcc.v30:3 
[59]  D. A. Hendrix, and C. Jarzynski, J. Chem. Phys. 114 , 5974 (2001). DOI:10.1063/1.1353552 
[60]  G. Hummer, J. Chem. Phys. 114 , 7330 (2001). DOI:10.1063/1.1363668 
[61]  F. M. Ytreberg, and D. M. Zuckerman, J. Chem. Phys. 120 , 10876 (2004). DOI:10.1063/1.1760511 
[62]  W. Lechner, H. Oberhofer, C. Dellago, and P. L. Geissler, J. Chem. Phys. 124 , 044113 (2006). DOI:10.1063/1.2162874 
[63]  P. Nicolini, D. Frezzato, and R. Chelli, J. Chem. Theory Comput. 7 , 582 (2011). DOI:10.1021/ct100568n 
[64]  H. Oberhofer, C. Dellago, and P. L. Geissler, J. Phys. Chem. B 109 , 6902 (2005). DOI:10.1021/jp044556a 
[65]  C. Oostenbrink, W. F. van Gunsteren, Chem. Phys. 323 , 102 (2006). 
[66]  M. Spichty, M. Cecchini, and M. Karplus, J. Phys. Chem. L 1 , 1922 (2010). DOI:10.1021/jz1005016 
[67]  P. S. Hudson, H. L. Woodcock, and S. Boresch, J. Phys. Chem. L 6 , 4850 (2015). DOI:10.1021/acs.jpclett.5b02164 
[68]  V. A. Ngo, I. Kim, T. W. Allen, and S. Y. Noskov, J. Chem. Theory Comput. 12 , 1000 (2016). DOI:10.1021/acs.jctc.5b01050 
[69]  I. Echeverria, and L. M. Amzel, J. Phys. Chem. B 116 , 10986 (2012). DOI:10.1021/jp300527q 
[70]  D. M. Zuckerman, and T. B. Woolf, Chem. Phys. Lett. 351 , 445 (2002). DOI:10.1016/S00092614(01)013975 
[71]  R. Chelli, S. Marsili, A. Barducci, and P. Procacci, Phys. Rev. E 75 , 050101 (2007). 
[72]  R. Chelli, S. Marsili, and P. Procacci, Phys. Rev. E 77 , 031104 (2008). DOI:10.1103/PhysRevE.77.031104 
[73]  R. Chelli, and P. Procacci, Phys. Chem. Chem. Phys. 11 , 1152 (2009). DOI:10.1039/b810914c 
[74]  P. Nicolini, P. Procacci, and R. Chelli, J. Phys. Chem. B 114 , 9546 (2010). DOI:10.1021/jp102263y 
[75]  P. Procacci, and C. Cardelli, J. Chem. Theory Comput. 10 , 2813 (2014). DOI:10.1021/ct500142c 
[76]  R. B. Sandberg, M. Banchelli, C. Guardiani, S. Menichetti, G. Caminati, and P. Procacci, J. Chem. Theory Comput. 11 , 423 (2015). DOI:10.1021/ct500964e 
[77]  C. Jarzynski, Phys. Rev. E 73 , 046105 (2006). DOI:10.1103/PhysRevE.73.046105 
[78]  G. Crooks, Phys. Rev. E 60 , 2721 (1999). DOI:10.1103/PhysRevE.60.2721 
[79]  G. Crooks, Phys. Rev. E 61 , 2361 (1999). 
[80]  M. R. Shirts, E. Bair, G. Hooker, and V. S. Pande, Phys. Rev. Lett. 91 , 140601 (2003). DOI:10.1103/PhysRevLett.91.140601 
[81]  C. I. Bayly, P. Cieplak, W. Cornell, and P. A. Kollman, J. Phys. Chem. 97 , 10269 (1993). DOI:10.1021/j100142a004 
[82]  P. Cieplak, W. D. Cornell, C. Bayly, and P. A. Kollman, J. Comput. Chem. 16 , 1357 (1995). DOI:10.1002/(ISSN)1096987X 
[83]  J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, and D. A. Case, J. Comput. Chem. 25 , 1157 (2004). DOI:10.1002/(ISSN)1096987X 
[84]  J. A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K. E. Hauser, and C. Simmerling, J. Chem. Theory Comput. 11 , 3696 (2015). DOI:10.1021/acs.jctc.5b00255 
[85]  R. W. Pastor, B. R. Brooks, and A. Szabo, Mol. Phys. 65 , 1409 (1988). DOI:10.1080/00268978800101881 
[86]  H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak, J. Chem. Phys. 81 , 3684 (1984). DOI:10.1063/1.448118 
[87]  T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98 , 10089 (1993). DOI:10.1063/1.464397 
[88]  D. A. Case, J. T. Berryman, R. M. Betz, D. S. Cerutti, T. E. Cheatham, Ⅲ, T. A. Darden, R. E. Duke, T. J. Giese, H. Gohlke, A. W. Goetz, N. Homeyer, N. Izadi, P. Janowski, J. Kaus, A. Kovalenko, T. S. Lee, S. LeGrand, P. Li, T. Luchko, R. Luo, B. Madej, K. M. Merz, G. Monard, P. Needham, H. Nguyen, H. T. Nguyen, I. Omelyan, A. Onufriev, D. R. Roe, A. Roitberg, R. SalomonFerrer, C. L. Simmerling, W. Smith, J. Swails, R. C. Walker, J. Wang, R. M. Wolf, X. Wu, D. M. York, and P. A. Kollman, AMBER 2014, San Francisco:University of California (2014). 
[89]  M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, Gaussian 09, Revision B.01, Wallingford, CT:Gaussian, Inc. (2010). 
[90]  H. Liu, F. Chen, H. Sun, D. Li, and T. Hou, J. Chem. Theory Comput. 13 , 1827 (2017). DOI:10.1021/acs.jctc.6b01139 
b. 华东师范大学上海纽约大学计算化学联合中心, 上海 200062;
c. 俄克拉荷马大学化学和生物化学系, 诺曼 73019