Search

Understanding the dynamics of ligands bound to proteins is an important task in medicinal chemistry and drug design. However, the dominant technique for determining protein-ligand structures, X-ray crystallography, does not fully account for dynamics and cannot accurately describe the movements of ligands in protein binding sites. In this article, an alternative method, ensemble refinement, is used on six protein-ligand complexes with the aim of understanding the conformational diversity of ligands in protein crystal structures. The results show that ensemble refinement sometimes indicates that the flexibility of parts of the ligand and some protein side chains is larger than that which can be described by a single conformation and atomic displacement parameters. However, since the electron-density maps are comparable and R free values are slightly increased, the original crystal structure is still a better model from a statistical point of view. On the other hand, it is shown that molecular-dynamics simulations and automatic generation of alternative conformations in crystallographic refinement confirm that the flexibility of these groups is larger than is observed in standard refinement. Moreover, the flexible groups in ensemble refinement coincide with groups that give high atomic displacement parameters or non-unity occupancy if optimized in standard refinement. Therefore, the conformational diversity indicated by ensemble refinement seems to be qualitatively correct, indicating that ensemble refinement can be an important complement to standard crystallographic refinement as a tool to discover which parts of crystal structures may show extensive flexibility and therefore are poorly described by a single conformation. However, the diversity of the ensembles is often exaggerated (probably partly owing to the rather poor force field employed) and the ensembles should not be trusted in detail., (open access.)

Conformational entropies are of great interest when studying the binding of small ligands to proteins or the interaction of proteins. Unfortunately, there are no experimental methods available to measure conformational entropies of all groups in a protein. Instead, they are normally estimated from molecular dynamics (MD) simulations, although such methods show problems with convergence and correlation of motions, and depend on the accuracy of the underlying potential-energy function. Crystallographic atomic displacement parameters (also known as B-factors) are available in all crystal structures and contain information about the atomic fluctuations, which can be converted to entropies. We have studied whether B-factors can be employed to extract conformational entropies for proteins by comparing such entropies to those measured by NMR relaxation experiments or obtained from MD simulations in solution or in the crystal. Unfortunately, our results show that B-factor entropies are unreliable, because they include the movement and rotation of the entire protein, they exclude correlation of the movements and they include contributions other than the fluctuations, e.g. static disorder, as well as errors in the model and the scattering factors. We have tried to reduce the first problem by employing translation-libration-screw refinement, the second by employing a description of the correlated movement from MD simulations, and the third by studying only the change in entropy when a pair of ligands binds to the same protein, thoroughly re-refining the structures in exactly the same way and using the same set of alternative conformations. However, the experimental B-factors seem to be incompatible with fluctuations from MD simulations and the precision is too poor to give any reliable entropies.

We have estimated free energies for the binding of eight carboxylate ligands to two variants of the octa-acid deep-cavity host in the SAMPL6 blind-test challenge (with or without endo methyl groups on the four upper-rim benzoate groups, OAM and OAH, respectively). We employed free-energy perturbation (FEP) for relative binding energies at the molecular mechanics (MM) and the combined quantum mechanical (QM) and MM (QM/MM) levels, the latter obtained with the reference-potential approach with QM/MM sampling for the MM → QM/MM FEP. The semiempirical QM method PM6-DH+ was employed for the ligand in the latter calculations. Moreover, binding free energies were also estimated from QM/MM optimised structures, combined with COSMO-RS estimates of the solvation energy and thermostatistical corrections from MM frequencies. They were performed at the PM6-DH+ level of theory with the full host and guest molecule in the QM system (and also four water molecules in the geometry optimisations) for 10-20 snapshots from molecular dynamics simulations of the complex. Finally, the structure with the lowest free energy was recalculated using the dispersion-corrected density-functional theory method TPSS-D3, for both the structure and the energy. The two FEP approaches gave similar results (PM6-DH+/MM slightly better for OAM), which were among the five submissions with the best performance in the challenge and gave the best results without any fit to data from the SAMPL5 challenge, with mean absolute deviations (MAD) of 2.4-5.2 kJ/mol and a correlation coefficient (R 2 ) of 0.77-0.93. This is the first time QM/MM approaches give binding free energies that are competitive to those obtained with MM for the octa-acid host. The QM/MM-optimised structures gave somewhat worse performance (MAD = 3-8 kJ/mol and R 2  = 0.1-0.9), but the results were improved compared to previous studies of this system with similar methods.

