Ulf Ryde

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.

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. 15N and 2H 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.

On the Convergence of QM/MM Energies

We have studied the convergence of QM/MM calculations with respect to the size of the QM system. We study a proton transfer between a first-sphere cysteine ligand and a second-sphere histidine group in [Ni,Fe] hydrogenase and use a 446-atom model of the protein, treated purely with QM methods as a reference. We have tested 12 different ways to redistribute charges close to the junctions (to avoid overpolarization of the QM system), but once the junctions are moved away from the active site, there is little need to redistribute the charges. We have tested 13 different variants of QM/MM approaches, including two schemes to correct errors caused by the truncation of the QM system. However, we see little gain from such correction schemes; on the contrary, they are sensitive to the charge-redistribution scheme and may cause large errors if charges are close to the junctions. In fact, the best results were obtained with a mechanical embedding approach that does not employ any correction scheme and ignores polarization. It gives a mean unsigned error for 40 QM systems of different sizes of 7 kJ/mol with a maximum error of 28 kJ/mol. The errors can be significantly decreased if bonds between the QM and MM system (junctions) are moved one residue away from all active-site residues. Then, most QM/MM variants give mean unsigned errors of 5-9 kJ/mol, maximum errors of 16-35 kJ/mol, and only five to seven residues give an error of over 5 kJ/mol. In general, QM/MM calculations converge faster with system size than pure QM calculations.

How to obtain statistically converged MM/GBSA results.

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

The coordination of the catalytic zinc ion in alcohol dehydrogenase studied by combined quantum-chemical and molecular mechanics calculations

SummaryThe 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.

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.

Accurate QM/MM free energy calculations of enzyme reactions: Methylation by catechol O-methyltransferase

We recently described a method to compute accurate quantum mechanical free energies [Rod, T. H.; Ryde, U. Phys. Rev. Lett. 2005, 94, 138302]. The method, which we term quantum mechanical thermodynamic cycle perturbation (QTCP), employs a molecular mechanics force field to sample phase space and, subsequently, a thermodynamic cycle to estimate QM/MM free energy changes. Here, we discuss the methodology in detail and test an approach based on a different thermodynamic cycle. We also show that a new way of treating hydrogen link atoms makes the free energy changes converge faster and that extrapolation to higher accuracy can be performed. We finally discuss the quantum mechanical free energy (QM/MM-FE) method in the framework of the QTCP method. All methods considered are applied to the methylation of catecholate catalyzed by catechol O-methyltransferase. We compute the free energy barrier for the reaction by computing free energy changes in steps between fixed QM regions along a predetermined reaction pathway. Using the QTCP approach, an extrapolated activation free energy of 69 kJ/mol for the forward reaction and 90 kJ/mol for the reverse reaction are obtained at the level of the B3LYP functional and the 6-311++G(2d,2p) basis set. The value for the forward reaction is in excellent agreement with the experimental value of 75 kJ/mol. Results based on the QM/MM-FE method differ by less than 10 kJ/mol from those values, indicating that QM/MM-FE may be a fairly accurate and cheap alternative to calculate QM/MM free energy changes. Moreover, the results are compared to barriers obtained with a fixed molecular mechanics environment as well as with structures optimized in a vacuum. All the computed free energy barriers are well converged. A major approximation in the current implementation of the QTCP method is that the QM region is fixed. The approximation leads to well-converged free energy barriers, which has been a problem in similar studies.

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.

How Accurate Can a Force Field Become? A Polarizable Multipole Model Combined with Fragment-wise Quantum-Mechanical Calculations

A new method to accurately estimate the interaction energy between a large molecule and a smaller ligand is presented. The method approximates the electrostatic and induction contributions classically by multipole and polarizability expansions, but uses explicit quantum-mechanical fragment calculations for the remaining (nonclassical) contributions, mainly dispersion and exchange repulsion. Thus, it represents a limit of how accurate a force field can ever become for interaction energies if pairwise additivity of the nonclassical term is assumed (e.g., all general-purpose force fields). The accuracy is tested by considering protein-ligand model systems for which the true MP2/6-31G* interaction energies can be computed. The method is shown to be more accurate than related fragmentation approaches. The remaining error (2-5 and approximately10 kJ/mol for neutral and charged ligands, respectively) can be decreased by including the polarizing effect from surrounding fragments in the quantum-mechanical calculations.