The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities

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.

Ligand Affinities Predicted with the MM/PBSA Method: Dependence on the Simulation Method and the Force Field

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.

Protein flexibility and conformational entropy in ligand design targeting the carbohydrate recognition domain of galectin-3

Rational drug design is predicated on knowledge of the three-dimensional structure of the protein-ligand complex and the thermodynamics of ligand binding. Despite the fundamental importance of both enthalpy and entropy in driving ligand binding, the role of conformational entropy is rarely addressed in drug design. In this work, we have probed the conformational entropy and its relative contribution to the free energy of ligand binding to the carbohydrate recognition domain of galectin-3. Using a combination of NMR spectroscopy, isothermal titration calorimetry, and X-ray crystallography, we characterized the binding of three ligands with dissociation constants ranging over 2 orders of magnitude. (15)N and (2)H spin relaxation measurements showed that the protein backbone and side chains respond to ligand binding by increased conformational fluctuations, on average, that differ among the three ligand-bound states. Variability in the response to ligand binding is prominent in the hydrophobic core, where a distal cluster of methyl groups becomes more rigid, whereas methyl groups closer to the binding site become more flexible. The results reveal an intricate interplay between structure and conformational fluctuations in the different complexes that fine-tunes the affinity. The estimated change in conformational entropy is comparable in magnitude to the binding enthalpy, demonstrating that it contributes favorably and significantly to ligand binding. We speculate that the relatively weak inherent protein-carbohydrate interactions and limited hydrophobic effect associated with oligosaccharide binding might have exerted evolutionary pressure on carbohydrate-binding proteins to increase the affinity by means of conformational entropy.

Performance of density functionals for first row transition metal systems

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 coordination of the catalytic zinc in alcohol dehydrogenase studied by combined quantum-chemical and molecular mechanics calculations

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.

The cupric geometry of blue copper proteins is not strained

The geometry of several realistic models of the metal coordination sphere in the blue copper proteins has been optimised using high-level quantum chemical methods. The results show that the optimal vacuum structure of the Cu(II) models is virtually identical to the crystal structure of oxidised blue copper proteins. For the reduced forms, the optimised structure seems to be more tetrahedral than the one found in the proteins, but the energy difference between the two geometries is less than 5 kJ/mol, i.e. within the error limits of the method. Thus, the results raise strong doubts against hypotheses (entatic state and the induced-rack theory) suggesting that blue copper proteins force the oxidised metal coordination sphere into a structure similar to that preferred by Cu(I) in order to minimise the reorganisation energy of the electron transfer reaction. Instead, a small reorganisation energy seems to be reached by an appropriate choice of metal ligands. In particular, the cysteine thiolate ligand appears to be crucial, changing the preferred geometry of the oxidised complexes from square-planar to a more trigonal geometry.

On the role of the axial ligand in heme proteins: a theoretical study

We present a systematic investigation of how the axial ligand in heme proteins influences the geometry, electronic structure, and spin states of the active site, and the energies of the reaction cycles. Using the density functional B3LYP method and medium-sized basis sets, we have compared models with His, His+Asp, Cys, Tyr, and Tyr+Arg as found in myoglobin and hemoglobin, peroxidases, cytochrome P450, and heme catalases, respectively. We have studied 12 reactants and intermediates of the reaction cycles of these enzymes, including complexes with H(2)O, OH(-), O(2-), CH(3)OH, O(2), H(2)O(2), and HO(2)(-) in various formal oxidation states of the iron ion (II to V). The results show that His gives ~0.6 V higher reduction potentials than the other ligands. In particular, it is harder to reduce and protonate the O(2) complex with His than with the other ligands, in accordance with the O(2) carrier function of globins and the oxidative chemistry of the other proteins. For most properties, the trend Cys<Tyr<Tyr+Arg<His+Asp<His is found, reflecting the donor capacity of the various ligands. Thus, it is easier to reduce compound I with a His+Asp ligand than with a Cys ligand, in accordance with the one-electron chemistry of peroxidases and the hydroxylation reactions of cytochromes P450. However, the Tyr complexes have an unusually low affinity for all neutral ligands, giving them a slightly enhanced driving force in the oxidation of H(2)O(2) by compound I.

How O2 binds to heme: reasons for rapid binding and spin inversion

