Proton Transfer at Metal Sites in Proteins Studied by Quantum Mechanical Free-Energy Perturbations

Catalytic Reaction Mechanism of Glyoxalase II: A Quantum Mechanics/Molecular Mechanics Study

Methylglyoxal (MG) is a reactive and toxic compound produced in carbohydrate, lipid, and amino acid metabolism. The glyoxalase system is the main detoxifying route for MG and consists of two enzymes, glyoxalase I (GlxI) and glyoxalase II (GlxII). GlxI catalyzes the formation of S-d-lactoylglutathione from hemithioacetal, and GlxII converts this intermediate to d-lactate. A relationship between the glyoxalase system and some diseases like diabetes has been shown, and inhibiting enzymes of this system may be an effective means of controlling certain diseases. A detailed understanding of the reaction mechanism of an enzyme is essential to the rational design of competitive inhibitors. In this work, we use quantum mechanics/molecular mechanics (QM/MM) calculations and energy refinement utilizing the big-QM and QM/MM thermodynamic cycle perturbation methods to propose a mechanism for the GlxII reaction that starts with a nucleophilic attack of the bridging OH^- group on the substrate. The coordination of the substrate to the Zn ions places its electrophilic center close to the hydroxide group, enabling the reaction to proceed. Our estimated reaction energies are in excellent agreement with experimental data, thus demonstrating the reliability of our approach and the proposed mechanism. Additionally, we examined alternative protonation states of Asp-29, Asp-58, Asp-134, and the bridging hydroxide ion in the catalytic process. However, these give less favorable reactions, a poorer reproduction of the crystal structure geometry of the active site, and higher root-mean-squared deviations of the active site residues in molecular dynamics simulations.

Metalloprotein catalysis: structural and mechanistic insights into oxidoreductases from neutron protein crystallography

Metalloproteins catalyze a range of reactions, with enhanced chemical functionality due to their metal cofactor. The reaction mechanisms of metalloproteins have been experimentally characterized by spectroscopy, macromolecular crystallography and cryo-electron microscopy. An important caveat in structural studies of metalloproteins remains the artefacts that can be introduced by radiation damage. Photoreduction, radiolysis and ionization deriving from the electromagnetic beam used to probe the structure complicate structural and mechanistic interpretation. Neutron protein diffraction remains the only structural probe that leaves protein samples devoid of radiation damage, even when data are collected at room temperature. Additionally, neutron protein crystallography provides information on the positions of light atoms such as hydrogen and deuterium, allowing the characterization of protonation states and hydrogen-bonding networks. Neutron protein crystallography has further been used in conjunction with experimental and computational techniques to gain insight into the structures and reaction mechanisms of several transition-state metal oxidoreductases with iron, copper and manganese cofactors. Here, the contribution of neutron protein crystallography towards elucidating the reaction mechanism of metalloproteins is reviewed.

Two-Substrate Glyoxalase I Mechanism: A Quantum Mechanics/Molecular Mechanics Study