The DMSO reductase family is the largest and most diverse family of mononuclear molybdenum oxygen-atom-transfer proteins. Their active sites contain a Mo ion coordinated to two molybdopterin ligands, one oxo group in the oxidised state, and one additional, often protein-derived ligand. We have used density-functional theory to evaluate how the fourth ligand (serine, cysteine, selenocysteine, OH - , O 2- , SH - , or S 2- ) affects the geometries, reaction mechanism, reaction energies, and reduction potentials of intermediates in the DMSO reductase reaction. Our results show that there are only small changes in the geometries of the reactant and product states, except from the elongation of the MoX bond as the ionic radius of XO, S, Se increases. The five ligands with a single negative charge gave an identical two-step reaction mechanism, in which DMSO first binds to the reduced active site, after which the SO bond is cleaved, concomitantly with the transfer of two electrons from Mo in a rate-determining second transition state. The five models gave similar activation energies of 69-85kJ/mol, with SH - giving the lowest barrier. In contrast, the O 2- and S 2- ligands gave much higher activation energies (212 and 168kJ/mol) and differing mechanisms (a more symmetric intermediate for O 2- and a one-step reaction without any intermediate for S 2- ). The high activation energies are caused by a less exothermic reaction energy, 13-25kJ/mol, and by a more stable reactant state owing to the strong MoO 2- or MoS 2- bonds., (Copyright © 2017 Elsevier Inc. All rights reserved.)

We have tried to calculate the free energy for the binding of six small ligands to two variants of the octa-acid deep cavitand host in the SAMPL5 blind challenge. We employed structures minimised with dispersion-corrected density-functional theory with small basis sets and energies were calculated using large basis sets. Solvation energies were calculated with continuum methods and thermostatistical corrections were obtained from frequencies calculated at the HF-3c level. Care was taken to minimise the effects of the flexibility of the host by keeping the complexes as symmetric and similar as possible. In some calculations, the large net charge of the host was reduced by removing the propionate and benzoate groups. In addition, the effect of a restricted molecular dynamics sampling of structures was tested. Finally, we tried to improve the energies by using the DLPNO-CCSD(T) approach. Unfortunately, results of quite poor quality were obtained, with no correlation to the experimental data, systematically too positive affinities (by ~50 kJ/mol) and a mean absolute error (after removal of the systematic error) of 11-16 kJ/mol. DLPNO-CCSD(T) did not improve the results, so the accuracy is not limited by the energy function. Instead, four likely sources of errors were identified: first, the minimised structures were often incorrect, owing to the omission of explicit solvent. They could be partly improved by performing the minimisations in a continuum solvent with four water molecules around the charged groups of the ligands. Second, some ligands could bind in several different conformations, requiring sampling of reasonable structures. Third, there is an indication the continuum-solvation model has problems to accurately describe the binding of both the negatively and positively charged guest molecules. Fourth, different methods to calculate the thermostatistical corrections gave results that differed by up to 30 kJ/mol and there is an indication that HF-3c overestimates the entropy term. In conclusion, it is a challenge to calculate binding affinities for this octa-acid system with quantum-mechanical methods.

We have examined how the self-interaction error in density-functional theory (DFT) calculations affects energies calculated on large systems (600-1000 atoms) involving several charged groups. We employ 18 different quantum mechanical (QM) methods, including Hartree-Fock, as well as pure, hybrid, and range-separated DFT methods. They are used to calculate reaction and activation energies for three different protein models in vacuum, in a point-charge surrounding, or with a continuum-solvent model. We show that pure DFT functionals give rise to a significant delocalization of the charges in charged groups in the protein, typically by ∼0.1 e, as evidenced from the Mulliken charges. This has a clear effect on how the surroundings affect calculated reaction and activation energies, indicating that these methods should be avoided for DFT calculations on large systems. Fortunately, methods such as CAM-B3LYP, BHLYP, and M06-2X give results that agree within a few kilojoules per mole, especially when the calculations are performed in a point-charge surrounding. Therefore, we recommend these methods to estimate the effect of the surroundings with large QM systems (but other QM methods may be used to study the intrinsic reaction and activation energies).