We have used density functional methods to calculate fully relaxed potential energy curves of the seven lowest electronic states during the binding of O(2) to a realistic model of ferrous deoxyheme. Beyond a Fe-O distance of approximately 2.5 A, we find a broad crossing region with five electronic states within 15 kJ/mol. The almost parallel surfaces strongly facilitate spin inversion, which is necessary in the reaction of O(2) with heme (deoxyheme is a quintet and O(2) a triplet, whereas oxyheme is a singlet). Thus, despite a small spin-orbit coupling in heme, the transition probability approaches unity. Using reasonable parameters, we estimate a transition probability of 0.06-1, which is at least 15 times larger than for the nonbiological Fe-O(+) system. Spin crossing is anticipated between the singlet ground state of bound oxyheme, the triplet and septet dissociation states, and a quintet intermediate state. The fact that the quintet state is close in energy to the dissociation couple is of biological importance, because it explains how both spin states of O(2) may bind to heme, thereby increasing the overall efficiency of oxygen binding. The activation barrier is estimated to be <15 kJ/mol based on our results and Mössbauer experiments. Our results indicate that both the activation energy and the spin-transition probability are tuned by the porphyrin as well as by the choice of the proximal heme ligand, which is a histidine in the globins. Together, they may accelerate O(2) binding to iron by approximately 10(11) compared with the Fe-O(+) system. A similar near degeneracy between spin states is observed in a ferric deoxyheme model with the histidine ligand hydrogen bonded to a carboxylate group, i.e. a model of heme peroxidases, which bind H(2)O(2) in this oxidation state.

The carbohydrate-binding site in galectin-3 is preorganized to recognize a sugarlike framework of oxygens: ultra-high-resolution structures and water dynamics

The recognition of carbohydrates by proteins is a fundamental aspect of communication within and between living cells. Understanding the molecular basis of carbohydrate–protein interactions is a prerequisite for the rational design of synthetic ligands. Here we report the high- to ultra-high-resolution crystal structures of the carbohydrate recognition domain of galectin-3 (Gal3C) in the ligand-free state (1.08 Å at 100 K, 1.25 Å at 298 K) and in complex with lactose (0.86 Å) or glycerol (0.9 Å). These structures reveal striking similarities in the positions of water and carbohydrate oxygen atoms in all three states, indicating that the binding site of Gal3C is preorganized to coordinate oxygen atoms in an arrangement that is nearly optimal for the recognition of β-galactosides. Deuterium nuclear magnetic resonance (NMR) relaxation dispersion experiments and molecular dynamics simulations demonstrate that all water molecules in the lactose-binding site exchange with bulk water on a time scale of nanoseconds or shorter. Nevertheless, molecular dynamics simulations identify transient water binding at sites that agree well with those observed by crystallography, indicating that the energy landscape of the binding site is maintained in solution. All heavy atoms of glycerol are positioned like the corresponding atoms of lactose in the Gal3C complexes. However, binding of glycerol to Gal3C is insignificant in solution at room temperature, as monitored by NMR spectroscopy or isothermal titration calorimetry under conditions where lactose binding is readily detected. These observations make a case for protein cryo-crystallography as a valuable screening method in fragment-based drug discovery and further suggest that identification of water sites might inform inhibitor design.

Structures of the high-valent metal-ion haem–oxygen intermediates in peroxidases, oxygenases and catalases

Peroxidases, oxygenases and catalases have similar high-valent metal-ion intermediates in their respective reaction cycles. In this review, haem-based examples will be discussed. The intermediates of the haem-containing enzymes have been extensively studied for many years by different spectroscopic methods like UV-Vis, EPR (electron paramagnetic resonance), resonance Raman, Mössbauer and MCD (magnetic circular dichroism). The first crystal structure of one of these high-valent intermediates was on cytochrome c peroxidase in 1987. Since then, structures have appeared for catalases in 1996, 2002, 2003, putatively for cytochrome P450 in 2000, for myoglobin in 2002, for horseradish peroxidase in 2002 and for cytochrome c peroxidase again in 1994 and 2003. This review will focus on the most recent structural investigations for the different intermediates of these proteins. The structures of these intermediates will also be viewed in light of quantum mechanical (QM) calculations on haem models. In particular quantum refinement, which is a combination of QM calculations and crystallography, will be discussed. Only small structural changes accompany the generation of these intermediates. The crystal structures show that the compound I state, with a so called pi-cation radical on the haem group, has a relatively short iron-oxygen bond (1.67-1.76A) in agreement with a double-bond character, while the compound II state or the compound I state with a radical on an amino acid residue have a relatively long iron-oxygen bond (1.86-1.92A) in agreement with a single-bond character where the oxygen-atom is protonated.