Glyoxalase I (GlxI) is an important enzyme that catalyzes the detoxification of methylglyoxal (MG) with the help of glutathione (H-SG). It is currently unclear whether MG and H-SG are substrates of GlxI or whether the enzyme processes hemithioacetal (HTA), which is nonenzymatically formed from MG and H-SG. Most previous studies have concentrated on the latter mechanism. Here, we study the two-substrate reaction mechanism of GlxI from humans (HuGlxI) and corn (ZmGlxI), which are Zn(II)-active and -inactive, respectively. Hybrid quantum mechanics/molecular mechanics calculations were used to obtain geometrical structures of the stationary points along reaction paths, and big quantum mechanical systems with more than 1000 atoms and free-energy perturbations were used to improve the quality of the calculated energies. We studied, on an equal footing, all reasonable reaction paths to the S- and R-enantiomers of HTA from MG and H-SG (the latter was considered in two different binding modes). The results indicate that the MG and H-SG reaction in both enzymes can follow the same path to reach S-HTA. However, the respective overall barriers and reaction energies are different for the two enzymes (6.1 and -9.8 kcal/mol for HuGlxI and 15.7 and -2.2 kcal/mol for ZmGlxI). The first reaction step to produce S-HTA is facilitated by a crystal water molecule that forms hydrogen bonds with a Glu and a Thr residue in the active site. The two enzymes also follow similar paths to R-HTA. However, the reactions reach a deprotonated and protonated R-HTA in the human and corn enzymes, respectively. The production of deprotonated R-HTA in HuGlxI is consistent with other theoretical and experimental works. However, our calculations show a different behavior for ZmGlxI (both S- and R-HTA can be formed in the enzyme with the alcoholic proton on HTA). This implies that Glu-144 of corn GlxI is not basic enough to keep the alcoholic proton. In HuGlxI, the two binding modes of H-SG that lead to S- and R-HTA are degenerate, but the barrier leading to R-HTA is lower than the barrier to S-HTA. On the other hand, ZmGlxI prefers the binding mode, which produces S-HTA; this observation is consistent with experiments. Based on the results, we present a modification for a previously proposed two-substrate reaction mechanism for ZmGlxI.

Functional role of an unusual tyrosine residue in the electron transfer chain of a prokaryotic (6-4) photolyase

Cryptochromes and photolyases form a flavoprotein family in which the FAD chromophore undergoes light induced changes of its redox state. During this process, termed photoreduction, electrons flow from the surface via conserved amino acid residues to FAD. The bacterial (6-4) photolyase PhrB belongs to a phylogenetically ancient group. Photoreduction of PhrB differs from the typical pattern because the amino acid of the electron cascade next to FAD is a tyrosine (Tyr391), whereas photolyases and cryptochromes of other groups have a tryptophan as direct electron donor of FAD. Mutagenesis studies have identified Trp342 and Trp390 as essential for charge transfer. Trp342 is located at the periphery of PhrB while Trp390 connects Trp342 and Tyr391. The role of Tyr391, which lies between Trp390 and FAD, is however unclear as its replacement by phenylalanine did not block photoreduction. Experiments reported here, which replace Tyr391 by Ala, show that photoreduction is blocked, underlining the relevance of Tyr/Phe at position 391 and indicating that charge transfer occurs via the triad 391-390-342. This raises the question, why PhrB positions a tyrosine at this location, having a less favourable ionisation potential than tryptophan, which occurs at this position in many proteins of the photolyase/cryptochrome family. Tunnelling matrix calculations show that tyrosine or phenylalanine can be involved in a productive bridged electron transfer between FAD and Trp390, in line with experimental findings. Since replacement of Tyr391 by Trp resulted in loss of FAD and DMRL chromophores, electron transfer cannot be studied experimentally in this mutant, but calculations on a mutant model suggest that Trp might participate in the electron transfer cascade. Charge transfer simulations reveal an unusual stabilization of the positive charge on site 391 compared to other photolyases or cryptochromes. Water molecules near Tyr391 offer a polar environment which stabilizes the positive charge on this site, thereby lowering the energetic barrier intrinsic to tyrosine. This opens a second charge transfer channel in addition to tunnelling through the tyrosine barrier, based on hopping and therefore transient oxidation of Tyr391, which enables a fast charge transfer similar to proteins utilizing a tryptophan-triad. Our results suggest that evolution of the first site of the redox chain has just been possible by tuning the protein structure and environment to manage a downhill hole transfer process from FAD to solvent.

Influence of Coupling and Embedding Schemes on QM Size Convergence in QM/MM Approaches for the Example of a Proton Transfer in DNA

