Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations

The Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) and the Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) methods calculate binding free energies for macromolecules by combining molecular mechanics calculations and continuum solvation models. To systematically evaluate the performance of these methods, we report here an extensive study of 59 ligands interacting with six different proteins. First, we explored the effects of the length of the molecular dynamics (MD) simulation, ranging from 400 to 4800 ps, and the solute dielectric constant (1, 2, or 4) on the binding free energies predicted by MM/PBSA. The following three important conclusions could be observed: (1) MD simulation length has an obvious impact on the predictions, and longer MD simulation is not always necessary to achieve better predictions. (2) The predictions are quite sensitive to the solute dielectric constant, and this parameter should be carefully determined according to the characteristics of the protein/ligand binding interface. (3) Conformational entropy often show large fluctuations in MD trajectories, and a large number of snapshots are necessary to achieve stable predictions. Next, we evaluated the accuracy of the binding free energies calculated by three Generalized Born (GB) models. We found that the GB model developed by Onufriev and Case was the most successful model in ranking the binding affinities of the studied inhibitors. Finally, we evaluated the performance of MM/GBSA and MM/PBSA in predicting binding free energies. Our results showed that MM/PBSA performed better in calculating absolute, but not necessarily relative, binding free energies than MM/GBSA. Considering its computational efficiency, MM/GBSA can serve as a powerful tool in drug design, where correct ranking of inhibitors is often emphasized.     

Fast and accurate predictions of binding free energies using MM‐PBSA and MM‐GBSA

In the drug discovery process, accurate methods of computing the affinity of small molecules with a biological target are strongly needed. This is particularly true for molecular docking and virtual screening methods, which use approximated scoring functions and struggle in estimating binding energies in correlation with experimental values. Among the various methods, MM‐PBSA and MM‐GBSA are emerging as useful and effective approaches. Although these methods are typically applied to large collections of equilibrated structures of protein‐ligand complexes sampled during molecular dynamics in water, the possibility to reliably estimate ligand affinity using a single energy‐minimized structure and implicit solvation models has not been explored in sufficient detail. Herein, we thoroughly investigate this hypothesis by comparing different methods for the generation of protein‐ligand complexes and diverse methods for free energy prediction for their ability to correlate with experimental values. The methods were tested on a series of structurally diverse inhibitors of Plasmodium falciparum DHFR with known binding mode and measured affinities. The results showed that correlations between MM‐PBSA or MM‐GBSA binding free energies with experimental affinities were in most cases excellent. Importantly, we found that correlations obtained with the use of a single protein‐ligand minimized structure and with implicit solvation models were similar to those obtained after averaging over multiple MD snapshots with explicit water molecules, with consequent save of computing time without loss of accuracy. When applied to a virtual screening experiment, such an approach proved to discriminate between true binders and decoy molecules and yielded significantly better enrichment curves. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Crystal molecular dynamics simulations to speed up MM/PB(GB)SA evaluation of binding free energies of di-mannose deoxy analogs with P51G-m4-Cyanovirin-N

Complexes of two Cyanovirin-N (CVN) mutants, m4-CVN and P51G-m4-CVN, with deoxy di-mannose analogs were employed as models to generate conformational ensembles using explicit water Molecular Dynamics (MD) simulations in solution and in crystal environment. The results were utilized for evaluation of binding free energies with the molecular mechanics Poisson-Boltzmann (or Generalized Born) surface area, MM/PB(GB)SA, methods. The calculations provided the ranking of deoxy di-mannose ligands affinity in agreement with available qualitative experimental evidences. This confirms the importance of the hydrogen-bond network between di-mannose 3'- and 4'-hydroxyl groups and the protein binding site B(M) as a basis of the CVN activity as an effective HIV fusion inhibitor. Comparison of binding free energies averaged over snapshots from the solution and crystal simulations showed high promises in the use of the crystal matrix for acceleration of the conformational ensemble generation, the most time consuming step in MM/PB(GB)SA approach. Correlation between energy values based on solution versus crystal ensembles is 0.95 for both MM/PBSA and MM/GBSA methods.     

新型二氟甲基磷酸类酪氨酸蛋白磷酸酯酶1B抑制剂的分子动力学模拟和结合自由能计算

