Methods Corner DOI: 10.1002/minf.201100135
Free Energy Calculations by the Molecular Mechanics Poisson Boltzmann Surface Area Method Nadine Homeyer[a] and Holger Gohlke*[a]
114
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Mol. Inf. 2012, 31, 114 – 122
Free Energy Calculations by the Molecular Mechanics Poisson
Boltzmann Surface Area Method
Abstract: Detailed knowledge of how molecules recognize interaction partners and of the conformational preferences of biomacromolecules is pivotal for understanding biochemical processes. Such knowledge also provides the foundation for the design of novel molecules, as undertaken in pharmaceutical research. Computer-based free energy calculations enable a detailed investigation of the energetic factors that are responsible for molecular stability or binding affinity. The Molecular Mechanics Poisson–Boltzmann Surface Area (MM-PBSA) approach is an efficient method for the calculation of free energies of diverse mo-
lecular systems. Here we describe the concepts of this approach and outline the practical proceeding. Furthermore we give an overview of the wide spectrum of problems that have been addressed with this method and of successful analyses carried out, thereby focussing on ambitious and recent studies. Limits of the approach in terms of accuracy and applicability are discussed. Despite these limitations MM-PBSA is a method with great potential that allows comparative free energy analyses for various molecular systems at low computational cost.
Keywords: MM-PBSA · Binding affinity · Implicit solvent · Molecular recognition · Drug design
1 Introduction Free energy calculations have become part of the standard repertoire of computational biochemical research due to their widespread applicability, reaching from the prediction of binding affinities of small drug-candidate compounds to the evaluation of relative stabilities of large biomacromolecular structures. The calculations can be performed by the rigorous and computationally expensive free energy perturbation (FEP) and thermodynamic integration (TI) methods, which have recently been described in this ”Methods Corner”.[1] The high computational costs of TI and FEP are caused, first, by the explicit treatment of solvent and, second, by determining the difference in free energies of two states based on simulations that are carried out at intermediate points along a transition path from one state to another. The computationally intensive evaluation of individual contributions of a large number of solvent molecules can be avoided by regarding solvent effects only implicitly, using continuum solvent methods. The computational cost can be further reduced by considering only the end-point states in the free energy calculations. End-point, implicit solvent approaches thus have the advantage of a relatively low computational cost, while important solvent effects are still regarded, so that they provide a trade-off between efficiency and accuracy. As a consequence, with today’s computational resources, the binding of large compound data sets to biomacromolecules can be analyzed within a reasonable time span, and these approaches are more accurate in binding pose and affinity predictions than simple scoring functions.[2] Probably the most well-known end-point, implicit solvent free energy method is the Molecular Mechanics Poisson– Boltzmann Surface Area (MM-PBSA) approach, for the first time so denoted in 1998.[3] The approach has been used in more than one hundred studies to determine the free energies of molecular systems. Topics addressed in these studies include evaluating docking poses, determining structural stability, and predicting binding affinities and hot-spots. This clearly demonstrates the wide applicability of the approach. As a further advantage, MM-PBSA allows analyzing Mol. Inf. 2012, 31, 114 – 122
individual energy contributions by means of a free energy decomposition, which gives additional insights into the energetics of the investigated system. Furthermore, MM-PBSA can be used to study very different molecular systems ranging from small chemical compounds binding to an RNA aptamer to large protein complexes with thousands of atoms.[4,5] In this review we provide an account of selected topics that have been addressed with the MM-PBSA approach, thereby focusing on more recent studies, and imparting knowledge to the reader about planning and evaluating MM-PBSA analyses.
2 State-of-the-Art 2.1 Methodology and Practical Proceeding
In MM-PBSA the free energy of a molecule is estimated as the sum of its gas-phase energy, the solvation free energy, and a contribution due to the configurational entropy of the solute (see Frontispiece[6]). The gas-phase energy is approximated by the molecular mechanics energy of the molecule, determined from a force field comprising terms for bond, angle, and torsion energies as well as van der Waals and electrostatic interactions. In the calculation of the solvation free energy, two contributions are considered, a polar and a nonpolar one. For the polar contribution, the change in the free energy upon transfer of a charged molecule from gas-phase (modeled as a homogeneous medium with dielectric constant e = 1) to solvent (modeled as a homogeneous medium with e = 78 or 80) is estimated by an implicit solvent model.[7,8] In the initial MM-PBSA implementation and still widely used, the polar contribution is calcu[a] N. Homeyer, H. Gohlke Heinrich-Heine-Universitt, Institut fr Pharmazeutische und Medizinische Chemie Dsseldorf, Germany *e-mail: [email protected] Supporting Information for this article is available on the WWW under http://dx.doi.org/10.1002/minf.201100135.
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.molinf.com
115
Methods Corner
N. Homeyer, H. Gohlke
lated via a finite-difference solution of the Poisson–Boltzmann equation described in detail elsewhere.[9,10] Alternatively, an implicit solvent model based on Generalized Born (GB) theory can be used, which is a computationally more efficient approximation to Poisson theory. This then leads to the so-called MM-GBSA variant.[11,12] The nonpolar contribution to the solvation free energy primarily arises from the free energy required for cavity formation for the solute within the solvent. Most frequently, this contribution is estimated by a term proportional to the solvent accessible surface area of the solute. As a more appropriate alternative, the nonpolar contribution is computed as the sum of a disfavorable energy resulting from cavity formation and a favorable energy resulting from attractive interactions between solute and solvent molecules.[12,13] Finally, the configurational entropy of the solute is usually estimated using a rigid-rotor harmonic oscillator approximation, applying either normal mode analysis[3] or quasi-harmonic analysis.[14] From the so-computed free energies of interacting molecules and by applying a thermodynamic cycle a binding free energy can be calculated as the difference between the free energy of a complex and the sum of the free energies of its components.[15] According to Figure 1, this binding free energy is equivalent to the sum of energy and configurational entropy contributions associated with complex formation in the gas-phase and the difference in solvation free energies between the complex and the unbound molecules. Changes in the configurational entropy are frequently neglected if only the relative binding free energy
Nadine Homeyer is a postdoctoral research fellow in the Institute for Pharmaceutical and Medicinal Chemistry at the Heinrich-Heine-University, Dsseldorf. She received her PhD from the University Erlangen-Nuremberg (Germany) in 2008 for her work on the effects of phosphorylation on protein structure, dynamics, and interaction. Her research interests lie in molecular recognition, lead optimization, and the dynamics of biomacromolecules. Holger Gohlke is a Professor of Pharmaceutical and Medicinal Chemistry at the Heinrich-Heine-University, Dsseldorf. His research aims at understanding and predicting receptor-ligand interactions and the modulation of biological processes by pharmacologically relevant molecules. His group develops and applies methods at the interface of computational pharmaceutical and biophysical chemistry and molecular bioinformatics.
116
www.molinf.com
of similar ligands shall be analyzed; this then results in an effective binding energy. Although MM-PBSA calculations can be performed based on single structures, e.g. minimized structures[2,16,17] or crystal structures,[18] they are typically conducted based on a conformational ensemble generated by molecular dynamics (MD) simulations. In the latter case, each energy component is determined by averaging over the respective energy contributions from all conformers in the ensemble. That way, effects of conformational changes taking place in the molecular system are considered in the free energy calculation, and an estimate of the precision of the calculations is obtained. It is common to run MD simulations with explicit solvent for generating the conformational ensemble, because continuum solvent simulations have been shown to yield less accurate results.[19] Usually, the molecular system is solvated in a box of water, and counter ions are added to enforce neutrality. In binding free energy studies it is also possible to just add a water sphere centered on the ligand and move only those atoms located in the core of this sphere during the simulation,[20] which decreases the computational costs. The simulations are typically conducted with the same force field that is used for the MM-PBSA calculations in order to avoid inconsistencies,[19] and long-range electrostatic interactions are computed by the particle-mesh Ewald procedure.[3,21] After the simulations, structures for the conformational ensemble are extracted from the generated MD trajectories. Water molecules and counter ions are removed, and the free energy is calculated as described above. In practice, depending on the information one wishes to gain from the free energy calculations, three main types of MM-PBSA calculations can be distinguished. If the stability of two conformations of a molecule shall be compared, simulations and free energy calculations are conducted separately for both conformers, and the free energies are directly evaluated (Figure 2A). If a binding free energy shall be computed the difference between the free energies of the complex and its components is calculated either from a single trajectory of the complex (“single-trajectory approach”, Figure 2B) or from separate trajectories of complex, receptor, and ligand (“three-trajectory approach”, Figure 2C). Although the singe-trajectory approach neglects the conformational flexibility of the unbound structures, it is usually applied in all cases where no large structural changes upon binding are expected, because it gives less noisy results due to the cancellation of intramolecular contributions[20,22] and, therefore, allows performing MM-PBSA analyses based on shorter simulations. Usual simulation times range from a few picoseconds to several nanoseconds depending on system flexibility, size of the data set, and desired precision of the results, although longer simulation times are recommended for obtaining properly thermalized conformational ensembles.
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Mol. Inf. 2012, 31, 114 – 122
Free Energy Calculations by the Molecular Mechanics Poisson
Boltzmann Surface Area Method
Figure 1. Thermodynamic cycle for calculating a binding free energy. In MM-PBSA, the binding free energy results as the sum of the energy and configurational entropy associated with complex formation in the gas-phase and the difference in solvation free energies between the complex and the unbound molecules.
2.2 Application in Drug Design
The lower computational costs, compared to TI and FEP, together with a more sophisticated computation of the free energy components, compared to common scoring functions, makes MM-PBSA an attractive method for all stages of drug design. In early stages, large compound data bases are screened for promising drug candidates using docking programs and scoring functions to assess the binding capability of the molecules to a target receptor. However, simple scoring functions may neglect important energetic contributions, such as the solvation free energy, and, thus, may have difficulties in correctly predicting the binding pose and binding affinity of a molecule. A combination of docking, molecular mechanics calculations, and MM-PBSA Mol. Inf. 2012, 31, 114 – 122
proved to be more successful in the identification of the correct ligand binding pose than pure docking approaches in several cases.[2,17,23,24] As such, in an extensive performance study, Huo et al.[2] analyzed the performance of MM-PBSA to discriminate correct from false ligand binding poses in 98 complexes. Initially, ligands were docked to receptors with the Autodock[25] program. For the 100 best-scored docking poses as well as the experimentally determined binding pose of each complex the positions of all atoms within a sphere of 9 centered at the ligand were optimized by molecular mechanics minimization. The optimized structures were reevaluated using implicit solvent free energy calculations. Although outperformed by another implicit solvent ap-
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.molinf.com
117
Methods Corner
N. Homeyer, H. Gohlke
Figure 2. Main types of MM-PBSA calculations: Comparison of the stability of two conformations (A) as well as binding free energy calculations according to the single-trajectory approach (B) and the three-trajectory approach (C).
proach, MM-PBSA reached a success rate of 79.6 % in the re-evaluation of the structures, if the top five scored poses were regarded. While, until recently, the use of MM-PBSA in virtual screening has been limited to evaluating a few hundred docking solutions,[26] it can nowadays be applied for scoring thousands of compounds due to the great increase in computer power. As the potential of MM-PBSA in discriminating true binders from a much larger number of decoys has been demonstrated in initial high-throughput virtual screening studies,[27,28] we expect the method to be used more extensively in lead identification in the future. However, the general benefit of MM-PBSA in this field still requires further validation on more protein systems and ligand data sets.[29] When leads have been identified they need to be further improved in order to obtain drug-candidates with good pharmacokinetic (absorption, distribution, metabolism, excretion, and toxicity) features that, at the same time, show a high binding affinity and selectivity for the target. In this stage of drug design, it is essential to be able to determine the relative binding affinity of a small number of drug-candidates, usually less than 50, with high accuracy. The main challenge for the applied method lies in the small range of binding affinities studied, which is typically below 10 kcal mol 1. In initial studies, promising compound rankings were described based on MM-PBSA relative binding affinities for single receptor systems.[20,30,31] However, it has been recognized recently that a more extensive testing is needed to validate the general applicability of the approach for lead optimization. With the required computational power now at hand, extensive performance studies investigating the 118
www.molinf.com
accuracy of MM-PBSA for determining relative binding affinities for several different protein-ligand systems can be conducted. In one such study, Huo et al.[22] calculated the binding free energies of 59 ligands interacting with six different proteins and found good correlations with r2 > 0.5 for three out of the six protein systems studied. In addition, Yang et al.[32] recently described the results of MM-PBSA analyses in which they determined the binding affinities of a total of 156 ligands to seven different protein families separated into six groups. Here, correlations of r2 > 0.5 were found for three out of six groups. However, it needs to be considered that the correlations might have come out weaker due to the sequence variation within the protein families in this case. Overall the results of the two studies suggest that the goodness of the agreement between experimentally determined binding affinities and those computed with MMPBSA does not only depend on the binding affinity range of the analyzed ligand data set, but is also determined by the specific features of the protein-ligand interaction and the degree of structural similarity between the studied compounds. While for sets of ligands with similar structures usually satisfactory results were obtained, the performance varied for those with diverse chemical structures.[22] Furthermore, the number of polar or charged groups involved in ligand binding seems to have a high influence on the calculated binding free energies. For systems with a large number of charged residues in the binding pocket and, hence, where a strong change in the polarization occurs upon binding, using e > 1 for the “gas-phase state” was shown to yield better results.[22,32] Thus, carefully choosing e according to the expected degree of polarization in the
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Mol. Inf. 2012, 31, 114 – 122
Free Energy Calculations by the Molecular Mechanics Poisson
Boltzmann Surface Area Method
binding pocket and/or based on test calculations with evaluation data may prove beneficial. Finally, applications of MM-PBSA in lead optimization have shown that for ligand data sets with a binding affinity range of less than 3 kcal mol 1 only unsatisfactory correlations of experimentally determined and computed relative binding affinities can be expected (Table S1). Although exceptions to this observation have been described for individual data sets,[22] a statistical uncertainty of at least 1.1 kcal mol 1 mean absolute deviation in the solvation free energy relative to a method-weighted average[18] and a standard deviation of 0.7–1.0 kcal mol 1 in the MM-PBSA binding free energy, found for 200 snapshots from a 2 ns simulation,[19] make it difficult to obtain good relative binding free energy predictions for ligand data sets with small binding affinity ranges. For data sets with wider binding affinity ranges the quality of the prediction depends on the specific features of the binding interaction and on the similarity of the ligands. The best results are obtained for ligand series with high structural similarity and/or high hydrophobicity as well as with a uniform total charge. 2.3 Application to Biomacromolecules
Free energy studies on biomacromolecules can give important insights into the structures and mechanisms underlying biological processes. As such, MM-PBSA has been applied to study the stability of DNA,[3,33] RNA,[34] and protein[35] conformers in order to identify the native fold among a set of decoy structures or the fold that is most stable in a given environment. Recently, it has been demonstrated that MM-PBSA is also able to predict disulfide connectivities in proteins.[36] When small cysteine-rich proteins of the knottin family with different possible patterns of disulfide bridges and similar overall structure were analyzed with MM-PBSA, the connectivities proposed by previous experimental studies were identified as being energetically most favorable. This study showed that, in those cases where determining the correct cysteine pairing is difficult with X-ray crystallography and nuclear magnetic resonance spectroscopy, MM-PBSA free energy calculations can help to identify the most favorable pattern. However, it has to be considered that the energetically most stable pattern of disulfide bridges does not need to be the connectivity found in nature, where the oxidative folding pathway might be kinetically controlled. Binding free energy calculations of biomacromolecules can yield information about the stability of complexes and the binding mechanisms. For example, additional knowledge about DNA binding of Cys2His2 zinc finger proteins could be gained from an energetic analysis of the interaction between the Zif268 zinc finger protein and DNA.[37] Binding affinity calculations for one-finger, two-finger, and three-finger Zif268 protein/DNA complexes revealed that the binding affinity does not increase linearly with the number of zinc fingers; rather, it is mainly enhanced by the Mol. Inf. 2012, 31, 114 – 122
addition of the third zinc finger domain. Molecular dynamics simulations and detailed energy decomposition analyses showed that the non-linear cooperativity in Zif268/DNA binding is primarily caused by a tighter interaction between these components, enabled by an enlargement of the DNA major grove and an unwinding of the DNA helix upon attachment of the third zinc finger. The gain in polar interaction energy caused by the tight DNA contact of the three finger protein is much larger than the cost of the DNA conformational change, so that the binding affinity is highly favorable for the protein with three zinc finger domains. The specificity of protein-protein, protein-DNA, and protein-RNA interactions can be studied by modifying residues involved in complex formation and analyzing the differences in binding free energies caused by these modifications. Such modifications can involve interconversions of DNA or RNA bases,[37,38] transformations of phospho-amino acid residues,[39] and amino acid mutations.[38] As for the latter, it is a popular approach to investigate protein-protein interactions by mutating residues in the contact interface consecutively to alanine and, then, calculating the binding free energy for wildtype and mutants. Since alanine can be assumed to not significantly contribute to binding due to its small nonpolar side chain, the change in binding free energy upon an alanine mutation reveals interaction “hotspots”, i.e., amino acids that are essential for complex formation. Already in 1999, Massova and Kollman proposed a fast alanine scanning procedure,[40] which has proven to be successful in hot-spot identification[5,41] and is still widely applied today[42–44] to reveal the specific “fingerprints” of interactions. The speed up compared to conventional MMPBSA calculations, in which the binding free energies are calculated based on separate trajectories for wildtype and mutants, is achieved by creating mutants directly from the structures sampled during a simulation of the wild type complex. Despite neglecting conformational adjustments that could occur upon a mutation, this approach has been shown to be as good as[40] or even better[45] than the separate trajectory alternative for all residues that are not involved in highly polar interactions. In contrast, binding free energy contributions of buried residues involved in salt bridges have often been overestimated.[41,45] Thus, the influence on binding of the latter type of residues should be determined by the conventional, separate trajectory approach.[41] Insights into the determinants of binding can be obtained with even less computational effort by applying a decomposition of the binding free energy in terms of substructural elements, such as residues, side chains, or backbone parts, both within the MM-PBSA[46,47] and the MMGBSA[11] frameworks. In the case of protein-protein interactions, these calculations have been shown to yield results in agreement with computational alanine scanning[48] and have even been used in a prospective manner.[49] Thus, MM-PBSA based screens of binding interfaces of biomacro-
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.molinf.com
119
Methods Corner
N. Homeyer, H. Gohlke
molecular complexes have been demonstrated to be very valuable for revealing the specific features of an interaction.
3 Current Limitations Although MM-PBSA proved its strengths in many drug design and biomacromolecular analyses, the studies also revealed weaknesses of the method. Sources of error include the prediction of solute entropies, the estimation of solvation free energies for charged, buried groups, conformational sampling, and parameter selection. Although absolute free energies of binding have been computed that are in good agreement with experimental data,[20,30,37,50] MM-PBSA performs better in determining relative free energies.[22] This is probably due to a cancellation of errors, especially when the relative free energies of similar molecules or conformers are analyzed. The largest uncertainty in MM-PBSA calculations seems to originate from the estimation of the solute entropy.[19] While this contribution can be neglected if relative binding free energies of very similar molecules shall be computed, it needs to be considered for absolute free energies or if significant conformational changes occur upon binding.[22,30] Estimating changes in the vibrational entropy by normal mode analysis probably leads to systematic errors, because anharmonic contributions are neglected.[22,30] Variations in the changes of the vibrational entropy upon complexation of ~ 5 kcal mol 1 have been observed for individually minimized structures from the same trajectory.[31] These significant fluctuations in the entropic contribution are potentially caused by a mismatch in the geometry between minimized structures of the complex and of the receptor or ligand.[22,51] In addition, it has been speculated that the entropy change upon ligand binding is overestimated by normal mode analysis. However, this overestimation seems to be compensated by a likewise overestimation of the effective binding energy of similar magnitude, so that the resulting binding free energies still agreed well with those from experimental studies.[52] In contrast to normal mode analysis, quasiharmonic analysis implicitly takes into account some anharmonic effects.[30] However, with this method, a convergence of the entropy estimates can hardly be reached within currently common simulation times, even for small biomacromolecular structures.[12] Thus, although efforts to improve the accuracy of entropy calculations have been undertaken,[53–55] determining entropic contributions remains a challenging task. Another contribution that may not be accurately captured by MM-PBSA is the solvation free energy of polar compounds and of charged groups, especially those located in the interior of biomacromolecules. As the uncertainty in the estimation of the solvation free energy by the Poisson–Boltzmann method is proportional to the size of this energy, the uncertainty increases with the polarity of the 120
www.molinf.com
studied compounds.[18] Furthermore, while the solvation free energy of solvent-exposed residues of biomolecules can be fairly predicted, computing the free energy of solvation of buried groups is a difficult task due to the inhomogeneity of the interior of the molecules.[52] By using a single dielectric constant for the solute in this case, the inhomogeneous screening of electrostatic interactions cannot be adequately accounted for. This explains to some extent why binding free energies of polar ligands and ion-ion interactions in biomacromolecular complexes are often wrongly estimated.[30,41] In some cases, this problem has been overcome by selecting a dielectric constant for the solute that properly describes the polarization effects within a specific binding pocket.[22,32] Furthermore, contributions of structural water molecules to the binding free energy are not accounted for by implicit solvation models. Explicitly considering such water molecules in the free energy calculations may provide a possibility to overcome this shortcoming of MM-PBSA.[2,56] Representative conformational ensembles are a prerequisite for accurate free energy estimates of a molecule. Therefore, an adequate sampling to cover as much of the conformational space of the solute as possible is essential for ensemble-based MM-PBSA calculations. Among others, the length of the MD simulation required for representative sampling depends on the flexibility of the investigated system. Generally, a stable trajectory, e.g., as judged by the root-mean-square deviation to the initial structure, and stable MM-PBSA energies should be used as minimal criteria for determining the necessary sampling time.[57] For binding free energy estimates of small compounds by the single-trajectory approach short simulation times of 400– 500 ps have been reported to yield converged results for individual systems.[32,58] In contrast, the three-trajectory approach shows usually larger fluctuations in the binding free energy and, hence, requires more extensive sampling.[58] Still, the higher computational costs of the latter approach can be justified in those cases where larger conformational changes take place upon complex formation, as then the conformational strain energy needs to be explicitly taken into account.[37] Whether longer sampling times can improve the accuracy of the free energy estimates also depends on the quality of the parameters used for the description of the molecular system. If the force field cannot correctly describe the system’s features, longer simulations will not lead to better results.[22] Therefore, a good force field is a prerequisite for generating representative conformational ensembles, too. In addition, intrinsic parameters of the method, such as the radii used for the calculation of the solvation free energy and the dielectric constant of the solute, also affect the quality of the results.[12,22,32] Consequently, a careful selection of parameters, force field, and sampling conditions is critical for the success of MM-PBSA free energy estimates.
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Mol. Inf. 2012, 31, 114 – 122
Free Energy Calculations by the Molecular Mechanics Poisson
Boltzmann Surface Area Method
4 Outlook Out of the wide spectrum of applications of MM-PBSA, two areas are especially promising and will probably play a major role in future research: Virtual screening for drugcandidate compounds and residue-wise energy analyses of binding sites to characterize molecular interactions. The usefulness of MM-PBSA for the latter type of application has been proven in numerous studies, and it is highly likely that MM-PBSA will play an important role in the detailed energetic investigation of complex formation in the future. However, as for the field of drug design, the focus of MMPBSA applications will presumably change. As the current accuracy of MM-PBSA hampers the investigation of compounds with small differences in the binding affinities, the main focus will be on re-scoring large compound sets covering a wider range of binding affinities during the lead identification stage of drug design. For lead optimization, the large increase in computer power will soon allow finetuning of drug candidates with the help of more rigorous free energy calculation approaches on a routine basis. In order to be able to compete with the less computationally demanding scoring functions in terms of efficiency, fast high-throughput procedures for MM-PBSA analyses are needed. First efforts in this direction have been undertaken.[27,59] The developed procedures allow screening thousands of compounds within one day, and a special adaptation of MM-PBSA for enterprise grid-computing allows using otherwise idle computers for these calculations.[59,60]
Acknowledgements We are grateful to Bayer Schering Pharma AG for financial support.
References [1] J. Michel, N. Foloppe, J. W. Essex, Mol. Inf. 2010, 29, 570 – 578. [2] T. Hou, J. Wang, Y. Li, W. Wang, J. Comput. Chem. 2011, 32, 866 – 877. [3] J. Srinivasan, T. E. Cheatham III, P. Cieplak, P. A. Kolman, D. A. Case, J. Am. Chem. Soc. 1998, 120, 9401 – 9409. [4] H. Gouda, I. D. Kuntz, D. A. Case, P. A. Kollman, Biopolymers 2003, 68, 16 – 34. [5] N. Homeyer, T. Essigke, H. Meiselbach, G. M. Ullmann, H. Sticht, J. Mol. Model. 2007, 13, 431 – 444. [6] A. Moll, A. Hildebrandt, H. P. Lenhof, O. Kohlbacher, J. Comput.-Aided Mol. Des. 2005, 19, 791 – 800. [7] M. K. Gilson, B. Honig, Proteins: Struct., Funct., Bioinf. 1988, 4, 7 – 18. [8] J. Wang, T. Hou, X. Xu, Curr. Comput.-Aided Drug Des. 2006, 2, 95 – 103. [9] B. Honig, A. Nicholls, Science 1995, 268, 1144 – 1149. [10] M. K. Gilson, K. A. Sharp, B. H. Honig, J. Comput. Chem. 1987, 9, 327 – 335. [11] H. Gohlke, C. Kiel, D. A. Case, J. Mol. Biol. 2003, 330, 891 – 913. [12] H. Gohlke, D. A. Case, J. Comput. Chem. 2004, 25, 238 – 250.
Mol. Inf. 2012, 31, 114 – 122
[13] C. Tan, Y. H. Tan, R. Luo, J. Phys. Chem. B 2007, 111, 12263 – 12274. [14] M. Karplus, J. N. Kushick, Macromolecules 1981, 14, 325 – 332. [15] N. Huang, C. Kalyanaraman, K. Bernacki, M. P. Jacobson, Phys. Chem. Chem. Phys. 2006, 8, 5166 – 5177. [16] A. M. Ferrari, G. Degliesposti, M. Sgobba, G. Rastelli, Bioorg. Med. Chem. 2007, 15, 7865 – 7877. [17] A. Lindstrçm, L. Edvinsson, A. Johansson, C. D. Andersson, I. E. Andersson, F. Raubacher, A. Linusson, J. Chem. Inf. Model. 2011, 51, 267 – 282. [18] J. Kongsted, P. Sçderhjelm, U. Ryde, J. Comput.-Aided Mol. Des. 2009, 23, 395 – 409. [19] A. Weis, K. Katebzadeh, P. Sçderhjelm, I. Nilsson, U. Ryde, J. Med. Chem. 2006, 49, 6596 – 6606. [20] S. Huo, J. Wang, P. Cieplak, P. A. Kollman, I. D. Kuntz, J. Med. Chem. 2002, 45, 1412 – 1419. [21] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, L. G. Pedersen, J. Chem. Phys. 1995, 103, 8577 – 8593. [22] T. Hou, J. Wang, Y. Li, W. Wang, J. Chem. Inf. Model. 2011, 51, 69 – 82. [23] T. Steinbrecher, D. A. Case, A. Labahn, J. Med. Chem. 2006, 49, 1837 – 1844. [24] D. C. Thompson, C. Humblet, D. Joseph-McCarthy, J. Chem. Inf. Model. 2008, 48, 1081 – 1091. [25] G. M. Morris, D. S. Goodsell, R. S. Halliday, R. Huey, W. E. Hart, R. K. Belew, A. J. Olson, J. Comput. Chem. 1998, 19, 1639 – 1662. [26] B. Kuhn, P. Gerber, T. Schulz-Gasch, M. Stahl, J. Med. Chem. 2005, 48, 4040 – 4048. [27] G. Rastelli, A. Del Rio, G. Degliesposti, M. Sgobba, J. Comput. Chem. 2010, 31, 797 – 810. [28] N. Okimoto, N. Futatsugi, H. Fuji, A. Suenaga, G. Morimoto, R. Yanai, Y. Ohno, T. Narumi, M. Taiji, PLoS Comput. Biol. 2009, 5, e1000528. [29] M. D. Parentia, G. Raste, Biotechnol. Adv. 2011, doi: 10.1016/ j.biotechadv2011.08.003. [30] J. Wang, P. Morin, W. Wang, P. A. Kollman, J. Am. Chem. Soc. 2001, 123, 5221 – 5230. [31] B. Kuhn, P. A. Kollman, J. Med. Chem. 2000, 43, 3786 – 3791. [32] T. Yang, J. C. Wu, C. Yan, Y. Wang, R. Luo, M. B. Gonzales, K. N. Dalby, P. Ren, Proteins: Struct., Funct., Bioinf. 2011, 79, 1940 – 1951. [33] A. R. Brice, B. N. Dominy, J. Comput. Chem. 2011, 32, 1431 – 1440. [34] K. Reblova, Z. Strelcova, P. Kulhanek, I. Besseova, D. H. Mathews, K. V. Nostrand, I. Yildirim, D. H. Turner, J. Sponer, J. Chem. Theory Comput. 2010, 6, 910 – 929. [35] M. R. Lee, J. Tsai, D. Baker, P. A. Kollman, J. Mol. Biol. 2001, 313, 417 – 430. [36] C. Combelles, J. Gracy, A. Heitz, D. J. Craik, L. Chiche, Proteins: Struct., Funct., Bioinf. 2008, 73, 87 – 103. [37] J. Lee, J. S. Kim, C. Seok, J. Phys. Chem. B 2010, 114, 7662 – 7671. [38] C. M. Reyes, P. A. Kollman, J. Mol. Biol. 2000, 295, 1 – 6. [39] Y. M. Huang, C. E. A. Chang, BMC Biophysics 2011, 4, 12. [40] I. Massova, P. A. Kollman, J. Am. Chem. Soc. 1999, 121, 8133 – 8143. [41] S. Huo, I. Massova, P. A. Kollman, J. Comput. Chem. 2002, 23, 15 – 27. [42] D. J. Cole, E. Rajendra, M. Roberts-Thomson, B. Hardwick, G. J. McKenzie, M. C. Payne, A. R. Venkitaraman, C. K. Skylaris, PLoS Comput. Biol. 2011, 7, e1002096. [43] N. J. Bruce, D. Chen, S. G. Dastidar, G. E. Marks, C. H. Schein, R. A. Bryce, Peptides 2010, 31, 2100 – 2108.
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.molinf.com
121
Methods Corner
N. Homeyer, H. Gohlke
[44] I. Iacovache, M. T. Degiacomi, L. Pernot, S. Ho, M. Schiltz, M. Dal Peraro, F. G. van der Goot, PLoS Pathog. 2011, 7, e1002135. [45] R. T. Bradshaw, B. H. Patel, E. W. Tate, R. J. Leatherbarrow, I. R. Gould, Protein Eng., Des. Sel. 2011, 24, 197 – 207. [46] A. Metz, C. Pfleger, H. Kopitz, S. Pfeiffer-Marek, K.-H. Baringhaus, H. Gohlke, J. Chem. Inf. Model 2011, doi: 10.1021/ ci200322s. [47] V. Lafont, M. Schaefer, R. H. Stote, D. Altschuh, A. Dejaegere, Proteins: Struct., Funct., Bioinf. 2007, 67, 418 – 434. [48] V. Zoete, O. Michielin, Proteins Struct., Funct., Bioinf. 2007, 67, 1026 – 1047. [49] C. Wichmann, Y. Becker, L. Chen-Wichmann, V. Vogel, A. Vojtkova, J. Herglotz, S. Moore, J. Koch, J. Lausen, W. Mntele, H. Gohlke, M. Grez, Blood 2010, 116, 603 – 613. [50] H. Gouda, Y. Yanai, A. Sugawara, T. Sunazuka, S. Omura, S. Hirono, Bioorg. Med. Chem. 2008, 16, 3565 – 3579. [51] C. S. Page, P. A. Bates, J. Comput. Chem. 2006, 27, 1990 – 2007. [52] N. Singh, A. Warshel, Proteins: Struct., Funct., Bioinf. 2010, 78, 1705 – 1723.
122
www.molinf.com
[53] R. Baron, P. H. Hnenberger, J. A. McCammon, J. Chem. Theory Comput. 2009, 5, 3150 – 3160. [54] J. Kongsted, U. Ryde, J. Comput.-Aided Mol. Des. 2009, 23, 63 – 71. [55] J. Numata, M. Wan, E. W. Knapp, Genome Inf. Ser. 2007, 18, 192 – 205. [56] A. Nurisso, B. Blanchard, A. Audfray, L. Rydner, S. Oscarson, A. Varrot, A. Imberty, J. Biol. Chem. 2010, 285, 20316 – 20327. [57] N. X. Wang, J. J. Zheng, Protein Sci. 2009, 18, 707 – 715. [58] T. Hou, R. Yu, J. Med. Chem. 2007, 50, 1177 – 1188. [59] S. P. Brown, S. W. Muchmore, J. Chem. Inf. Model. 2006, 46, 999 – 1005. [60] S. P. Brown, S. W. Muchmore, J. Chem. Inf. Model. 2007, 47, 1493 – 1503.
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Received: June 27, 2011 Accepted: November 26, 2011 Published online: January 10, 2012
Mol. Inf. 2012, 31, 114 – 122
Free Energy Calculations by the Molecular Mechanics Poisson Boltzmann Surface Area Method Nadine Homeyer[a] and Holger Gohlke*[a]
114
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Mol. Inf. 2012, 31, 114 – 122
Free Energy Calculations by the Molecular Mechanics Poisson
Boltzmann Surface Area Method
Abstract: Detailed knowledge of how molecules recognize interaction partners and of the conformational preferences of biomacromolecules is pivotal for understanding biochemical processes. Such knowledge also provides the foundation for the design of novel molecules, as undertaken in pharmaceutical research. Computer-based free energy calculations enable a detailed investigation of the energetic factors that are responsible for molecular stability or binding affinity. The Molecular Mechanics Poisson–Boltzmann Surface Area (MM-PBSA) approach is an efficient method for the calculation of free energies of diverse mo-
lecular systems. Here we describe the concepts of this approach and outline the practical proceeding. Furthermore we give an overview of the wide spectrum of problems that have been addressed with this method and of successful analyses carried out, thereby focussing on ambitious and recent studies. Limits of the approach in terms of accuracy and applicability are discussed. Despite these limitations MM-PBSA is a method with great potential that allows comparative free energy analyses for various molecular systems at low computational cost.
Keywords: MM-PBSA · Binding affinity · Implicit solvent · Molecular recognition · Drug design
1 Introduction Free energy calculations have become part of the standard repertoire of computational biochemical research due to their widespread applicability, reaching from the prediction of binding affinities of small drug-candidate compounds to the evaluation of relative stabilities of large biomacromolecular structures. The calculations can be performed by the rigorous and computationally expensive free energy perturbation (FEP) and thermodynamic integration (TI) methods, which have recently been described in this ”Methods Corner”.[1] The high computational costs of TI and FEP are caused, first, by the explicit treatment of solvent and, second, by determining the difference in free energies of two states based on simulations that are carried out at intermediate points along a transition path from one state to another. The computationally intensive evaluation of individual contributions of a large number of solvent molecules can be avoided by regarding solvent effects only implicitly, using continuum solvent methods. The computational cost can be further reduced by considering only the end-point states in the free energy calculations. End-point, implicit solvent approaches thus have the advantage of a relatively low computational cost, while important solvent effects are still regarded, so that they provide a trade-off between efficiency and accuracy. As a consequence, with today’s computational resources, the binding of large compound data sets to biomacromolecules can be analyzed within a reasonable time span, and these approaches are more accurate in binding pose and affinity predictions than simple scoring functions.[2] Probably the most well-known end-point, implicit solvent free energy method is the Molecular Mechanics Poisson– Boltzmann Surface Area (MM-PBSA) approach, for the first time so denoted in 1998.[3] The approach has been used in more than one hundred studies to determine the free energies of molecular systems. Topics addressed in these studies include evaluating docking poses, determining structural stability, and predicting binding affinities and hot-spots. This clearly demonstrates the wide applicability of the approach. As a further advantage, MM-PBSA allows analyzing Mol. Inf. 2012, 31, 114 – 122
individual energy contributions by means of a free energy decomposition, which gives additional insights into the energetics of the investigated system. Furthermore, MM-PBSA can be used to study very different molecular systems ranging from small chemical compounds binding to an RNA aptamer to large protein complexes with thousands of atoms.[4,5] In this review we provide an account of selected topics that have been addressed with the MM-PBSA approach, thereby focusing on more recent studies, and imparting knowledge to the reader about planning and evaluating MM-PBSA analyses.
2 State-of-the-Art 2.1 Methodology and Practical Proceeding
In MM-PBSA the free energy of a molecule is estimated as the sum of its gas-phase energy, the solvation free energy, and a contribution due to the configurational entropy of the solute (see Frontispiece[6]). The gas-phase energy is approximated by the molecular mechanics energy of the molecule, determined from a force field comprising terms for bond, angle, and torsion energies as well as van der Waals and electrostatic interactions. In the calculation of the solvation free energy, two contributions are considered, a polar and a nonpolar one. For the polar contribution, the change in the free energy upon transfer of a charged molecule from gas-phase (modeled as a homogeneous medium with dielectric constant e = 1) to solvent (modeled as a homogeneous medium with e = 78 or 80) is estimated by an implicit solvent model.[7,8] In the initial MM-PBSA implementation and still widely used, the polar contribution is calcu[a] N. Homeyer, H. Gohlke Heinrich-Heine-Universitt, Institut fr Pharmazeutische und Medizinische Chemie Dsseldorf, Germany *e-mail: [email protected] Supporting Information for this article is available on the WWW under http://dx.doi.org/10.1002/minf.201100135.
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.molinf.com
115
Methods Corner
N. Homeyer, H. Gohlke
lated via a finite-difference solution of the Poisson–Boltzmann equation described in detail elsewhere.[9,10] Alternatively, an implicit solvent model based on Generalized Born (GB) theory can be used, which is a computationally more efficient approximation to Poisson theory. This then leads to the so-called MM-GBSA variant.[11,12] The nonpolar contribution to the solvation free energy primarily arises from the free energy required for cavity formation for the solute within the solvent. Most frequently, this contribution is estimated by a term proportional to the solvent accessible surface area of the solute. As a more appropriate alternative, the nonpolar contribution is computed as the sum of a disfavorable energy resulting from cavity formation and a favorable energy resulting from attractive interactions between solute and solvent molecules.[12,13] Finally, the configurational entropy of the solute is usually estimated using a rigid-rotor harmonic oscillator approximation, applying either normal mode analysis[3] or quasi-harmonic analysis.[14] From the so-computed free energies of interacting molecules and by applying a thermodynamic cycle a binding free energy can be calculated as the difference between the free energy of a complex and the sum of the free energies of its components.[15] According to Figure 1, this binding free energy is equivalent to the sum of energy and configurational entropy contributions associated with complex formation in the gas-phase and the difference in solvation free energies between the complex and the unbound molecules. Changes in the configurational entropy are frequently neglected if only the relative binding free energy
Nadine Homeyer is a postdoctoral research fellow in the Institute for Pharmaceutical and Medicinal Chemistry at the Heinrich-Heine-University, Dsseldorf. She received her PhD from the University Erlangen-Nuremberg (Germany) in 2008 for her work on the effects of phosphorylation on protein structure, dynamics, and interaction. Her research interests lie in molecular recognition, lead optimization, and the dynamics of biomacromolecules. Holger Gohlke is a Professor of Pharmaceutical and Medicinal Chemistry at the Heinrich-Heine-University, Dsseldorf. His research aims at understanding and predicting receptor-ligand interactions and the modulation of biological processes by pharmacologically relevant molecules. His group develops and applies methods at the interface of computational pharmaceutical and biophysical chemistry and molecular bioinformatics.
116
www.molinf.com
of similar ligands shall be analyzed; this then results in an effective binding energy. Although MM-PBSA calculations can be performed based on single structures, e.g. minimized structures[2,16,17] or crystal structures,[18] they are typically conducted based on a conformational ensemble generated by molecular dynamics (MD) simulations. In the latter case, each energy component is determined by averaging over the respective energy contributions from all conformers in the ensemble. That way, effects of conformational changes taking place in the molecular system are considered in the free energy calculation, and an estimate of the precision of the calculations is obtained. It is common to run MD simulations with explicit solvent for generating the conformational ensemble, because continuum solvent simulations have been shown to yield less accurate results.[19] Usually, the molecular system is solvated in a box of water, and counter ions are added to enforce neutrality. In binding free energy studies it is also possible to just add a water sphere centered on the ligand and move only those atoms located in the core of this sphere during the simulation,[20] which decreases the computational costs. The simulations are typically conducted with the same force field that is used for the MM-PBSA calculations in order to avoid inconsistencies,[19] and long-range electrostatic interactions are computed by the particle-mesh Ewald procedure.[3,21] After the simulations, structures for the conformational ensemble are extracted from the generated MD trajectories. Water molecules and counter ions are removed, and the free energy is calculated as described above. In practice, depending on the information one wishes to gain from the free energy calculations, three main types of MM-PBSA calculations can be distinguished. If the stability of two conformations of a molecule shall be compared, simulations and free energy calculations are conducted separately for both conformers, and the free energies are directly evaluated (Figure 2A). If a binding free energy shall be computed the difference between the free energies of the complex and its components is calculated either from a single trajectory of the complex (“single-trajectory approach”, Figure 2B) or from separate trajectories of complex, receptor, and ligand (“three-trajectory approach”, Figure 2C). Although the singe-trajectory approach neglects the conformational flexibility of the unbound structures, it is usually applied in all cases where no large structural changes upon binding are expected, because it gives less noisy results due to the cancellation of intramolecular contributions[20,22] and, therefore, allows performing MM-PBSA analyses based on shorter simulations. Usual simulation times range from a few picoseconds to several nanoseconds depending on system flexibility, size of the data set, and desired precision of the results, although longer simulation times are recommended for obtaining properly thermalized conformational ensembles.
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Mol. Inf. 2012, 31, 114 – 122
Free Energy Calculations by the Molecular Mechanics Poisson
Boltzmann Surface Area Method
Figure 1. Thermodynamic cycle for calculating a binding free energy. In MM-PBSA, the binding free energy results as the sum of the energy and configurational entropy associated with complex formation in the gas-phase and the difference in solvation free energies between the complex and the unbound molecules.
2.2 Application in Drug Design
The lower computational costs, compared to TI and FEP, together with a more sophisticated computation of the free energy components, compared to common scoring functions, makes MM-PBSA an attractive method for all stages of drug design. In early stages, large compound data bases are screened for promising drug candidates using docking programs and scoring functions to assess the binding capability of the molecules to a target receptor. However, simple scoring functions may neglect important energetic contributions, such as the solvation free energy, and, thus, may have difficulties in correctly predicting the binding pose and binding affinity of a molecule. A combination of docking, molecular mechanics calculations, and MM-PBSA Mol. Inf. 2012, 31, 114 – 122
proved to be more successful in the identification of the correct ligand binding pose than pure docking approaches in several cases.[2,17,23,24] As such, in an extensive performance study, Huo et al.[2] analyzed the performance of MM-PBSA to discriminate correct from false ligand binding poses in 98 complexes. Initially, ligands were docked to receptors with the Autodock[25] program. For the 100 best-scored docking poses as well as the experimentally determined binding pose of each complex the positions of all atoms within a sphere of 9 centered at the ligand were optimized by molecular mechanics minimization. The optimized structures were reevaluated using implicit solvent free energy calculations. Although outperformed by another implicit solvent ap-
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.molinf.com
117
Methods Corner
N. Homeyer, H. Gohlke
Figure 2. Main types of MM-PBSA calculations: Comparison of the stability of two conformations (A) as well as binding free energy calculations according to the single-trajectory approach (B) and the three-trajectory approach (C).
proach, MM-PBSA reached a success rate of 79.6 % in the re-evaluation of the structures, if the top five scored poses were regarded. While, until recently, the use of MM-PBSA in virtual screening has been limited to evaluating a few hundred docking solutions,[26] it can nowadays be applied for scoring thousands of compounds due to the great increase in computer power. As the potential of MM-PBSA in discriminating true binders from a much larger number of decoys has been demonstrated in initial high-throughput virtual screening studies,[27,28] we expect the method to be used more extensively in lead identification in the future. However, the general benefit of MM-PBSA in this field still requires further validation on more protein systems and ligand data sets.[29] When leads have been identified they need to be further improved in order to obtain drug-candidates with good pharmacokinetic (absorption, distribution, metabolism, excretion, and toxicity) features that, at the same time, show a high binding affinity and selectivity for the target. In this stage of drug design, it is essential to be able to determine the relative binding affinity of a small number of drug-candidates, usually less than 50, with high accuracy. The main challenge for the applied method lies in the small range of binding affinities studied, which is typically below 10 kcal mol 1. In initial studies, promising compound rankings were described based on MM-PBSA relative binding affinities for single receptor systems.[20,30,31] However, it has been recognized recently that a more extensive testing is needed to validate the general applicability of the approach for lead optimization. With the required computational power now at hand, extensive performance studies investigating the 118
www.molinf.com
accuracy of MM-PBSA for determining relative binding affinities for several different protein-ligand systems can be conducted. In one such study, Huo et al.[22] calculated the binding free energies of 59 ligands interacting with six different proteins and found good correlations with r2 > 0.5 for three out of the six protein systems studied. In addition, Yang et al.[32] recently described the results of MM-PBSA analyses in which they determined the binding affinities of a total of 156 ligands to seven different protein families separated into six groups. Here, correlations of r2 > 0.5 were found for three out of six groups. However, it needs to be considered that the correlations might have come out weaker due to the sequence variation within the protein families in this case. Overall the results of the two studies suggest that the goodness of the agreement between experimentally determined binding affinities and those computed with MMPBSA does not only depend on the binding affinity range of the analyzed ligand data set, but is also determined by the specific features of the protein-ligand interaction and the degree of structural similarity between the studied compounds. While for sets of ligands with similar structures usually satisfactory results were obtained, the performance varied for those with diverse chemical structures.[22] Furthermore, the number of polar or charged groups involved in ligand binding seems to have a high influence on the calculated binding free energies. For systems with a large number of charged residues in the binding pocket and, hence, where a strong change in the polarization occurs upon binding, using e > 1 for the “gas-phase state” was shown to yield better results.[22,32] Thus, carefully choosing e according to the expected degree of polarization in the
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Mol. Inf. 2012, 31, 114 – 122
Free Energy Calculations by the Molecular Mechanics Poisson
Boltzmann Surface Area Method
binding pocket and/or based on test calculations with evaluation data may prove beneficial. Finally, applications of MM-PBSA in lead optimization have shown that for ligand data sets with a binding affinity range of less than 3 kcal mol 1 only unsatisfactory correlations of experimentally determined and computed relative binding affinities can be expected (Table S1). Although exceptions to this observation have been described for individual data sets,[22] a statistical uncertainty of at least 1.1 kcal mol 1 mean absolute deviation in the solvation free energy relative to a method-weighted average[18] and a standard deviation of 0.7–1.0 kcal mol 1 in the MM-PBSA binding free energy, found for 200 snapshots from a 2 ns simulation,[19] make it difficult to obtain good relative binding free energy predictions for ligand data sets with small binding affinity ranges. For data sets with wider binding affinity ranges the quality of the prediction depends on the specific features of the binding interaction and on the similarity of the ligands. The best results are obtained for ligand series with high structural similarity and/or high hydrophobicity as well as with a uniform total charge. 2.3 Application to Biomacromolecules
Free energy studies on biomacromolecules can give important insights into the structures and mechanisms underlying biological processes. As such, MM-PBSA has been applied to study the stability of DNA,[3,33] RNA,[34] and protein[35] conformers in order to identify the native fold among a set of decoy structures or the fold that is most stable in a given environment. Recently, it has been demonstrated that MM-PBSA is also able to predict disulfide connectivities in proteins.[36] When small cysteine-rich proteins of the knottin family with different possible patterns of disulfide bridges and similar overall structure were analyzed with MM-PBSA, the connectivities proposed by previous experimental studies were identified as being energetically most favorable. This study showed that, in those cases where determining the correct cysteine pairing is difficult with X-ray crystallography and nuclear magnetic resonance spectroscopy, MM-PBSA free energy calculations can help to identify the most favorable pattern. However, it has to be considered that the energetically most stable pattern of disulfide bridges does not need to be the connectivity found in nature, where the oxidative folding pathway might be kinetically controlled. Binding free energy calculations of biomacromolecules can yield information about the stability of complexes and the binding mechanisms. For example, additional knowledge about DNA binding of Cys2His2 zinc finger proteins could be gained from an energetic analysis of the interaction between the Zif268 zinc finger protein and DNA.[37] Binding affinity calculations for one-finger, two-finger, and three-finger Zif268 protein/DNA complexes revealed that the binding affinity does not increase linearly with the number of zinc fingers; rather, it is mainly enhanced by the Mol. Inf. 2012, 31, 114 – 122
addition of the third zinc finger domain. Molecular dynamics simulations and detailed energy decomposition analyses showed that the non-linear cooperativity in Zif268/DNA binding is primarily caused by a tighter interaction between these components, enabled by an enlargement of the DNA major grove and an unwinding of the DNA helix upon attachment of the third zinc finger. The gain in polar interaction energy caused by the tight DNA contact of the three finger protein is much larger than the cost of the DNA conformational change, so that the binding affinity is highly favorable for the protein with three zinc finger domains. The specificity of protein-protein, protein-DNA, and protein-RNA interactions can be studied by modifying residues involved in complex formation and analyzing the differences in binding free energies caused by these modifications. Such modifications can involve interconversions of DNA or RNA bases,[37,38] transformations of phospho-amino acid residues,[39] and amino acid mutations.[38] As for the latter, it is a popular approach to investigate protein-protein interactions by mutating residues in the contact interface consecutively to alanine and, then, calculating the binding free energy for wildtype and mutants. Since alanine can be assumed to not significantly contribute to binding due to its small nonpolar side chain, the change in binding free energy upon an alanine mutation reveals interaction “hotspots”, i.e., amino acids that are essential for complex formation. Already in 1999, Massova and Kollman proposed a fast alanine scanning procedure,[40] which has proven to be successful in hot-spot identification[5,41] and is still widely applied today[42–44] to reveal the specific “fingerprints” of interactions. The speed up compared to conventional MMPBSA calculations, in which the binding free energies are calculated based on separate trajectories for wildtype and mutants, is achieved by creating mutants directly from the structures sampled during a simulation of the wild type complex. Despite neglecting conformational adjustments that could occur upon a mutation, this approach has been shown to be as good as[40] or even better[45] than the separate trajectory alternative for all residues that are not involved in highly polar interactions. In contrast, binding free energy contributions of buried residues involved in salt bridges have often been overestimated.[41,45] Thus, the influence on binding of the latter type of residues should be determined by the conventional, separate trajectory approach.[41] Insights into the determinants of binding can be obtained with even less computational effort by applying a decomposition of the binding free energy in terms of substructural elements, such as residues, side chains, or backbone parts, both within the MM-PBSA[46,47] and the MMGBSA[11] frameworks. In the case of protein-protein interactions, these calculations have been shown to yield results in agreement with computational alanine scanning[48] and have even been used in a prospective manner.[49] Thus, MM-PBSA based screens of binding interfaces of biomacro-
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.molinf.com
119
Methods Corner
N. Homeyer, H. Gohlke
molecular complexes have been demonstrated to be very valuable for revealing the specific features of an interaction.
3 Current Limitations Although MM-PBSA proved its strengths in many drug design and biomacromolecular analyses, the studies also revealed weaknesses of the method. Sources of error include the prediction of solute entropies, the estimation of solvation free energies for charged, buried groups, conformational sampling, and parameter selection. Although absolute free energies of binding have been computed that are in good agreement with experimental data,[20,30,37,50] MM-PBSA performs better in determining relative free energies.[22] This is probably due to a cancellation of errors, especially when the relative free energies of similar molecules or conformers are analyzed. The largest uncertainty in MM-PBSA calculations seems to originate from the estimation of the solute entropy.[19] While this contribution can be neglected if relative binding free energies of very similar molecules shall be computed, it needs to be considered for absolute free energies or if significant conformational changes occur upon binding.[22,30] Estimating changes in the vibrational entropy by normal mode analysis probably leads to systematic errors, because anharmonic contributions are neglected.[22,30] Variations in the changes of the vibrational entropy upon complexation of ~ 5 kcal mol 1 have been observed for individually minimized structures from the same trajectory.[31] These significant fluctuations in the entropic contribution are potentially caused by a mismatch in the geometry between minimized structures of the complex and of the receptor or ligand.[22,51] In addition, it has been speculated that the entropy change upon ligand binding is overestimated by normal mode analysis. However, this overestimation seems to be compensated by a likewise overestimation of the effective binding energy of similar magnitude, so that the resulting binding free energies still agreed well with those from experimental studies.[52] In contrast to normal mode analysis, quasiharmonic analysis implicitly takes into account some anharmonic effects.[30] However, with this method, a convergence of the entropy estimates can hardly be reached within currently common simulation times, even for small biomacromolecular structures.[12] Thus, although efforts to improve the accuracy of entropy calculations have been undertaken,[53–55] determining entropic contributions remains a challenging task. Another contribution that may not be accurately captured by MM-PBSA is the solvation free energy of polar compounds and of charged groups, especially those located in the interior of biomacromolecules. As the uncertainty in the estimation of the solvation free energy by the Poisson–Boltzmann method is proportional to the size of this energy, the uncertainty increases with the polarity of the 120
www.molinf.com
studied compounds.[18] Furthermore, while the solvation free energy of solvent-exposed residues of biomolecules can be fairly predicted, computing the free energy of solvation of buried groups is a difficult task due to the inhomogeneity of the interior of the molecules.[52] By using a single dielectric constant for the solute in this case, the inhomogeneous screening of electrostatic interactions cannot be adequately accounted for. This explains to some extent why binding free energies of polar ligands and ion-ion interactions in biomacromolecular complexes are often wrongly estimated.[30,41] In some cases, this problem has been overcome by selecting a dielectric constant for the solute that properly describes the polarization effects within a specific binding pocket.[22,32] Furthermore, contributions of structural water molecules to the binding free energy are not accounted for by implicit solvation models. Explicitly considering such water molecules in the free energy calculations may provide a possibility to overcome this shortcoming of MM-PBSA.[2,56] Representative conformational ensembles are a prerequisite for accurate free energy estimates of a molecule. Therefore, an adequate sampling to cover as much of the conformational space of the solute as possible is essential for ensemble-based MM-PBSA calculations. Among others, the length of the MD simulation required for representative sampling depends on the flexibility of the investigated system. Generally, a stable trajectory, e.g., as judged by the root-mean-square deviation to the initial structure, and stable MM-PBSA energies should be used as minimal criteria for determining the necessary sampling time.[57] For binding free energy estimates of small compounds by the single-trajectory approach short simulation times of 400– 500 ps have been reported to yield converged results for individual systems.[32,58] In contrast, the three-trajectory approach shows usually larger fluctuations in the binding free energy and, hence, requires more extensive sampling.[58] Still, the higher computational costs of the latter approach can be justified in those cases where larger conformational changes take place upon complex formation, as then the conformational strain energy needs to be explicitly taken into account.[37] Whether longer sampling times can improve the accuracy of the free energy estimates also depends on the quality of the parameters used for the description of the molecular system. If the force field cannot correctly describe the system’s features, longer simulations will not lead to better results.[22] Therefore, a good force field is a prerequisite for generating representative conformational ensembles, too. In addition, intrinsic parameters of the method, such as the radii used for the calculation of the solvation free energy and the dielectric constant of the solute, also affect the quality of the results.[12,22,32] Consequently, a careful selection of parameters, force field, and sampling conditions is critical for the success of MM-PBSA free energy estimates.
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Mol. Inf. 2012, 31, 114 – 122
Free Energy Calculations by the Molecular Mechanics Poisson
Boltzmann Surface Area Method
4 Outlook Out of the wide spectrum of applications of MM-PBSA, two areas are especially promising and will probably play a major role in future research: Virtual screening for drugcandidate compounds and residue-wise energy analyses of binding sites to characterize molecular interactions. The usefulness of MM-PBSA for the latter type of application has been proven in numerous studies, and it is highly likely that MM-PBSA will play an important role in the detailed energetic investigation of complex formation in the future. However, as for the field of drug design, the focus of MMPBSA applications will presumably change. As the current accuracy of MM-PBSA hampers the investigation of compounds with small differences in the binding affinities, the main focus will be on re-scoring large compound sets covering a wider range of binding affinities during the lead identification stage of drug design. For lead optimization, the large increase in computer power will soon allow finetuning of drug candidates with the help of more rigorous free energy calculation approaches on a routine basis. In order to be able to compete with the less computationally demanding scoring functions in terms of efficiency, fast high-throughput procedures for MM-PBSA analyses are needed. First efforts in this direction have been undertaken.[27,59] The developed procedures allow screening thousands of compounds within one day, and a special adaptation of MM-PBSA for enterprise grid-computing allows using otherwise idle computers for these calculations.[59,60]
Acknowledgements We are grateful to Bayer Schering Pharma AG for financial support.
References [1] J. Michel, N. Foloppe, J. W. Essex, Mol. Inf. 2010, 29, 570 – 578. [2] T. Hou, J. Wang, Y. Li, W. Wang, J. Comput. Chem. 2011, 32, 866 – 877. [3] J. Srinivasan, T. E. Cheatham III, P. Cieplak, P. A. Kolman, D. A. Case, J. Am. Chem. Soc. 1998, 120, 9401 – 9409. [4] H. Gouda, I. D. Kuntz, D. A. Case, P. A. Kollman, Biopolymers 2003, 68, 16 – 34. [5] N. Homeyer, T. Essigke, H. Meiselbach, G. M. Ullmann, H. Sticht, J. Mol. Model. 2007, 13, 431 – 444. [6] A. Moll, A. Hildebrandt, H. P. Lenhof, O. Kohlbacher, J. Comput.-Aided Mol. Des. 2005, 19, 791 – 800. [7] M. K. Gilson, B. Honig, Proteins: Struct., Funct., Bioinf. 1988, 4, 7 – 18. [8] J. Wang, T. Hou, X. Xu, Curr. Comput.-Aided Drug Des. 2006, 2, 95 – 103. [9] B. Honig, A. Nicholls, Science 1995, 268, 1144 – 1149. [10] M. K. Gilson, K. A. Sharp, B. H. Honig, J. Comput. Chem. 1987, 9, 327 – 335. [11] H. Gohlke, C. Kiel, D. A. Case, J. Mol. Biol. 2003, 330, 891 – 913. [12] H. Gohlke, D. A. Case, J. Comput. Chem. 2004, 25, 238 – 250.
Mol. Inf. 2012, 31, 114 – 122
[13] C. Tan, Y. H. Tan, R. Luo, J. Phys. Chem. B 2007, 111, 12263 – 12274. [14] M. Karplus, J. N. Kushick, Macromolecules 1981, 14, 325 – 332. [15] N. Huang, C. Kalyanaraman, K. Bernacki, M. P. Jacobson, Phys. Chem. Chem. Phys. 2006, 8, 5166 – 5177. [16] A. M. Ferrari, G. Degliesposti, M. Sgobba, G. Rastelli, Bioorg. Med. Chem. 2007, 15, 7865 – 7877. [17] A. Lindstrçm, L. Edvinsson, A. Johansson, C. D. Andersson, I. E. Andersson, F. Raubacher, A. Linusson, J. Chem. Inf. Model. 2011, 51, 267 – 282. [18] J. Kongsted, P. Sçderhjelm, U. Ryde, J. Comput.-Aided Mol. Des. 2009, 23, 395 – 409. [19] A. Weis, K. Katebzadeh, P. Sçderhjelm, I. Nilsson, U. Ryde, J. Med. Chem. 2006, 49, 6596 – 6606. [20] S. Huo, J. Wang, P. Cieplak, P. A. Kollman, I. D. Kuntz, J. Med. Chem. 2002, 45, 1412 – 1419. [21] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, L. G. Pedersen, J. Chem. Phys. 1995, 103, 8577 – 8593. [22] T. Hou, J. Wang, Y. Li, W. Wang, J. Chem. Inf. Model. 2011, 51, 69 – 82. [23] T. Steinbrecher, D. A. Case, A. Labahn, J. Med. Chem. 2006, 49, 1837 – 1844. [24] D. C. Thompson, C. Humblet, D. Joseph-McCarthy, J. Chem. Inf. Model. 2008, 48, 1081 – 1091. [25] G. M. Morris, D. S. Goodsell, R. S. Halliday, R. Huey, W. E. Hart, R. K. Belew, A. J. Olson, J. Comput. Chem. 1998, 19, 1639 – 1662. [26] B. Kuhn, P. Gerber, T. Schulz-Gasch, M. Stahl, J. Med. Chem. 2005, 48, 4040 – 4048. [27] G. Rastelli, A. Del Rio, G. Degliesposti, M. Sgobba, J. Comput. Chem. 2010, 31, 797 – 810. [28] N. Okimoto, N. Futatsugi, H. Fuji, A. Suenaga, G. Morimoto, R. Yanai, Y. Ohno, T. Narumi, M. Taiji, PLoS Comput. Biol. 2009, 5, e1000528. [29] M. D. Parentia, G. Raste, Biotechnol. Adv. 2011, doi: 10.1016/ j.biotechadv2011.08.003. [30] J. Wang, P. Morin, W. Wang, P. A. Kollman, J. Am. Chem. Soc. 2001, 123, 5221 – 5230. [31] B. Kuhn, P. A. Kollman, J. Med. Chem. 2000, 43, 3786 – 3791. [32] T. Yang, J. C. Wu, C. Yan, Y. Wang, R. Luo, M. B. Gonzales, K. N. Dalby, P. Ren, Proteins: Struct., Funct., Bioinf. 2011, 79, 1940 – 1951. [33] A. R. Brice, B. N. Dominy, J. Comput. Chem. 2011, 32, 1431 – 1440. [34] K. Reblova, Z. Strelcova, P. Kulhanek, I. Besseova, D. H. Mathews, K. V. Nostrand, I. Yildirim, D. H. Turner, J. Sponer, J. Chem. Theory Comput. 2010, 6, 910 – 929. [35] M. R. Lee, J. Tsai, D. Baker, P. A. Kollman, J. Mol. Biol. 2001, 313, 417 – 430. [36] C. Combelles, J. Gracy, A. Heitz, D. J. Craik, L. Chiche, Proteins: Struct., Funct., Bioinf. 2008, 73, 87 – 103. [37] J. Lee, J. S. Kim, C. Seok, J. Phys. Chem. B 2010, 114, 7662 – 7671. [38] C. M. Reyes, P. A. Kollman, J. Mol. Biol. 2000, 295, 1 – 6. [39] Y. M. Huang, C. E. A. Chang, BMC Biophysics 2011, 4, 12. [40] I. Massova, P. A. Kollman, J. Am. Chem. Soc. 1999, 121, 8133 – 8143. [41] S. Huo, I. Massova, P. A. Kollman, J. Comput. Chem. 2002, 23, 15 – 27. [42] D. J. Cole, E. Rajendra, M. Roberts-Thomson, B. Hardwick, G. J. McKenzie, M. C. Payne, A. R. Venkitaraman, C. K. Skylaris, PLoS Comput. Biol. 2011, 7, e1002096. [43] N. J. Bruce, D. Chen, S. G. Dastidar, G. E. Marks, C. H. Schein, R. A. Bryce, Peptides 2010, 31, 2100 – 2108.
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
www.molinf.com
121
Methods Corner
N. Homeyer, H. Gohlke
[44] I. Iacovache, M. T. Degiacomi, L. Pernot, S. Ho, M. Schiltz, M. Dal Peraro, F. G. van der Goot, PLoS Pathog. 2011, 7, e1002135. [45] R. T. Bradshaw, B. H. Patel, E. W. Tate, R. J. Leatherbarrow, I. R. Gould, Protein Eng., Des. Sel. 2011, 24, 197 – 207. [46] A. Metz, C. Pfleger, H. Kopitz, S. Pfeiffer-Marek, K.-H. Baringhaus, H. Gohlke, J. Chem. Inf. Model 2011, doi: 10.1021/ ci200322s. [47] V. Lafont, M. Schaefer, R. H. Stote, D. Altschuh, A. Dejaegere, Proteins: Struct., Funct., Bioinf. 2007, 67, 418 – 434. [48] V. Zoete, O. Michielin, Proteins Struct., Funct., Bioinf. 2007, 67, 1026 – 1047. [49] C. Wichmann, Y. Becker, L. Chen-Wichmann, V. Vogel, A. Vojtkova, J. Herglotz, S. Moore, J. Koch, J. Lausen, W. Mntele, H. Gohlke, M. Grez, Blood 2010, 116, 603 – 613. [50] H. Gouda, Y. Yanai, A. Sugawara, T. Sunazuka, S. Omura, S. Hirono, Bioorg. Med. Chem. 2008, 16, 3565 – 3579. [51] C. S. Page, P. A. Bates, J. Comput. Chem. 2006, 27, 1990 – 2007. [52] N. Singh, A. Warshel, Proteins: Struct., Funct., Bioinf. 2010, 78, 1705 – 1723.
122
www.molinf.com
[53] R. Baron, P. H. Hnenberger, J. A. McCammon, J. Chem. Theory Comput. 2009, 5, 3150 – 3160. [54] J. Kongsted, U. Ryde, J. Comput.-Aided Mol. Des. 2009, 23, 63 – 71. [55] J. Numata, M. Wan, E. W. Knapp, Genome Inf. Ser. 2007, 18, 192 – 205. [56] A. Nurisso, B. Blanchard, A. Audfray, L. Rydner, S. Oscarson, A. Varrot, A. Imberty, J. Biol. Chem. 2010, 285, 20316 – 20327. [57] N. X. Wang, J. J. Zheng, Protein Sci. 2009, 18, 707 – 715. [58] T. Hou, R. Yu, J. Med. Chem. 2007, 50, 1177 – 1188. [59] S. P. Brown, S. W. Muchmore, J. Chem. Inf. Model. 2006, 46, 999 – 1005. [60] S. P. Brown, S. W. Muchmore, J. Chem. Inf. Model. 2007, 47, 1493 – 1503.
2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim
Received: June 27, 2011 Accepted: November 26, 2011 Published online: January 10, 2012
Mol. Inf. 2012, 31, 114 – 122
