A novel quantum mechanical/molecular mechanical approach to the free energy calculation for isomerization of glycine in aqueous solution

The free energy change associated with the isomerization reaction of glycine in water solution has been studied by a hybrid quantum mechanical/molecular mechanical (QM/MM) approach combined with the theory of energy representation (QM/MM-ER) recently developed. The solvation free energies for both neutral and zwitterionic form of glycine have been determined by means of the QM/MM-ER simulation. The contributions of the electronic polarization and the fluctuation of the QM solute to the solvation free energy have been investigated. It has been found that the contribution of the density fluctuation of the zwitterionic solute is estimated as -4.2 kcal/mol in the total solvation free energy of -46.1 kcal/mol, while that of the neutral form is computed as -3.0 kcal/mol in the solvation free energy of -15.6 kcal/mol. The resultant free energy change associated with the isomerization of glycine in water has been obtained as -7.8 kcal/mol, in excellent agreement with the experimental data of -7.3 or -7.7 kcal/mol, implying the accuracy of the QM/MM-ER approach. The results have also been compared with those computed by other methodologies such as the polarizable continuum model and the classical molecular simulation. The efficiency and advantage of the QM/MM-ER method has been discussed.     

Converging ligand‐binding free energies obtained with free‐energy perturbations at the quantum mechanical level

In this article, the convergence of quantum mechanical (QM) free‐energy simulations based on molecular dynamics simulations at the molecular mechanics (MM) level has been investigated. We have estimated relative free energies for the binding of nine cyclic carboxylate ligands to the octa‐acid deep‐cavity host, including the host, the ligand, and all water molecules within 4.5 Å of the ligand in the QM calculations (158–224 atoms). We use single‐step exponential averaging (ssEA) and the non‐Boltzmann Bennett acceptance ratio (NBB) methods to estimate QM/MM free energy with the semi‐empirical PM6‐DH2X method, both based on interaction energies. We show that ssEA with cumulant expansion gives a better convergence and uses half as many QM calculations as NBB, although the two methods give consistent results. With 720,000 QM calculations per transformation, QM/MM free‐energy estimates with a precision of 1 kJ/mol can be obtained for all eight relative energies with ssEA, showing that this approach can be used to calculate converged QM/MM binding free energies for realistic systems and large QM partitions. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Investigation of the dominant hydration structures among the ionic species in aqueous solution: novel quantum mechanics/molecular mechanics simulations combined with the theory of energy representation

In the present work, we have performed quantum chemical calculations to determine preferable species among the ionic complexes that are present in ambient water due to the autodissociation of water molecule. First, we have formulated the relative population of the hydrated complexes with respect to the bare ion (H(3)O(+) or OH(-)) in terms of the solvation free energies of the relevant molecules. The solvation free energies for various ionic species (H(3)O(+), H(5)O(2) (+), H(7)O(3) (+), H(9)O(4) (+) or OH(-), H(3)O(2) (-), H(5)O(3) (-), H(7)O(4) (-), H(9)O(5) (-)), categorized as proton or hydroxide ion in solution, have been computed by employing the QM/MM-ER method recently developed by combining the quantum mechanical/molecular mechanical (QM/MM) approach with the theory of energy representation (ER). Then, the computed solvation free energies have been used to evaluate the ratio of the populations of the ionic complexes to that of the bare ion (H(3)O(+) or OH(-)). Our results suggest that the Zundel form, i.e., H(5)O(2) (+), is the most preferable in the solution among the cationic species listed above though the Eigen form (H(9)O(4) (+)) is very close to the Zundel complex in the free energy, while the anionic fragment from water molecules mostly takes the form of OH(-). It has also been found that the loss of the translational entropy of water molecules associated with the formation of the complex plays a role in determining the preferable size of the cluster.     

QM/MM-PBSA method to estimate free energies for reactions in proteins

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.     

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.     

Influence of inter- and intramolecular hydrogen bonding on kemp decarboxylations from QM/MM simulations