通过分子对接建立了一系列含二氟甲基磷酸基团(DFMP)或二氟甲基硫酸基团(DFMS)的抑制剂与酪氨酸蛋白磷酸酯酶1B(PTP1B)的相互作用模式, 并通过1 ns的分子动力学模拟和molecular mechanics/generalized Born surface area (MM/GBSA)方法计算了其结合自由能. 计算获得的结合自由能排序和抑制剂与靶酶间结合能力排序一致; 通过基于主方程的自由能计算方法, 获得了抑制剂与靶酶残基间相互作用的信息, 这些信息显示DFMP/DFMS基团的负电荷中心与PTP1B的221位精氨酸正电荷中心之间的静电相互作用强弱决定了此类抑制剂的活性, 进一步的分析还显示位于DFMP/DFMS基团中的氟原子或其他具有适当原子半径的氢键供体原子会增进此类抑制剂与PTP1B活性位点的结合能力.     

Importance of accurate charges in molecular docking: quantum mechanical/molecular mechanical (QM/MM) approach

The extent to which accuracy of electric charges plays a role in protein-ligand docking is investigated through development of a docking algorithm, which incorporates quantum mechanical/molecular mechanical (QM/MM) calculations. In this algorithm, fixed charges of ligands obtained from force field parameterization are replaced by QM/MM calculations in the protein environment, treating only the ligands as the quantum region. The algorithm is tested on a set of 40 cocrystallized structures taken from the Protein Data Bank (PDB) and provides strong evidence that use of nonfixed charges is important. An algorithm, dubbed "Survival of the Fittest" (SOF) algorithm, is implemented to incorporate QM/MM charge calculations without any prior knowledge of native structures of the complexes. Using an iterative protocol, this algorithm is able in many cases to converge to a nativelike structure in systems where redocking of the ligand using a standard fixed charge force field exhibits nontrivial errors. The results demonstrate that polarization effects can play a significant role in determining the structures of protein-ligand complexes, and provide a promising start towards the development of more accurate docking methods for lead optimization applications.     

Binding affinities of factor Xa inhibitors estimated by thermodynamic integration and MM/GBSA

We present free energy estimates of nine 3-amidinobenzyl-1H-indole-2-carboxamide inhibitors of factor Xa. Using alchemical thermodynamic integration (TI) calculations, we estimate the difference in binding free energies with high accuracy and precision, except for mutations involving one of the amidinobenzyl rings. Crystal studies show that the inhibitors may bind in two distinct conformations, and using TI, we show that the two conformations give a similar binding affinity. Furthermore, we show that we can reduce the computational demand, while still retaining a high accuracy and precision, by using fewer integration points and shorter protein-ligand simulations. Finally, we have compared the TI results to those obtained with the simpler MM/GBSA method (molecular-mechanics with generalized Born surface-area solvation). MM/GBSA gives better results for the mutations that involve a change of net charge, but if a precision similar to that of the TI method is required, the MM/GBSA method is actually slightly more expensive. Thus, we have shown that TI could be a valuable tool in drug design.     

基于分子动力学模拟和连续介质模型的自由能计算方法*

近些年,基于分子动力学模拟和连续介质模型的自由能计算方法受到了越来越多的关注,其中MM/PBSA就是最具代表性的方法.在MM/PBSA中,体系的焓变采用分子力学(MM)的方法计算得到;溶剂效应中极性部分对自由能的贡献通过解Poisson-Boltzmann(PB)方程的方法计算得到;溶液效应中非极性部分对自由能的贡献则通过分子表面积(SA)计算得到.本文结合我们科研组的工作,就近几年MM/PBSA方法的最新进展做了较为详细的阐述,同时对MM/PBSA的发展前景进行了展望.     

Further development and validation of empirical scoring functions for structure-based binding affinity prediction