Quantum mechanical free energy barrier for an enzymatic reaction

We discuss problems related to in silico studies of enzymes and show that accurate and converged free energy changes for complex chemical reactions can be computed if a method based on a thermodynamic cycle is employed. The method combines the sampling speed of molecular mechanics with the accuracy of a high-level quantum mechanics method. We use the method to compute the free energy barrier for a methyl transfer reaction catalyzed by the enzyme catechol O-methyltransferase at the level of density functional theory. The surrounding protein and solvent are found to have a profound effect on the reaction, and we show that energies can be extrapolated easily from one basis set and exchange-correlation functional to another. Using this procedure we calculate a barrier of 69 kJ/mol, in excellent agreement with the experimental value of 75 kJ/mol.

An MM/3D-RISM approach for ligand binding affinities

We have modified the popular MM/PBSA or MM/GBSA approaches (molecular mechanics for a biomolecule, combined with a Poisson-Boltzmann or generalized Born electrostatic and surface area nonelectrostatic solvation energy) by employing instead the statistical-mechanical, three-dimensional molecular theory of solvation (also known as 3D reference interaction site model, or 3D-RISM-KH) coupled with molecular mechanics or molecular dynamics ( Blinov , N. ; et al. Biophys. J. 2010 ; Luchko , T. ; et al. J. Chem. Theory Comput. 2010 ). Unlike the PBSA or GBSA semiempirical approaches, the 3D-RISM-KH theory yields a full molecular picture of the solvation structure and thermodynamics from the first principles, with proper account of chemical specificities of both solvent and biomolecules, such as hydrogen bonding, hydrophobic interactions, salt bridges, etc. We test the method on the binding of seven biotin analogues to avidin in aqueous solution and show it to work well in predicting the ligand-binding affinities. We have compared the results of 3D-RISM-KH with four different generalized Born and two Poisson-Boltzmann methods. They give absolute binding energies that differ by up to 208 kJ/mol and mean absolute deviations in the relative affinities of 10-43 kJ/mol.

Quantum Chemistry Can Locally Improve Protein Crystal Structures

We have re-refined the X-ray structure of the heme site in cytochrome c553, supplementing the crystallographic data with quantum chemical geometry optimizations, instead of the molecular mechanics force field used in standard crystallographic refinement. By comparing the resulting structure, obtained using medium-resolution data (170 pm), with an atomic-resolution structure (95 pm) of the same protein, we show that the inclusion of quantum chemical information into the refinement procedure improves the structure significantly. Thus, errors in the Fe-ligand distances are reduced from 3 to 32 pm in the low-resolution structure to 0-5 pm in the re-refined structure, one side-chain atom changes its conformation (a movement by 214 pm toward its position in the high-resolution structure), and the R factors are improved by up to 0.018. Thus, quantum refinement may be a powerful method to obtain an accurate structure for interesting parts of a protein.

Metal ion enhanced binding of AMD3100 to Asp262 in the CXCR4 receptor

The affinity of AMD3100, a symmetrical nonpeptide antagonist composed of two 1,4,8,11-tetraazacyclotetradecane (cyclam) rings connected through a 1,4-dimethylene(phenylene) linker to the CXCR4 chemokine receptor was increased 7, 36, and 50-fold, respectively, by incorporation of the following: Cu(2+), Zn(2+), or Ni(2+) into the cyclam rings of the compound. The rank order of the transition metal ions correlated with the calculated binding energy between free acetate and the metal ions coordinated in a cyclam ring. Construction of AMD3100 substituted with only a single Cu(2+) or Ni(2+) ion demonstrated that the increase in binding affinity of the metal ion substituted bicyclam is achieved through an enhanced interaction of just one of the ring systems. Mutational analysis of potential metal ion binding residues in the main ligand binding crevice of the CXCR4 receptor showed that although binding of the bicyclam is dependent on both Asp(171) and Asp(262), the enhancing effect of the metal ion was selectively eliminated by substitution of Asp(262) located at the extracellular end of TM-VI. It is concluded that the increased binding affinity of the metal ion substituted AMD3100 is obtained through enhanced interaction of one of the cyclam ring systems with the carboxylate group of Asp(262). It is suggested that this occurs through a strong concomitant interaction of one of the oxygen's directly with the metal ion and the other oxygen to one of the nitrogens of the cyclam ring through a hydrogen bond.