The influence of embedding and coupling schemes on the convergence of the QM size in the QM/MM approach is investigated for the transfer of a proton in a DNA base pair. We find that the embedding scheme (mechanical or electrostatic) has a much greater impact on the convergence behavior than the coupling scheme (additive QM/MM or subtractive ONIOM). To achieve size convergence, QM regions with up to 6000 atoms are necessary for pure QM or mechanical embedding. In contrast, electrostatic embedding converges faster: for the example of the transfer of a proton between DNA base pairs, we recommend including at least five base pairs and 5 Å of solvent (including counterions) into the QM region, i.e., a total of 1150 atoms.

Human Ferrochelatase: Insights for the Mechanism of Ferrous Iron Approaching Protoporphyrin IX by QM/MM and QTCP Free Energy Studies

Ferrochelatase catalyzes the insertion of ferrous iron into protoporphyrin IX, the terminal step in heme biosynthesis. Some disputes in its mechanism remain unsolved, especially for human ferrochelatase. In this paper, high-level quantum mechanical/molecular mechanics (QM/MM) and free-energy studies were performed to address these controversial issues including the iron-binding site, the optimal reaction path, the substrate porphyrin distortion, and the presence of the sitting-atop (SAT) complex. Our results reveal that the ferrous iron is probably at the binding site coordinating with Met76, and His263 plays the role of proton acceptor. The rate-determining step is either the first proton removed by His263 or the proton transition within the porphyrin with an energy barrier of 14.99 or 14.87 kcal/mol by the quantum mechanical thermodynamic cycle perturbation (QTCP) calculations, respectively. The fast deprotonation step with the conservative residues rather than porphyrin deformation found in solution provides the driving force for biochelation. The SAT complex is not a necessity for the catalysis though it induces a modest distortion on the porphyrin ring.

Efficiently computing pathway free energies: New approaches based on chain-of-replica and Non-Boltzmann Bennett reweighting schemes

BACKGROUND

Accurately modeling condensed phase processes is one of computation's most difficult challenges. Include the possibility that conformational dynamics may be coupled to chemical reactions, where multiscale (i.e., QM/MM) methods are needed, and this task becomes even more daunting.

METHODS

Free energy simulations (i.e., molecular dynamics), multiscale modeling, and reweighting schemes.

RESULTS

Herein, we present two new approaches for mitigating the aforementioned challenges. The first is a new chain-of-replica method (off-path simulations, OPS) for computing potentials of mean force (PMFs) along an easily defined reaction coordinate. This development is coupled with a new distributed, highly-parallel replica framework (REPDstr) within the CHARMM package. Validation of these new schemes is carried out on two processes that undergo conformational changes. First is the simple torsional rotation of butane, while a much more challenging glycosidic rotation (in vacuo and solvated) is the second. Additionally, a new approach that greatly improves (i.e., possibly an order of magnitude) the efficiency of computing QM/MM PMFs is introduced and compared to standard schemes. Our efforts are grounded in the recently developed method for efficiently computing QM-based free energies (i.e., QM-Non-Boltzmann Bennett, QM-NBB). Again, we validate this new technique by computing the QM/MM PMF of butane's torsional rotation.

CONCLUSIONS

The OPS-REPDstr method is a promising new approach that overcomes many limitations of standard pathway simulations in CHARMM. The combination of QM-NBB with pathway techniques is very promising as it offers significant advantages over current procedures.

GENERAL SIGNIFICANCE

Efficiently computing potentials of mean force is a major, unresolved, area of interest. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

Evidence for substrate binding-induced zwitterion formation in the catalytic Cys-His dyad of the SARS-CoV main protease