New empirical scoring functions have been developed to estimate the binding affinity of a given protein-ligand complex with known three-dimensional structure. These scoring functions include terms accounting for van der Waals interaction, hydrogen bonding, deformation penalty, and hydrophobic effect. A special feature is that three different algorithms have been implemented to calculate the hydrophobic effect term, which results in three parallel scoring functions. All three scoring functions are calibrated through multivariate regression analysis of a set of 200 protein-ligand complexes and they reproduce the binding free energies of the entire training set with standard deviations of 2.2 kcal/mol, 2.1 kcal/mol, and 2.0 kcal/mol, respectively. These three scoring functions are further combined into a consensus scoring function, X-CSCORE. When tested on an independent set of 30 protein-ligand complexes, X-CSCORE is able to predict their binding free energies with a standard deviation of 2.2 kcal/mol. The potential application of X-CSCORE to molecular docking is also investigated. Our results show that this consensus scoring function improves the docking accuracy considerably when compared to the conventional force field computation used for molecular docking.     

An improved method to predict the entropy term with the MM/PBSA approach

A method is suggested to calculate improved entropies within the MM/PBSA approach (molecular mechanics combined with Poisson-Boltzmann and surface area calculations) to estimate protein-ligand binding affinities. In the conventional approach, the protein is truncated outside ~8 A from the ligand. This system is freely minimised using a distance-dependent dielectric constant (to simulate the removed protein and solvent). However, this can lead to extensive changes in the molecular geometry, giving rise to a large standard deviation in this term. In our new approach, we introduce a buffer region approximately 4 A outside the truncated protein (including solvent molecules) and keep it fixed during the minimisation. Thereby, we reduce the standard deviation by a factor of 2-4, ensuring that the entropy term no longer limits the precision of the MM/PBSA predictions. The new method is tested for the binding of seven biotin analogues to avidin, eight amidinobenzyl-indole-carboxamide inhibitors to factor Xa, and two substrates to cytochrome P450 3A4 and 2C9. It is shown that it gives more stable results and often improved predictions of the relative binding affinities.     

EGFR和4-苯胺喹唑啉类抑制剂之间相互作用模式的研究

采用分子动力学和MM/PBSA相结合的方法预测了表皮生长因子受体和4-苯胺喹 啉类抑制剂的相互作用模式。在分子动力学采样的基础上,采用MM/PBSA的方法分 别预测了四种可能结合模式下表皮生长因子受体和4-苯胺喹唑啉类抑制剂间的结合 自由能。在MM/PBSA计算中,受体和抑制剂之间的非键相互作用能采用分子力学 (MM)的方法得到;溶剂效应中极性部分对自由能的贡献通过解Possion- Boltzmanne (PB)方程的方法得到;溶液效应中非极性部分对自由能的贡献则通过 分子表面积计算(SA)的方法得到。计算表明,在四种结合模式下,表皮生长因子受 体和4-苯胺喹唑啉类抑制剂之间的结合自由能有较大的差别。在最佳的相互作用模 式中,抑制剂的苯胺部分位于活性口袋的底部,能够与受体残基的非极性侧链产生 很强的范德华和疏水相互作用。抑制剂喹唑啉环上的N(1)原子能够和Met-769上的 NH形成稳定的氢键,而抑制剂上的N(3)原子则和周围的一个水分子形成氢键。同时 ,抑制剂双环上的取代基团也能和活性口袋外部的部分残基形成一定的范德华和疏 水相互作用。最佳结合模式能够很好地解释已有抑制剂结构和活性间的关系。     

Solvent interaction energy calculations on molecular dynamics trajectories: increasing the efficiency using systematic frame selection

End-point methods such as linear interaction energy (LIE) analysis, molecular mechanics generalized Born solvent-accessible surface (MM/GBSA), and solvent interaction energy (SIE) analysis have become popular techniques to calculate the free energy associated with protein-ligand binding. Such methods typically use molecular dynamics (MD) simulations to generate an ensemble of protein structures that encompasses the bound and unbound states. The energy evaluation method (LIE, MM/GBSA, or SIE) is subsequently used to calculate the energy of each member of the ensemble, thus providing an estimate of the average free energy difference between the bound and unbound states. The workflow requiring both MD simulation and energy calculation for each frame and each trajectory proves to be computationally expensive. In an attempt to reduce the high computational cost associated with end-point methods, we study several methods by which frames may be intelligently selected from the MD simulation including clustering and address the question of how the number of selected frames influences the accuracy of the SIE calculations.     

Protein/ligand binding free energies calculated with quantum mechanics/molecular mechanics