The Kemp decarboxylation reaction for benzisoxazole-3-carboxylic acid derivatives has been investigated using QM/MM calculations in protic and dipolar aprotic solvents. Aprotic solvents have been shown to accelerate the rates of reaction by 7-8 orders of magnitude over water; however, the inclusion of an internal hydrogen bond effectively inhibits the reaction with near solvent independence. The effects of solvation and intramolecular hydrogen bonding on the reactants, transition structures, and the rate of reaction are elucidated using two-dimensional potentials of mean force (PMF) derived from free energy perturbation calculations in Monte Carlo simulations (MC/FEP). Free energies of activation in six solvents have been computed to be in close agreement with experiment. Solute-solvent interaction energies show that poorer solvation of the reactant anion in the dipolar aprotic solvents is primarily responsible for the observed rate enhancements over protic media. In addition, a discrepancy for the experimental rate in chloroform has been studied in detail with the conclusion that ion-pairing between the reactant anion and tetramethylguanidinium counterion is responsible for the anomalously slow reaction rate. The overall quantitative success of the computations supports the present QM/MM/MC approach, which features PDDG/PM3 as the QM method.     

Performance assessment of semiempirical molecular orbital methods in describing halogen bonding: quantum mechanical and quantum mechanical/molecular mechanical-molecular dynamics study

The performance of semiempirical molecular-orbital methods--MNDO, MNDO-d, AM1, RM1, PM3 and PM6--in describing halogen bonding was evaluated, and the results were compared with molecular mechanical (MM) and quantum mechanical (QM) data. Three types of performance were assessed: (1) geometrical optimizations and binding energy calculations for 27 halogen-containing molecules complexed with various Lewis bases (Two of the tested methods, AM1 and RM1, gave results that agree with the QM data.); (2) charge distribution calculations for halobenzene molecules, determined by calculating the solvation free energies of the molecules relative to benzene in explicit and implicit generalized Born (GB) solvents (None of the methods gave results that agree with the experimental data.); and (3) appropriateness of the semiempirical methods in the hybrid quantum-mechanical/molecular-mechanical (QM/MM) scheme, investigated by studying the molecular inhibition of CK2 protein by eight halobenzimidazole and -benzotriazole derivatives using hybrid QM/MM molecular-dynamics (MD) simulations with the inhibitor described at the QM level by the AM1 method and the rest of the system described at the MM level. The pure MM approach with inclusion of an extra point of positive charge on the halogen atom approach gave better results than the hybrid QM/MM approach involving the AM1 method. Also, in comparison with the pure MM-GBSA (generalized Born surface area) binding energies and experimental data, the calculated QM/MM-GBSA binding energies of the inhibitors were improved by replacing the G(GB,QM/MM) solvation term with the corresponding G(GB,MM) term.     

Coupled semiempirical quantum mechanics and molecular mechanics (QM/MM) calculations on the aqueous solvation free energies of ionized molecules

The aqueous solvation free energies of ionized molecules were computed using a coupled quantum mechanical and molecular mechanical (QM/MM) model based on the AM1, MNDO, and PM3 semiempirical molecular orbital methods for the solute molecule and the TIP3P molecular mechanics model for liquid water. The present work is an extension of our model for neutral solutes where we assumed that the total free energy is the sum of components derived from the electrostatic/polarization terms in the Hamiltonian plus an empirical “nonpolar” term. The electrostatic/polarization contributions to the solvation free energies were computed using molecular dynamics (MD) simulation and thermodynamic integration techniques, while the nonpolar contributions were taken from the literature. The contribution to the electrostatic/polarization component of the free energy due to nonbonded interactions outside the cutoff radii used in the MD simulations was approximated by a Born solvation term. The experimental free energies were reproduced satisfactorily using variational parameters from the vdW terms as in the original model, in addition to a parameter from the one-electron integral terms. The new one-electron parameter was required to account for the short-range effects of overlapping atomic charge densities. The radial distribution functions obtained from the MD simulations showed the expected H-bonded structures between the ionized solute molecule and solvent molecules. We also obtained satisfactory results by neglecting both the empirical nonpolar term and the electronic polarization of the solute, i.e., by implementing a nonpolarization model. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1028–1038, 1999     