In this chapter, I discuss combined quantum mechanics (QM) and molecular mechanics (MM; QM/MM) calculations for proteins. In QM/MM, a small but interesting part of the protein is treated by accurate QM methods, whereas the remainder is treated by faster MM methods. The prime problems with QM/MM calculations are bonds between the QM and MM systems, the selection of the QM system, and the local-minima problem. The two first problems can be solved by the big-QM approach, including in the QM calculation all groups within 4.5-6Å of the active site and all buried charges in the protein. The third problem can be solved by calculating free energies. It is important to study QM/MM energy components to ensure that the results are stable and reliable. They can also be used to understand the reaction and the effect of the surroundings, eg, by dividing the catalytic effect into bonded, van der Waals, electrostatic, and geometric components and to deduce which parts of the protein contribute most to the catalysis. It should be ensured that the QM calculations are reliable and converged by extending the basis set to quadruple-zeta quality, including a proper treatment of dispersion, as well as years experience and method development calculations with both pure and hybrid density functional theory methods. If the latter give differing results, calibration with high-level QM methods is needed. Reactions that change the net charge should be avoided. QM/MM calculations can be combined with experimental methods., (© 2016 Elsevier Inc. All rights reserved.)

We have performed a large-scale test of alchemical perturbation calculations with the Bennett acceptance-ratio (BAR) approach to estimate relative affinities for the binding of 107 ligands to 10 different proteins. Employing 20-Å truncated spherical systems and only one intermediate state in the perturbations, we obtain an error of less than 4 kJ/mol for 54% of the studied relative affinities and a precision of 0.5 kJ/mol on average. However, only four of the proteins gave acceptable errors, correlations, and rankings. The results could be improved by using nine intermediate states in the simulations or including the entire protein in the simulations using periodic boundary conditions. However, 27 of the calculated affinities still gave errors of more than 4 kJ/mol, and for three of the proteins the results were not satisfactory. This shows that the performance of BAR calculations depends on the target protein and that several transformations gave poor results owing to limitations in the molecular-mechanics force field or the restricted sampling possible within a reasonable simulation time. Still, the BAR results are better than docking calculations for most of the proteins.

Molecular mechanics with generalised Born and surface area solvation (MM/GBSA) is a popular method to calculate the free energy of the binding of ligands to proteins. It involves molecular dynamics (MD) simulations with an explicit solvent of the protein-ligand complex to give a set of snapshots for which energies are calculated with an implicit solvent. This change in the solvation method (explicit → implicit) would strictly require that the energies are reweighted with the implicit-solvent energies, which is normally not done. In this paper we calculate MM/GBSA energies with two generalised Born models for snapshots generated by the same methods or by explicit-solvent simulations for five synthetic N-acetyllactosamine derivatives binding to galectin-3. We show that the resulting energies are very different both in absolute and relative terms, showing that the change in the solvent model is far from innocent and that standard MM/GBSA is not a consistent method. The ensembles generated with the various solvent models are quite different with root-mean-square deviations of 1.2-1.4 Å. The ensembles can be converted to each other by performing short MD simulations with the new method, but the convergence is slow, showing mean absolute differences in the calculated energies of 6-7 kJ mol(-1) after 2 ps simulations. Minimisations show even slower convergence and there are strong indications that the energies obtained from minimised structures are different from those obtained by MD.

We have performed a systematic study of the entropy term in the MM/GBSA (molecular mechanics combined with generalized Born and surface-area solvation) approach to calculate ligand-binding affinities. The entropies are calculated by a normal-mode analysis of harmonic frequencies from minimized snapshots of molecular dynamics simulations. For computational reasons, these calculations have normally been performed on truncated systems. We have studied the binding of eight inhibitors of blood clotting factor Xa, nine ligands of ferritin, and two ligands of HIV-1 protease and show that removing protein residues with distances larger than 8-16 Å to the ligand, including a 4 Å shell of fixed protein residues and water molecules, change the absolute entropies by 1-5 kJ/mol on average. However, the change is systematic, so relative entropies for different ligands change by only 0.7-1.6 kJ/mol on average. Consequently, entropies from truncated systems give relative binding affinities that are identical to those obtained for the whole protein within statistical uncertainty (1-2 kJ/mol). We have also tested to use a distance-dependent dielectric constant in the minimization and frequency calculation (ε = 4r), but it typically gives slightly different entropies and poorer binding affinities. Therefore, we recommend entropies calculated with the smallest truncation radius (8 Å) and ε =1. Such an approach also gives an improved precision for the calculated binding free energies.