The coronavirus main protease (M(pro)) represents an attractive drug target for antiviral therapy of coronavirus (CoV) infections, including severe acute respiratory syndrome (SARS). The SARS-CoV M(pro) and related CoV proteases have several distinct features, such as an uncharged Cys-His catalytic dyad embedded in a chymotrypsin-like protease fold, that clearly separate these enzymes from archetypical cysteine proteases. To further characterize the catalytic system of CoV main proteases and to obtain information about improved inhibitors, we performed comprehensive simulations of the proton-transfer reactions in the SARS-CoV M(pro) active site that lead to the Cys(-)/His(+) zwitterionic state required for efficient proteolytic activity. Our simulations, comprising the free enzyme as well as substrate-enzyme and inhibitor-enzyme complexes, lead us to predict that zwitterion formation is fostered by substrate binding but not inhibitor binding. This indicates that M(pro) employs a substrate-induced catalytic mechanism that further enhances its substrate specificity. Our computational data are in line with available experimental results, such as X-ray geometries, measured pKa values, mutagenesis experiments, and the measured differences between the kinetic parameters of substrates and inhibitors. The data also provide an atomistic picture of the formerly postulated electrostatic trigger involved in SARS-CoV M(pro) activity. Finally, they provide information on how a specific microenvironment may finely tune the activity of M(pro) toward specific viral protein substrates, which is known to be required for efficient viral replication. Our simulations also indicate that the low inhibition potencies of known covalently interacting inhibitors may, at least in part, be attributed to insufficient fostering of the proton-transfer reaction. These findings suggest ways to achieve improved inhibitors.

Recent advances in QM/MM free energy calculations using reference potentials

Background

Recent years have seen enormous progress in the development of methods for modeling (bio)molecular systems. This has allowed for the simulation of ever larger and more complex systems. However, as such complexity increases, the requirements needed for these models to be accurate and physically meaningful become more and more difficult to fulfill. The use of simplified models to describe complex biological systems has long been shown to be an effective way to overcome some of the limitations associated with this computational cost in a rational way.

Scope of review

Hybrid QM/MM approaches have rapidly become one of the most popular computational tools for studying chemical reactivity in biomolecular systems. However, the high cost involved in performing high-level QM calculations has limited the applicability of these approaches when calculating free energies of chemical processes. In this review, we present some of the advances in using reference potentials and mean field approximations to accelerate high-level QM/MM calculations. We present illustrative applications of these approaches and discuss challenges and future perspectives for the field.

Major conclusions

The use of physically-based simplifications has shown to effectively reduce the cost of high-level QM/MM calculations. In particular, lower-level reference potentials enable one to reduce the cost of expensive free energy calculations, thus expanding the scope of problems that can be addressed.

General significance

As was already demonstrated 40 years ago, the usage of simplified models still allows one to obtain cutting edge results with substantially reduced computational cost. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

•

We present some of the advances to accelerate high-level QM/MM calculations.

•

Quantitative limitations of low-level methods can be overcome by these approaches.

•

Reference potentials make free energy simulations feasible for large systems.

•

Automated fitting reduces the need of expensive sampling of high-level approaches.

•

Application of reference potentials can be extended to a wide range of processes.

Benchmark Study for the Cysteine-Histidine Proton Transfer Reaction in a Protein Environment: Gas Phase, COSMO, QM/MM Approaches

Proton transfer reactions are of crucial interest for the investigation of proteins. We have investigated the accuracy of commonly used quantum chemical methods for the description of proton transfer reactions in different environments (gas phase, COSMO, QM/MM) using the proton transfer between the catalytic dyad residues cysteine 145 and histidine 41 of SARS coronavirus main protease as a case study. The test includes thermodynamic, kinetic, and structural properties. The study comprises computationally demanding ab initio approaches (HF, CC2, MP2, SCS-CC2, SCS-MP2, CCSD(T)), popular density functional theories (BLYP, B3LYP, M06-2X), and semiempirical methods (MNDO/d, AM1, RM1, PM3, PM6). The approximated coupled cluster approach LCCSD(T) is taken as a reference method. We find that the robustness of the tested methods with respect to the environment correlates well with the level of theory. As an example HF, CC2, MP2, and their SCS variants show similar errors for gas phase, COSMO, or QM/MM computations. In contrast for semiempirical methods, the errors strongly diversify if one goes from gas phase to COSMO or QM/MM. Particular problems are observed for the recent semiempirical methods PM6 and RM1, which show the best performance for gas phase calculations but possess larger errors in conjunction with COSMO. Finally, a combination of SCS-MP2 and B3LYP or M06-2X allows reliable estimates about remaining errors.

