Introduction: The molecular mechanics energies combined with the Poisson–Boltzmann or generalized Born and surface area continuum solvation (MM/PBSA and MM/GBSA) methods are popular approaches to estimate the free energy of the binding of small ligands to biological macromolecules. They are typically based on molecular dynamics simulations of the receptor–ligand complex and are therefore intermediate in both accuracy and computational effort between empirical scoring and strict alchemical perturbation methods. They have been applied to a large number of systems with varying success. Areas covered: The authors review the use of MM/PBSA and MM/GBSA methods to calculate ligand-binding affinities, with an emphasis on calibration, testing and validation, as well as attempts to improve the methods, rather than on specific applications. Expert opinion: MM/PBSA and MM/GBSA are attractive approaches owing to their modular nature and that they do not require calculations on a training set. They have been used successfully to reproduce and rationalize experimental findings and to improve the results of virtual screening and docking. However, they contain several crude and questionable approximations, for example, the lack of conformational entropy and information about the number and free energy of water molecules in the binding site. Moreover, there are many variants of the method and their performance varies strongly with the tested system. Likewise, most attempts to ameliorate the methods with more accurate approaches, for example, quantum-mechanical calculations, polarizable force fields or improved solvation have deteriorated the results.

No abstract

The free energy of binding between avidin and seven biotin analogues has been calculated with the molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) method. We have studied how the force field and the method to generate geometries affect the calculated binding free energies. Four different force fields were compared, but we saw no significant difference in the results. However, mixing the force fields used for the geometry generation and energy calculations is not recommended. In the molecular dynamics simulations, explicit water molecules must be used, but the size of the simulated system and the boundary conditions are less important. In fact, nonperiodic simulations with a fixed protein outside a relatively small simulated system (18 A) seem to be a proper approach. The mean absolute error was 9-19 kJ/mol, with a standard error of 5-15 kJ/mol, which arises mainly from the entropy term.

The coordination number of the catalytic zinc ion in alcohol dehydrogenase has been studied by integrated ab initio quantum-chemical and molecular mechanics geometry optimisations involving the whole enzyme. A four-coordinate active-site zinc ion is 100-200 kJ/mol more stable than a five-coordinate one, depending on the ligands. The only stable binding site for a fifth ligand at the zinc ion is opposite to the normal substrate site, in a small cavity buried behind the zinc ion. The zinc coordination sphere has to be strongly distorted to accommodate a ligand in this site, and the ligand makes awkward contacts with surrounding atoms. Thus, the results do not support proposals attributing an important role to five-coordinate zinc complexes in the catalytic mechanism of alcohol dehydrogenase. The present approach makes it possible also to quantify the strain induced by the enzyme onto the zinc ion and its ligands; it amounts to 42-87 kJ/mol for four-coordinate active-site zinc ion complexes and 131-172 kJ/mol for five-coordinate ones. The four-coordinate structure with a water molecule bound to the zinc ion is about 20 kJ/mol less strained than the corresponding structure with a hydroxide ion, indicating that the enzyme does not speed up the reaction by forcing the zinc coordination sphere into a structure similar to the reaction intermediates.

This article investigates the performance of five commonly used density functionals, B3LYP, BP86, PBE0, PBE, and BLYP, for studying diatomic molecules consisting of a first row transition metal bonded to H, F, Cl, Br, N, C, O, or S. Results have been compared with experiment wherever possible. Open-shell configurations are found more often in the order PBE0>B3LYP>PBE approximately BP86>BLYP. However, on average, 58 of 63 spins are correctly predicted by any functional, with only small differences. BP86 and PBE are slightly better for obtaining geometries, with errors of only 0.020 A. Hybrid functionals tend to overestimate bond lengths by a few picometers and underestimate bond strengths by favoring open shells. Nonhybrid functionals usually overestimate bond energies. All functionals exhibit similar errors in bond energies, between 42 and 53 kJmol. Late transition metals are found to be better modeled by hybrid functionals, whereas nonhybrid functionals tend to have less of a preference. There are systematic errors in predicting certain properties that could be remedied. BLYP performs the best for ionization potentials studied here, PBE0 the worst. In other cases, errors are similar. Finally, there is a clear tendency for hybrid functionals to give larger dipole moments than nonhybrid functionals. These observations may be helpful in choosing and improving existing functionals for tasks involving transition metals, and for designing new, improved functionals.

The copper coordination geometry in the blue copper proteins plastocyanin, nitrite reductase, cucumber basic protein, and azurin has been studied by combined density functional (B3LYP) and molecular mechanical methods. Compared to quantum chemical vacuum calculations, a significant improvement of the geometry is seen (toward the experimental structures) not only for the dihedral angles of the ligands but also for the bond lengths and angles around the copper ion. The flexible Cu–SMet bond is well reproduced in the oxidized structures, whereas it is too long in some of the reduced complexes (too short in vacuum). The change in the geometry compared to the vacuum state costs 33–66 kJ/mol. If the covalent bonds between the ligands and the protein are broken, this energy decreases by ∼25 kJ/mol, which is an estimate of the covalent strain. This is similar to what is found for other proteins, so the blue copper proteins are not more strained than other metalloproteins. The inner‐sphere self‐exchange reorganization energy of all four proteins are ∼30 kJ/mol. This is 30–50 kJ/mol lower than in vacuum. The decrease is caused by dielectric and electrostatic effects in the protein, especially the hydrogen bond(s) to the cysteine copper ligands and not by covalent strain. © 2001 John Wiley & Sons, Inc. Int J Quant Chem 81: 335–347, 2001

Four methods for deriving partial atomic charges from the quantum chemical electrostatic potential (CHELP, CHELPG, Merz-Kollman, and RESP) have been compared and critically evaluated. It is shown the charges strongly depend on how and where the potential points are selected. Two alternative methods are suggested to avoid the arbitrariness in the point-selection schemes and Van der Waals exclusion radii: CHELP-BOW, which also estimate the charges from the electrostatic potential, but with potential points that are Boltzmann-weighted after their occurrence in actual simulations using the energy function of the program in which the charges will be used, and CHELMO, which estimates the charges directly from the electrostatic multipole moments. Different criteria for the quality of the charges are discussed. The CHELMO method gives the best multipole moments for small and medium-sized polar systems, whereas the CHELP-BOW charges reproduce best the total interaction energy in actual simulations. Among the standard methods, the Merz-Kollman charges give the best moments and potentials, but they show an appreciable dependence on the orientation of the molecule.We have also examined the recent warning that charges derived by a least-squares fit to the electrostatic potential normally are not statistically valid. It is shown that no rank-deficiency problems are encountered for molecules with up to 84 atoms if the least-squares fit is performed using pseudoinverses calculated by singular value decomposition and if constraints are treated by elimination.

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.