Simple and exact approach to the electronic polarization effect on the solvation free energy: formulation for quantum-mechanical/molecular-mechanical system and its applications to aqueous solutions

We present a simple and exact numerical approach to compute the free energy contribution δμ in solvation due to the electron density polarization and fluctuation of a quantum-mechanical solute in the quantum-mechanical/molecular-mechanical (QM/MM) simulation combined with the theory of the energy representation (QM/MM-ER). Since the electron density fluctuation is responsible for the many-body QM-MM interactions, the standard version of the energy representation method cannot be applied directly. Instead of decomposing the QM-MM polarization energy into the pairwise additive and non-additive contributions, we take sum of the polarization energies in the QM-MM interaction and adopt it as a new energy coordinate for the method of energy representation. Then, it is demonstrated that the free energy δμ can be exactly formulated in terms of the energy distribution functions for the solution and reference systems with respect to this energy coordinate. The benchmark tests were performed to examine the numerical efficiency of the method with respect to the changes in the individual properties of the solvent and the solute. Explicitly, we computed the solvation free energy of a QM water molecule in ambient and supercritical water, and also the free-energy change associated with the isomerization reaction of glycine from neutral to zwitterionic structure in aqueous solution. In all the systems examined, it was demonstrated that the computed free energy δμ agrees with the experimental value, irrespective of the choice of the reference electron density of the QM solute. The present method was also applied to a prototype reaction of adenosine 5'-triphosphate hydrolysis where the effect of the electron density fluctuation is substantial due to the excess charge. It was demonstrated that the experimental free energy of the reaction has been accurately reproduced with the present approach.     

Conformational dependence of charges in protein simulations

We have studied the conformational dependence of molecular mechanics atomic charges for proteins by calculating the charges fitted to the quantum mechanical (QM) electrostatic potential (ESP) for all atoms in complexes between avidin and seven biotin analogues for 20 snapshots from molecular dynamics simulations. We have studied how various other charge sets reproduce those charges. The QM charges, even if averaged over all snapshots or all residues, in general have a larger magnitude than standard Amber charges, indicating that the restraint toward zero in the restrained ESP method is too strong. This has a significant influence on the electrostatic conformational energies and the interaction energy between the biotin ligand and the protein, giving a difference between the QM and Amber charges of 43 and 8 kJ/mol for the negatively charged and neutral biotin analogues, respectively (3-4%). However, this energy difference is strongly reduced if the solvation energy (calculated by the Poisson-Boltzmann or Generalized Born methods) is added, viz., to 7 kJ/mol for charged and 3 kJ/mol for uncharged ligand. In fact, charges need to be recalculated with a QM method only for residues within 7 or 4 A of the ligand, if the error should be less than 4 kJ/mol. Unfortunately, the QM charges do not give significantly better MM/PBSA estimates of ligand-binding affinities than standard Amber charges.     

Solvation free energies of molecules. The most stable anionic tautomers of uracil

Anionic states of nucleic acid bases are suspected to play a role in the radiation damage processes of DNA. Our recent studies suggested that the excess electron attachment to the nucleic acid bases can stabilize some rare tautomers, i.e. imine-enamine tautomers and other tautomers with a proton being transferred from nitrogen sites to carbon sites (with respect to the canonical tautomer). So far, these new anionic tautomers have been characterized by the gas-phase electronic structure calculations and photoelectron spectroscopy experiments. In the current contribution we explore the effect of water solvation on the stability of the new anionic tautomers of uracil. The accurate free energies of solvation are calculated in a two step approach. The major contribution was calculated using the classical free-energy perturbation adiabatic-charging approach, where it is assumed that the solvated molecule has the charge distribution given by the polarizable continuum model. In the second step the free energy of solvation is refined by taking into account the real, average solvent charge distribution. This is done using our accelerated QM/MM simulations, where the QM energy of the solute is calculated in the mean potential averaged over many MD steps. We found that in water solution three of the recently identified anionic tautomers are 6.5-3.6 kcal mol(-1) more stable than the anion of the canonical tautomer.     

Aqueous solvation free energies of ions and ion-water clusters based on an accurate value for the absolute aqueous solvation free energy of the proton