Accurate Reaction Energies in Proteins Obtained by Combining QM/MM and Large QM Calculations

We here suggest and test a new method to obtain stable energies in proteins for charge-neutral reactions by running large quantum mechanical (QM) calculations on structures obtained by combined QM and molecular mechanics (QM/MM) geometry optimization on several snapshots from molecular dynamics simulations. As a test case, we use a proton transfer between a metal-bound cysteine residue and a second-sphere histidine residue in the active site of [Ni,Fe] hydrogenase, which has been shown to be very sensitive to the surroundings. We include in the QM calculations all residues within 4.5 Å of the active site, two capped residues on each side of the active-site residues, and all charged groups that are buried inside the protein, which for this enzyme includes three iron-sulfur clusters, in total, 930 atoms. These calculations are performed at the BP86/def2-SV(P) level, but the energies are then extrapolated to the B3LYP/def2-TZVP level with a smaller QM system, and zero-point energy, entropy, and thermal effects are added. We test three approaches to model the remaining atoms of the protein solvent, viz., by standard QM/MM approaches using either mechanical or electrostatic embedding or by using a continuum solvation model for the large QM systems. Quite encouragingly, the three approaches give the same results within 14 kJ/mol, and variations in the size of the QM system do not change the energies by more than 8 kJ/mol, provided that the QM/MM junctions are not moved closer to the QM system. The statistical precision for the average over 10 snapshots is 1-3 kJ/mol.

Convergence of QM/MM free-energy perturbations based on molecular-mechanics or semiempirical simulations

Lately, there has been great interest in performing free-energy perturbation (FEP) at the combined quantum mechanics and molecular mechanics (QM/MM) level, e.g. for enzyme reactions. Such calculations require extensive sampling of phase space, which typically is prohibitive with density-functional theory or ab initio methods. Therefore, such calculations have mostly been performed with semiempirical QM (SQM) methods, or by using a thermodynamic cycle involving sampling at the MM level and perturbations between the MM and QM/MM levels of theory. However, the latter perturbations typically have convergence problems, unless the QM system is kept fixed during the simulations, because the MM and QM/MM descriptions of the internal degrees of freedom inside the QM system are too dissimilar. We have studied whether the convergence of the MM → QM/MM perturbation can be improved by using a thoroughly parameterised force field or by using SQM/MM methods. As a test case we use the first half-reaction of haloalkane dehalogenase and the QM calculations are performed with the PBE, B3LYP, and TPSSH density-functional methods. We show that the convergence can be improved with a tailored force field, but only locally around the parameterised state. Simulations based on SQM/MM methods using the MNDO, AM1, PM3, RM1, PDDG-MNDO, and PDDG-PM3 Hamiltonians have slightly better convergence properties, but very long simulations are still needed (~10 ns) and convergence is obtained only if electrostatic interactions between the QM system and the surroundings are ignored. This casts some doubts on the common practice to base QM/MM FEPs on semiempirical simulations without any reweighting of the trajectories.

Reorganization energy for internal electron transfer in multicopper oxidases

We have calculated the reorganization energy for the intramolecular electron transfer between the reduced type 1 copper site and the peroxy intermediate of the trinuclear cluster in the multicopper oxidase CueO. The calculations are performed at the combined quantum mechanics and molecular mechanics (QM/MM) level, based on molecular dynamics simulations with tailored potentials for the two copper sites. We obtain a reorganization energy of 91-133 kJ/mol, depending on the theoretical treatment. The two Cu sites contribute by 12 and 22 kJ/mol to this energy, whereas the solvent contribution is 34 kJ/mol. The rest comes from the protein, involving small contributions from many residues. We have also estimated the energy difference between the two electron-transfer states and show that the reduction of the peroxy intermediate is exergonic by 43-87 kJ/mol, depending on the theoretical method. Both the solvent and the protein contribute to this energy difference, especially charged residues close to the two Cu sites. We compare these estimates with energies obtained from QM/MM optimizations and QM calculations in a vacuum and discuss differences between the results obtained at various levels of theory.