We show that conformational entropies calculated for five proteins and protein-ligand complexes with dihedral-distribution histogramming, the von Mises approach, or quasi-harmonic analysis do not converge to any useful precision even if molecular dynamics (MD) simulations of 380-500 ns length are employed (the uncertainty is 12-89 kJ mol(-1)). To explain this, we suggest a simple protein model involving dihedrals with effective barriers forming a uniform distribution and show that for such a model, the entropy increases logarithmically with time until all significantly populated dihedral states have been sampled, in agreement with the simulations (during the simulations, 52-70% of the available dihedral phase space has been visited). This is also confirmed by the analysis of the trajectories of a 1 ms simulation of bovine pancreatic trypsin inhibitor (31 kJ mol(-1) difference in the entropy between the first and second part of the simulation). Strictly speaking, this means that it is practically impossible to equilibrate MD simulations of proteins. We discuss the implications of such a lack of strict equilibration of protein MD simulations and show that ligand-binding free energies estimated with the MM/GBSA method (molecular mechanics with generalised Born and surface-area solvation) vary by 3-15 kJ mol(-1) during a 500 ns simulation (the higher estimate is caused by rare conformational changes), although they involve a questionable but well-converged normal-mode entropy estimate, whereas free energies estimated by free-energy perturbation vary by less than 0.6 kJ mol(-1) for the same simulation.

We present a combination of semiempirical quantum-mechanical (SQM) calculations in the conductor-like screening model with the MM/GBSA (molecular-mechanics with generalized Born and surface-area solvation) method for ligand-binding affinity calculations. We test three SQM Hamiltonians, AM1, RM1, and PM6, as well as hydrogen-bond corrections and two different dispersion corrections. As test cases, we use the binding of seven biotin analogues to avidin, nine inhibitors to factor Xa, and nine phenol-derivatives to ferritin. The results vary somewhat for the three test cases, but a dispersion correction is mandatory to reproduce experimental estimates. On average, AM1 with the DH2 hydrogen-bond and dispersion corrections gives the best results, which are similar to those of standard MM/GBSA calculations for the same systems. The total time consumption is only 1.3-1.6 times larger than for MM/GBSA., (Copyright © 2012 Wiley Periodicals, Inc.)

We have compared the predictions of ligand-binding affinities from several methods based on end-point molecular dynamics simulations and continuum solvation, that is, methods related to MM/PBSA (molecular mechanics combined with Poisson-Boltzmann and surface area solvation). Two continuum-solvation models were considered, viz., the Poisson-Boltzmann (PB) and generalised Born (GB) approaches. The nonelectrostatic energies were also obtained in two different ways, viz., either from the sum of the bonded, van der Waals, nonpolar solvation energies, and entropy terms (as in MM/PBSA), or from the scaled protein-ligand van der Waals interaction energy (as in the linear interaction energy approach, LIE). Three different approaches to calculate electrostatic energies were tested, viz., the sum of electrostatic interaction energies and polar solvation energies, obtained either from a single simulation of the complex or from three independent simulations of the complex, the free protein, and the free ligand, or the linear-response approximation (LRA). Moreover, we investigated the effect of scaling the electrostatic interactions by an effective internal dielectric constant of the protein (ϵ(int) ). All these methods were tested on the binding of seven biotin analogues to avidin and nine 3-amidinobenzyl-1H-indole-2-carboxamide inhibitors to factor Xa. For avidin, the best results were obtained with a combination of the LIE nonelectrostatic energies with the MM+GB electrostatic energies from a single simulation, using ϵ(int) = 4. For fXa, standard MM/GBSA, based on one simulation and using ϵ(int) = 4-10 gave the best result. The optimum internal dielectric constant seems to be slightly higher with PB than with GB solvation., (Copyright © 2012 Wiley Periodicals, Inc.)

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

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.

A common approach for the computational modeling of enzyme reactions is to study a rather small model of the active site (20-200 atoms) with quantum mechanical (QM) methods, modeling the rest of the surroundings by a featureless continuum with a dielectric constant of approximately 4. In this paper, we discuss how the residues included in the QM model should be selected and how many residues need to be included before reaction energies converge. As a test case, we use a proton-transfer reaction between a first-sphere cysteine ligand and a second-sphere histidine group in the active site of [Ni,Fe] hydrogenase. We show that it is not a good approach to add groups according to their distance to the active site. A better approach is to add groups according to their contributions to the QM/MM energy difference. However, the energies can still vary by up to 50 kJ/mol for QM systems of sizes up to 230 atoms. In fact, the QM-only approach is based on the hope that a large number of sizable contributions will cancel. Interactions with neutral groups are, in general, short-ranged, with net energy contributions of less than 4 kJ/mol at distances above 5 A from the active site. Interactions with charged groups are much more long-ranged, and interactions with buried charges 20 A from the active site can still contribute by 5 kJ/mol to the reaction energy. Thus, to accurately model the influence of the surroundings on enzyme reaction energies, a detailed and unbiased atomistic account of the surroundings needs to be included.