A QM/MM investigation of the activation and catalytic mechanism of Fe-only hydrogenases

Fe-only hydrogenases are enzymes that catalyze dihydrogen production or oxidation, due to the presence of an unusual Fe(6)S(6) cluster (the so-called H-cluster) in their active site, which is composed of a Fe(2)S(2) subsite, directly involved in catalysis, and a classical Fe(4)S(4) cubane cluster. Here, we present a hybrid quantum mechanical and molecular mechanical (QM/MM) investigation of the Fe-only hydrogenase from Desulfovibrio desulfuricans, in order to unravel key issues regarding the activation of the enzyme from its completely oxidized inactive state (Hoxinact) and the influence of the protein environment on the structural and catalytic properties of the H-cluster. Our results show that the Fe(2)S(2) subcluster in the Fe(II)Fe(II) redox state - which is experimentally observed for the completely oxidized form of the enzyme - binds a water molecule to one of its metal centers. The computed QM/MM energy values for water binding to the diferrous subsite are in fact over 70 kJ mol(-1); however, the affinity toward water decreases by 1 order of magnitude after a one-electron reduction of H(ox)(inact), thus leading to the release of coordinated water from the H-cluster. The investigation of a catalytic cycle of the Fe-only hydrogenase that implies formation of a terminal hydride ion and a di(thiomethyl)amine (DTMA) molecule acting as an acid/base catalyst indicates that all steps have reasonable reaction energies and that the influence of the protein on the thermodynamic profile of H(2) production catalysis is not negligible. QM/MM results show that the interactions between the Fe(2)S(2) subsite and the protein environment could give place to structural rearrangements of the H-cluster functional for catalysis, provided that the bidentate ligand that bridges the iron atoms in the binuclear subsite is actually a DTMA residue.

Carboxylate binding modes in zinc proteins: a theoretical study

The relative energies of different coordination modes (bidentate, monodentate, syn, and anti) of a carboxylate group bound to a zinc ion have been studied by the density functional method B3LYP with large basis sets on realistic models of the active site of several zinc proteins. In positively charged four-coordinate complexes, the mono- and bidentate coordination modes have almost the same energy (within 10 kJ/mol). However, if there are negatively charged ligands other than the carboxylate group, the monodentate binding mode is favored. In general, the energy difference between monodentate and bidentate coordination is small, 4-24 kJ/mol, and it is determined more by hydrogen-bond interactions with other ligands or second-sphere groups than by the zinc-carboxylate interaction. Similarly, the activation energy for the conversion between the two coordination modes is small, approximately 6 kJ/mol, indicating a very flat Zn-O potential surface. The energy difference between syn and anti binding modes of the monodentate carboxylate group is larger, 70-100 kJ/mol, but this figure again strongly depends on interactions with second-sphere molecules. Our results also indicate that the pK(a) of the zinc-bound water ligand in carboxypeptidase and thermolysin is 8-9.

Molecular dynamics simulations of alcohol dehydrogenase with a four- or five-coordinate catalytic zinc ion

A detailed parameterization is presented of a zinc ion with one histidine and two cysteinate ligands, together with one or two water, hydroxide, aldehyde, alcohol, or alkoxide ligands. The parameterization is tailored for the active site of alcohol dehydrogenase and is obtained entirely from quantum chemical computations. The force-field reproduces excellently the geometry of quantum chemically optimized zinc complexes as well as the crystallographic geometry of the active site of alcohol dehydrogenase and small organic structures. The parameterization is used in molecular dynamics simulations and molecular mechanical energy minimizations of alcohol dehydrogenase with a four- or five-coordinate catalytic zinc ion. The active-site zinc ion seems to prefer four-coordination over five-coordination by at least 36 kJ/mol. The only stable binding site of a fifth ligand at the active-site zinc ion is opposite to the normal substrate site, in a narrow cavity behind the zinc ion. Only molecules of the size of water or smaller may occupy this site. There are large fluctuations in the geometry of the zinc coordination sphere. A four-coordinate water molecule alternates frequently (every 7 ps) between the substrate site and the fifth binding site and even two five-coordinate water molecules may interchange ligation sites without prior dissociation. Ligand exchange at the zinc ion probably proceeds by a dissociative mechanism.(ABSTRACT TRUNCATED AT 250 WORDS)