The calculation of binding affinities for flexible ligands has hitherto required the availability of reliable molecular mechanics parameters for the ligands, a restriction that can in principle be lifted by using a mixed quantum mechanics/molecular mechanics (QM/MM) representation in which the ligand is treated quantum mechanically. The feasibility of this approach is evaluated here, combining QM/MM with the Poisson-Boltzmann/surface area model of continuum solvation and testing the method on a set of 47 benzamidine derivatives binding to trypsin. The experimental range of the absolute binding energy (DeltaG = -3.9 to -7.6 kcal/mol) is reproduced well, with a root-mean-square (RMS) error of 1.2 kcal/mol. When QM/MM is applied without reoptimization to the very different ligands of FK506 binding protein the RMS error is only 0.7 kcal/mol. The results show that QM/MM is a promising new avenue for automated docking and scoring of flexible ligands. Suggestions are made for further improvements in accuracy.     

Alpha7 nicotinic acetylcholine receptor agonists: prediction of their binding affinity through a molecular mechanics Poisson-Boltzmann surface area approach

A group of agonists for the alpha7 neuronal nicotinic acetylcholine receptors (nAChRs) was investigated, and their free energies of binding DeltaG(bind) were calculated by applying the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) approach. This method, based on molecular dynamics simulations of fully solvated protein-ligand complexes, allowed us to estimate the contribution of both polar and nonpolar terms as well as the entropy to the overall free energy of binding. The calculated results were in a good agreement with the experimentally determined DeltaG(bind) values, thereby pointing to the MM-PBSA protocol as a valuable computational tool for the rational design of specific agents targeting the neuronal alpha7 nAChR subtypes.     

Ligand-specific scoring functions: improved ranking of docking solutions

Molecular docking is a powerful computational method that has been widely used in many biomolecular studies to predict geometry of a protein-ligand complex. However, while its conformational search algorithms are usually able to generate correct conformation of a ligand in the binding site, the scoring methods often fail to discriminate it among many false variants. We propose to treat this problem by applying more precise ligand-specific scoring filters to re-rank docking solutions. In this way specific features of interactions between protein and different types of compounds can be implicitly taken into account. New scoring functions were constructed including hydrogen bonds, hydrophobic and hydrophilic complementarity terms. These scoring functions also discriminate ligands by the size of the molecule, the total hydrophobicity, and the number of peptide bonds for peptide ligands. Weighting coefficients of the scoring functions were adjusted using a training set of 60 protein-ligand complexes. The proposed method was then tested on the results of docking obtained for an additional 70 complexes. In both cases the success rate was 5-8% better compared to the standard functions implemented in popular docking software.     

A comparison of different initialization protocols to obtain statistically independent molecular dynamics simulations

We study how the results of molecular dynamics (MD) simulations are affected by various choices during the setup, e.g., the starting velocities, the solvation, the location of protons, the conformation of His, Asn, and Gln residues, the protonation and titration of His residues, and the treatment of alternative conformations. We estimate the binding affinity of ligands to four proteins calculated with the MM/GBSA method (molecular mechanics combined with a generalized Born and surface area solvation energy). For avidin and T4 lysozyme, all variations gave similar results within 2 kJ/mol. For factor Xa, differences in the solvation or in the selection of alternative conformations gave results that are significantly different from those of the other approaches by 4-6 kJ/mol, whereas for galectin-3, changes in the conformations, rotations, and protonation gave results that differed by 10 kJ/mol, but only if residues close to the binding site were modified. This shows that the results of MM/GBSA calculations are reasonably reproducible even if the MD simulations are set up with different software. Moreover, we show that the sampling of phase space can be enhanced by solvating the systems with different equilibrated water boxes, in addition to the common use of different starting velocities. If different conformations are available in the crystal structure, they should also be employed to enhance the sampling. Protonation, ionization, and conformations of Asn, Gln, and His may also be used to enhance sampling, but great effort should be spent to obtain as reliable predictions as possible close to the active site.     

How to obtain statistically converged MM/GBSA results