Reduction potentials and acidity constants of Mn superoxide dismutase calculated by QM/MM free-energy methods

We used two theoretical methods to estimate reduction potentials and acidity constants in Mn superoxide dismutase (MnSOD), namely combined quantum mechanical and molecular mechanics (QM/MM) thermodynamic cycle perturbation (QTCP) and the QM/MM-PBSA approach. In the latter, QM/MM energies are combined with continuum solvation energies calculated by solving the Poisson-Boltzmann equation (PB) or by the generalised Born approach (GB) and non-polar solvation energies calculated from the solvent-exposed surface area. We show that using the QTCP method, we can obtain accurate and precise estimates of the proton-coupled reduction potential for MnSOD, 0.30±0.01 V, which compares favourably with experimental estimates of 0.26-0.40 V. However, the calculated potentials depend strongly on the DFT functional used: The B3LYP functional gives 0.6 V more positive potentials than the PBE functional. The QM/MM-PBSA approach leads to somewhat too high reduction potentials for the coupled reaction and the results depend on the solvation model used. For reactions involving a change in the net charge of the metal site, the corresponding results differ by up to 1.3 V or 24 pK(a) units, rendering the QM/MM-PBSA method useless to determine absolute potentials. However, it may still be useful to estimate relative shifts, although the QTCP method is expected to be more accurate.

Formation of the iron-oxo hydroxylating species in the catalytic cycle of aromatic amino acid hydroxylases

The first part of the catalytic cycle of the pterin-dependent, dioxygen-using nonheme-iron aromatic amino acid hydroxylases, leading to the Fe(IV)=O hydroxylating intermediate, has been investigated by means of density functional theory. The starting structure in the present investigation is the water-free Fe-O(2) complex cluster model that represents the catalytically competent form of the enzymes. A model for this structure was obtained in a previous study of water-ligand dissociation from the hexacoordinate model complex of the X-ray crystal structure of the catalytic domain of phenylalanine hydroxylase in complex with the cofactor (6R)-L-erythro-5,6,7,8-tetrahydrobiopterin (BH(4)) (PAH-Fe(II)-BH(4)). The O-O bond rupture and two-electron oxidation of the cofactor are found to take place via a Fe-O-O-BH(4) bridge structure that is formed in consecutive radical reactions involving a superoxide ion, O(2)(-). The overall effective free-energy barrier to formation of the Fe(IV)=O species is calculated to be 13.9 kcal  mol(-1), less than 2 kcal  mol(-1) lower than that derived from experiment. The rate-limiting step is associated with a one-electron transfer from the cofactor to dioxygen, whereas the spin inversion needed to arrive at the quintet state in which the O-O bond cleavage is finalized, essentially proceeds without activation.

Computational modeling of carbohydrate-recognition process in E-selectin complex: structural mapping of sialyl Lewis X onto ab initio QM/MM free energy surface

To advance our knowledge of carbohydrate recognition by lectins, we propose a systematic computational modeling strategy to identify complex sugar-chain conformations on the reduced free energy surface (FES). We selected the complex of E-selectin with sialyl Lewis X (denoted E-selectin/SLe(x) complex) as a first target molecule. First, we introduced the reduced 2D-FES that characterizes conformational changes in carbohydrate structure as well as the degree of solvation stability of the carbohydrate ligand, and evaluated the overall free energy profile by classical molecular dynamics simulation combined with ab initio QM/MM energy corrections. Second, we mapped flexible carbohydrate structures onto the reduced QM/MM 2D-FES, and identified the details of molecular interactions between each monosaccharide component and the amino acid residues at the carbohydrate-recognition domain. Finally, we confirmed the validity of our modeling strategy by evaluating the chemical shielding tensor by ab initio QM/MM-GIAO computations for several QM/MM-refined geometries sampled from the minimum free energy region in the 2D-FES, and compared this theoretical averaging data with the experimental 1D-NMR profile. The model clearly shows that the binding geometries of the E-selectin/SLe(x) complex are determined not by one single, rigid carbohydrate structure but rather by the sum of averaged conformations fluctuating around the minimum free energy region. For the E-selectin/SLe(x) complex, the major molecular interactions are hydrogen bonds between Fuc and the Ca(2+) binding site in the carbohydrate-recognition domain, and Gal is important in determining the ligand specificity.