Thermochemical cycles that involve pKa, gas-phase acidities, aqueous solvation free energies of neutral species, and gas-phase clustering free energies have been used with the cluster pair approximation to determine the absolute aqueous solvation free energy of the proton. The best value obtained in this work is in good agreement with the value reported by Tissandier et al. (Tissandier, M. D.; Cowen, K. A.; Feng, W. Y.; Gundlach, E.; Cohen, M. J.; Earhart, A. D.; Coe, J. V. J. Phys. Chem. A 1998, 102, 7787), who applied the cluster pair approximation to a less diverse and smaller data set of ions. We agree with previous workers who advocated the value of -265.9 kcal/mol for the absolute aqueous solvation free energy of the proton. Considering the uncertainties associated with the experimental gas-phase free energies of ions that are required to use the cluster pair approximation as well as analyses of various subsets of data, we estimate an uncertainty for the absolute aqueous solvation free energy of the proton of no less than 2 kcal/mol. Using a value of -265.9 kcal/mol for the absolute aqueous solvation free energy of the proton, we expand and update our previous compilation of absolute aqueous solvation free energies; this new data set contains conventional and absolute aqueous solvation free energies for 121 unclustered ions (not including the proton) and 147 conventional and absolute aqueous solvation free energies for 51 clustered ions containing from 1 to 6 water molecules. When tested against the same set of ions that was recently used to develop the SM6 continuum solvation model, SM6 retains its previously determined high accuracy; indeed, in most cases the mean unsigned error improves when it is tested against the more accurate reference data.     

QM/MM study of aqueous solvation of the uranyl fluoride [UO2F4(2-)] complex

The aqueous solvation of the uranylfluoride complex [UO(2)F(4) (2-)] was studied using full quantum mechanical (QM) and hybrid QM/molecular mechanics (MM) methods. Inclusion of a complete first solvation shell was found necessary to reproduce the experimentally observed heptacoordination of uranium. An efficient and accurate computational model is proposed that consists of structure optimization of the coordinated uranium complex as QM region, followed by single-point full QM calculations to compute relative energies. This method is proven feasible for studies of large solvated actinide complexes.     

Binding free energies in the SAMPL6 octa-acid host–guest challenge calculated with MM and QM methods

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²) 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²?=?0.1–0.9), but the results were improved compared to previous studies of this system with similar methods.     

Fragment‐based quantum mechanical calculation of protein–protein binding affinities

The electrostatically embedded generalized molecular fractionation with conjugate caps (EE‐GMFCC) method has been successfully utilized for efficient linear‐scaling quantum mechanical (QM) calculation of protein energies. In this work, we applied the EE‐GMFCC method for calculation of binding affinity of Endonuclease colicin–immunity protein complex. The binding free energy changes between the wild‐type and mutants of the complex calculated by EE‐GMFCC are in good agreement with experimental results. The correlation coefficient (R) between the predicted binding energy changes and experimental values is 0.906 at the B3LYP/6‐31G*‐D level, based on the snapshot whose binding affinity is closest to the average result from the molecular mechanics/Poisson–Boltzmann surface area (MM/PBSA) calculation. The inclusion of the QM effects is important for accurate prediction of protein–protein binding affinities. Moreover, the self‐consistent calculation of PB solvation energy is required for accurate calculations of protein–protein binding free energies. This study demonstrates that the EE‐GMFCC method is capable of providing reliable prediction of relative binding affinities for protein–protein complexes. © 2018 Wiley Periodicals, Inc.     

Computational methods for the study of enzymic reaction mechanisms III: a perturbation plus QM/MM approach for calculating relative free energies of protonation