The molecular mechanics/generalized Born surface area (MM/GBSA) method has been investigated with the aim of achieving a statistical precision of 1 kJ/mol for the results. We studied the binding of seven biotin analogues to avidin, taking advantage of the fact that the protein is a tetramer with four independent binding sites, which should give the same estimated binding affinities. We show that it is not enough to use a single long simulation (10 ns), because the standard error of such a calculation underestimates the difference between the four binding sites. Instead, it is better to run several independent simulations and average the results. With such an approach, we obtain the same results for the four binding sites, and any desired precision can be obtained by running a proper number of simulations. We discuss how the simulations should be performed to optimize the use of computer time. The correlation time between the MM/GBSA energies is ～5 ps and an equilibration time of 100 ps is needed. For MM/GBSA, we recommend a sampling time of 20–200 ps for each separate simulation, depending on the protein. With 200 ps production time, 5–50 separate simulations are required to reach a statistical precision of 1 kJ/mol (800–8000 energy calculations or 1.5–15 ns total simulation time per ligand) for the seven avidin ligands. This is an order of magnitude more than what is normally used, but such a number of simulations is needed to obtain statistically valid results for the MM/GBSA method. © 2009 Wiley Periodicals, Inc. J Comput Chem 2010     

Analyzing the robustness of the MM/PBSA free energy calculation method: application to DNA conformational transitions

The ability to predict and characterize free energy differences associated with conformational equilibria or the binding of biomolecules is vital to understanding the molecular basis of many important biological functions. As biological studies focus on larger molecular complexes and properties of the genome, proteome, and interactome, the development and characterization of efficient methods for calculating free energy becomes increasingly essential. The aim of this study is to examine the robustness of the end-point free energy method termed the molecular mechanics Poisson-Boltzmann solvent accessible surface area (MM/PBSA) method. Specifically, applications of MM/PBSA to the conformational equilibria of nucleic acid (NA) systems are explored. This is achieved by comparing A to B form DNA conformational free energy differences calculated using MM/PBSA with corresponding free energy differences determined with a more rigorous and time-consuming umbrella sampling algorithm. In addition, the robustness of NA MM/PBSA calculations is also evaluated in terms of the sensitivity towards the choice of force field and the choice of solvent model used during conformational sampling. MM/PBSA calculations of the free energy difference between A-form and B-form DNA are shown to be in very close agreement with the PMF result determined using an umbrella sampling approach. Further, it is found that the MM/PBSA conformational free energy differences were also in agreement using either the CHARMM or AMBER force field. The influence of ionic strength on conformational stability was particularly insensitive to the choice of force field. Finally, it is also shown that the use of a generalized Born implicit solvent during conformational sampling results in free energy estimates that deviate slightly from those obtained using explicitly solvated MD simulations in these NA systems.     

Computational Characterization of Binding of Small Molecule Inhibitors to HIV‐1 gp41

Developing orally available small molecule inhibitors of HIV‐1 fusion has attracted significant interest over many years. Frey had recently reported several synthetic compounds which are experimentally shown to inhibit cell‐cell fusion in the low micromolar range. We carried out computational study to help identify possible binding modes by docking these compounds onto the hydrophobic pocket on gp41 and to characterize structures of binding complexes. The detailed gp41‐molecule binding interactions and free energies of binding are obtained through molecular dynamics simulation and MM‐PBSA calculation. Specific molecular interactions in the gp41‐inhibitor complexes are identified. The present computational study complements the corresponding experimental investigation and helps establish a good starting point for further refinement of small molecular gp41 inhibitors.     

Extension of QM/MM docking and its applications to metalloproteins

To overcome the limitation of conventional docking methods which assume fixed charge model from force field parameters, combined quantum mechanics/molecular mechanics (QM/MM) method has been applied to docking as a variable charge model and shown to exhibit improvement on the docking accuracy over fixed charge based methods. However, it has also been shown that there are a number of examples for which adoption of variable‐charge model fails to reproduce the native binding modes. In particular, for metalloproteins, previously implemented method of QM/MM docking failed most often. This class of proteins has highly polarized binding sites at which high‐coordinate‐numbered metal ions reside. We extend the QM/MM docking method so that protein atoms surrounding the binding site along with metal ions are included as quantum region, as opposed to only ligand atoms. This extension facilitates the required scaling of partial charges on metal ions leading to prediction of correct binding modes in metalloproteins. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009