We have developed a method to estimate accurate interaction energies between a full protein and a bound ligand. It is based on the recently proposed PMISP (polarizable multipole interaction with supermolecular pairs) method (Soderhjelm, P.; Ryde, U. J. Phys. Chem. A 2009, 113, 617), which treats electrostatic interaction by multipoles up to quadrupoles, induction by anisotropic polarizabilities, and nonclassical interactions by explicit quantum mechanical (QM) calculations, using a fragmentation approach. For a whole protein, electrostatics and induction are treated the same way, but for the nonclassical interactions, a Lennard-Jones term from a standard molecular mechanics (MM) force field (e.g., Amber) is used outside a certain distance from the ligand (4-7 A). This QM/MM variant of the PMISP method is carefully tested by varying this distance. Several approximations related to the classical interactions are also evaluated. It is found that one can speed up the calculation by using density functional theory to compute multipoles and polarizabilities but that a proper treatment of polarization is important. As a demonstration of the method, the interaction energies of two ligands bound to avidin are calculated at the MP2/aug-cc-pVTZ level, with an expected relative error of 1-2%.

We have developed a method to estimate free energies of reactions in proteins, called QM/MM-PBSA. It estimates the internal energy of the reactive site by quantum mechanical (QM) calculations, whereas bonded, electrostatic, and van der Waals interactions with the surrounding protein are calculated at the molecular mechanics (MM) level. The electrostatic part of the solvation energy of the reactant and the product is estimated by solving the Poisson-Boltzmann (PB) equation, and the nonpolar part of the solvation energy is estimated from the change in solvent-accessible surface area (SA). Finally, the change in entropy is estimated from the vibrational frequencies. We test this method for five proton-transfer reactions in the active sites of [Ni,Fe] hydrogenase and copper nitrite reductase. We show that QM/MM-PBSA reproduces the results of a strict QM/MM free-energy perturbation method with a mean absolute deviation (MAD) of 8-10 kJ/mol if snapshots from molecular dynamics simulations are used and 4-14 kJ/mol if a single QM/MM structure is used. This is appreciably better than the original QM/MM results or if the QM energies are supplemented with a point-charge model, a self-consistent reaction field, or a PB model of the protein and the solvent, which give MADs of 22-36 kJ/mol for the same test set.

The four-electron reduction of dioxygen to water in multicopper oxidases takes place in a trinuclear copper cluster, which is linked to a mononuclear blue copper site, where the substrates are oxidized. Recently, several intermediates in the catalytic cycle have been spectroscopically characterized, and two possible structural models have been suggested for both the peroxy and native intermediates. In this study, these spectroscopic results are complemented by hybrid quantum and molecular mechanical (QM/MM) calculations, taking advantage of recently available crystal structures with a full complement of copper ions. Thereby, we obtain optimized molecular structures for all of the experimentally studied intermediates involved in the reductive cleavage of the O(2) molecule and energy profiles for individual reaction steps. This allows identification of the experimentally observed intermediates and further insight into the reaction mechanism that is probably relevant for the whole class of multicopper oxidases. We suggest that the peroxy intermediate contains an O(2)(2-) ion, in which one oxygen atom bridges the type 2 copper ion and one of the type 3 copper ions, whereas the other one coordinates to the other type 3 copper ion. One-electron reduction of this intermediate triggers the cleavage of the O-O bond, which involves the uptake of a proton. The product of this cleavage is the observed native intermediate, which we suggest to contain a O(2-) ion coordinated to all three of the copper ions in the center of the cluster.

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.

A method is developed for the combination of quantum chemical geometry optimizations and crystallographic structure refinement. The method is implemented by integrating the quantum chemical software Turbomole with the crystallographic software Crystallography and NMR System (CNS), using three small procedures transferring information between the two programs. The program (COMQUM-X)is used to study the binding of the inhibitor N-methylmesoporphyrin to ferrochelatase, and we show that the method behaves properly and leads to an improvement of the structure of the inhibitor. It allows us to directly quantify in energy terms how much the protein distort the structure of the bound inhibitor compared to the optimum vacuum structure (4-6 kJ/mol). The approach improves the standard combined quantum chemical and molecular mechanics (QC/MM) approach by guaranteeing that the final structure is in accordance with experimental data (the reflections) and avoiding the risk of propagating errors in the crystal coordinates. The program can also be seen as an improvement of standard crystallographic refinement, providing an accurate empirical potential function for any group of interest. The results can be directly interpreted in standard crystallographic terms (e.g., R factors or electron density maps). The method can be used to interpret crystal structures (e.g., the protonation status of metal-bound water molecules) and even to locally improve them., (Copyright 2002 Wiley Periodicals, Inc. J Comput Chem 23: 1058-1070, 2002)

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.