Electronic and geometric structures of the blue copper site of azurin investigated by QM/MM hybrid calculations

The electronic and geometric structures of the copper-binding site in a fully solvated azurin were investigated using quantum mechanics (QM) and molecular mechanics (MM) hybrid calculations. Two types of computational models were applied to evaluate the effects of the environment surrounding the active site. In model I, long-distance electrostatic interactions between QM region atoms and partial point charges of the surrounding protein moieties and solvent water were calculated in a QM Hamiltonian, for which the spin-unrestricted Hartree-Fock (UHF)/density functional theory (DFT) hybrid all-electron calculation with the B3LYP functional was adopted. In model II, the QM Hamiltonian was not allowed to be polarized by those partial point charges. Models I and II provided different descriptions of the copper coordination structure, particularly for the coordinative bonds including a large dipole. In fact, the Cu-O(Gly45) and Cu-S(Cys112) bonds are sensitive to the treatment of long-distance electrostatic interactions in the QM Hamiltonian. This suggests that biological processes occurring in the active site are regulated by the surrounding structures of protein and solvent, and therefore the effects of long-range electrostatic interactions involved in the QM Hamiltonian are crucial for accurate descriptions of electronic structures of the copper active site of metalloenzymes.

Progress in ab initio QM/MM free-energy simulations of electrostatic energies in proteins: accelerated QM/MM studies of pKa, redox reactions and solvation free energies

Hybrid quantum mechanical/molecular mechanical (QM/MM) approaches have been used to provide a general scheme for chemical reactions in proteins. However, such approaches still present a major challenge to computational chemists, not only because of the need for very large computer time in order to evaluate the QM energy but also because of the need for proper computational sampling. This review focuses on the sampling issue in QM/MM evaluations of electrostatic energies in proteins. We chose this example since electrostatic energies play a major role in controlling the function of proteins and are key to the structure-function correlation of biological molecules. Thus, the correct treatment of electrostatics is essential for the accurate simulation of biological systems. Although we will be presenting different types of QM/MM calculations of electrostatic energies (and related properties) here, our focus will be on pKa calculations. This reflects the fact that pKa's of ionizable groups in proteins provide one of the most direct benchmarks for the accuracy of electrostatic models of macromolecules. While pKa calculations by semimacroscopic models have given reasonable results in many cases, existing attempts to perform pKa calculations using QM/MM-FEP have led to discrepancies between calculated and experimental values. In this work, we accelerate our QM/MM calculations using an updated mean charge distribution and a classical reference potential. We examine both a surface residue (Asp3) of the bovine pancreatic trypsin inhibitor and a residue buried in a hydrophobic pocket (Lys102) of the T4-lysozyme mutant. We demonstrate that, by using this approach, we are able to reproduce the relevant side chain pKa's with an accuracy of 3 kcal/mol. This is well within the 7 kcal/mol energy difference observed in studies of enzymatic catalysis, and is thus sufficient accuracy to determine the main contributions to the catalytic energies of enzymes. We also provide an overall perspective of the potential of QM/MM calculations in general evaluations of electrostatic free energies, pointing out that our approach should provide a very powerful and accurate tool to predict the electrostatics of not only solution but also enzymatic reactions, as well as the solvation free energies of even larger systems, such as nucleic acid bases incorporated into DNA.