We describe a coupling parameter, that is, perturbation, approach to effectively create and annihilate atoms in the quantum mechanical Hamiltonian within the closed shell restricted Hartree-Fock formalism. This perturbed quantum mechanical atom (PQA) method is combined with molecular mechanics (MM) methods (PQA/MM) within a molecular dynamics simulation, to model the protein environment (MM region) effects that also make a contribution to the overall free energy change. Using the semiempirical PM3 method to model the QM region, the application of this PQA/MM method is illustrated by calculation of the relative protonation free energy of the conserved OD2 (Asp27) and the N5 (dihydrofolate) proton acceptor sites in the active site of Escherichia coli dihydrofolate reductase (DHFR) with the bound nicotinamide adenine dinucleotide phosphate (NADPH) cofactor. For a number of choices for the QM region, the relative protonation free energy was calculated as the sum of contributions from the QM region and the interaction between the QM and MM regions via the thermodynamic integration (TI) method. The results demonstrate the importance of including the whole substrate molecule in the QM region, and the overall protein (MM) environment in determining the relative stabilities of protonation sites in the enzyme active site. The PQA/MM free energies obtained by TI were also compared with those estimated by a less computationally demanding nonperturbative method based on the linear response approximation (LRA). For some choices of QM region, the total free energies calculated using the LRA method were in very close agreement with the PQA/MM values. However, the QM and QM/MM component free energies were found to differ significantly between the two methods.     

A QM/MM study of the binding of RAPTA ligands to cathepsin B

We have carried out quantum mechanical (QM) and QM/MM (combined QM and molecular mechanics) calculations, as well as molecular dynamics (MD) simulations to study the binding of a series of six RAPTA (Ru(II)-arene-1,3,5-triaza-7-phosphatricyclo-[3.3.1.1] decane) complexes with different arene substituents to cathepsin B. The recently developed QM/MM-PBSA approach (QM/MM combined with Poisson–Boltzmann solvent-accessible surface area solvation) has been used to estimate binding affinities. The QM calculations reproduce the antitumour activities of the complexes with a correlation coefficient (r ²) of 0.35–0.86 after a conformational search. The QM/MM-PBSA method gave a better correlation (r ² = 0.59) when the protein was fixed to the crystal structure, but more reasonable ligand structures and absolute binding energies were obtained if the protein was allowed to relax, indicating that the ligands are strained when the protein is kept fixed. In addition, the best correlation (r ² = 0.80) was obtained when only the QM energies were used, which suggests that the MM and continuum solvation energies are not accurate enough to predict the binding of a charged metal complex to a charged protein. Taking into account the protein flexibility by means of MD simulations slightly improves the correlation (r ² = 0.91), but the absolute energies are still too large and the results are sensitive to the details in the calculations, illustrating that it is hard to obtain stable predictions when full flexible protein is included in the calculations.     

Use of a QM/MM-based FEP method to evaluate the anomalous hydration behavior of simple alkyl amines and amides: application to the design of FBPase inhibitors for the treatment of type-2 diabetes

Standard molecular mechanics (MM) force fields predict a nearly linear decrease in hydration free energy with each successive addition of a methyl group to ammonia or acetamide, whereas a nonadditive relationship is observed experimentally. In contrast, the non-additive hydration behavior is reproduced directly using a quantum mechanics (QM)/MM-based free-energy perturbation (FEP) method wherein the solute partial atomic charges are updated at every window. Decomposing the free energies into electrostatic and van der Waals contributions and comparing the results with the corresponding free energies obtained using a conventional FEP method and a QM/MM method wherein the charges are not updated suggests that inaccuracies in the electrostatic free energies are the primary reason for the inability of the conventional FEP method to predict the experimental findings. The QM/MM-based FEP method was subsequently used to evaluate inhibitors of the diabetes drug target fructose-1,6-bisphosphatase adenosine 5'-monophosphate and 6-methylamino purine riboside 5'-monophosphate. The predicted relative binding free energy was consistent with the experimental findings, whereas the relative binding free energy predicted using the conventional FEP method differed from the experimental finding by an amount consistent with the overestimated relative solvation free energies calculated for alkylamines. Accordingly, the QM/MM-based FEP method offers potential advantages over conventional FEP methods, including greater accuracy and reduced user input. Moreover, since drug candidates often contain either functionality that is inadequately treated by MM (e.g., simple alkylamines and alkylamides) or new molecular scaffolds that require time-consuming development of MM parameters, these advantages could enable future automation of FEP calculations as well as greatly increase the use and impact of FEP calculations in drug discovery.