Reaction path potential for complex systems derived from combined ab initio quantum mechanical and molecular mechanical calculations

Combined ab initio quantum mechanical and molecular mechanical calculations have been widely used for modeling chemical reactions in complex systems such as enzymes, with most applications being based on the determination of a minimum energy path connecting the reactant through the transition state to the product in the enzyme environment. However, statistical mechanics sampling and reaction dynamics calculations with a combined ab initio quantum mechanical (QM) and molecular mechanical (MM) potential are still not feasible because of the computational costs associated mainly with the ab initio quantum mechanical calculations for the QM subsystem. To address this issue, a reaction path potential energy surface is developed here for statistical mechanics and dynamics simulation of chemical reactions in enzymes and other complex systems. The reaction path potential follows the ideas from the reaction path Hamiltonian of Miller, Handy and Adams for gas phase chemical reactions but is designed specifically for large systems that are described with combined ab initio quantum mechanical and molecular mechanical methods. The reaction path potential is an analytical energy expression of the combined quantum mechanical and molecular mechanical potential energy along the minimum energy path. An expansion around the minimum energy path is made in both the nuclear and the electronic degrees of freedom for the QM subsystem internal energy, while the energy of the subsystem described with MM remains unchanged from that in the combined quantum mechanical and molecular mechanical expression and the electrostatic interaction between the QM and MM subsystems is described as the interaction of the MM charges with the QM charges. The QM charges are polarizable in response to the changes in both the MM and the QM degrees of freedom through a new response kernel developed in the present work. The input data for constructing the reaction path potential are energies, vibrational frequencies, and electron density response properties of the QM subsystem along the minimum energy path, all of which can be obtained from the combined quantum mechanical and molecular mechanical calculations. Once constructed, it costs much less for its evaluation. Thus, the reaction path potential provides a potential energy surface for rigorous statistical mechanics and reaction dynamics calculations of complex systems. As an example, the method is applied to the statistical mechanical calculations for the potential of mean force of the chemical reaction in triosephosphate isomerase.     

Quantum mechanics/molecular mechanics minimum free-energy path for accurate reaction energetics in solution and enzymes: sequential sampling and optimization on the potential of mean force surface

To accurately determine the reaction path and its energetics for enzymatic and solution-phase reactions, we present a sequential sampling and optimization approach that greatly enhances the efficiency of the ab initio quantum mechanics/molecular mechanics minimum free-energy path (QM/MM-MFEP) method. In the QM/MM-MFEP method, the thermodynamics of a complex reaction system is described by the potential of mean force (PMF) surface of the quantum mechanical (QM) subsystem with a small number of degrees of freedom, somewhat like describing a reaction process in the gas phase. The main computational cost of the QM/MM-MFEP method comes from the statistical sampling of conformations of the molecular mechanical (MM) subsystem required for the calculation of the QM PMF and its gradient. In our new sequential sampling and optimization approach, we aim to reduce the amount of MM sampling while still retaining the accuracy of the results by first carrying out MM phase-space sampling and then optimizing the QM subsystem in the fixed-size ensemble of MM conformations. The resulting QM optimized structures are then used to obtain more accurate sampling of the MM subsystem. This process of sequential MM sampling and QM optimization is iterated until convergence. The use of a fixed-size, finite MM conformational ensemble enables the precise evaluation of the QM potential of mean force and its gradient within the ensemble, thus circumventing the challenges associated with statistical averaging and significantly speeding up the convergence of the optimization process. To further improve the accuracy of the QM/MM-MFEP method, the reaction path potential method developed by Lu and Yang [Z. Lu and W. Yang, J. Chem. Phys. 121, 89 (2004)] is employed to describe the QM/MM electrostatic interactions in an approximate yet accurate way with a computational cost that is comparable to classical MM simulations. The new method was successfully applied to two example reaction processes, the classical SN2 reaction of Cl-+CH3Cl in solution and the second proton transfer step of the reaction catalyzed by the enzyme 4-oxalocrotonate tautomerase. The activation free energies calculated with this new sequential sampling and optimization approach to the QM/MM-MFEP method agree well with results from other simulation approaches such as the umbrella sampling technique with direct QM/MM dynamics sampling, demonstrating the accuracy of the iterative QM/MM-MFEP method.     

Reaction path determination for quantum mechanical/molecular mechanical modeling of enzyme reactions by combining first order and second order "chain-of-replicas" methods

A two-step procedure for the determination of reaction paths in enzyme systems is presented. This procedure combines two chain-of-states methods: a quantum mechanical/molecular mechanical (QM/MM) implementation of the nudged elastic band (NEB) method and a second order parallel path optimizer method both recently developed in our laboratory. In the first step, a reaction path determination is performed with the NEB method, along with a restrained minimization procedure for the MM environment to obtain a first approximation to the reaction path. In the second step, the calculated path is refined with the parallel path optimizer method. By combining these two methods the reaction paths are determined accurately, and in addition, the number of path optimization iterations are significantly reduced. This procedure is tested by calculating both steps of the isomerization of 2-oxo-4-hexenedioate by 4-oxalocrotonate tautomerase, which have been previously determined by our group. The calculated paths agree with the previously reported results and we obtain a reduction of 45%-55% in the number of path optimization cycles.     

Multistructural microiteration technique for geometry optimization and reaction path calculation in large systems

We propose a multistructural microiteration (MSM) method for geometry optimization and reaction path calculation in large systems. MSM is a simple extension of the geometrical microiteration technique. In conventional microiteration, the structure of the non‐reaction‐center (surrounding) part is optimized by fixing atoms in the reaction‐center part before displacements of the reaction‐center atoms. In this method, the surrounding part is described as the weighted sum of multiple surrounding structures that are independently optimized. Then, geometric displacements of the reaction‐center atoms are performed in the mean field generated by the weighted sum of the surrounding parts. MSM was combined with the QM/MM‐ONIOM method and applied to chemical reactions in aqueous solution or enzyme. In all three cases, MSM gave lower reaction energy profiles than the QM/MM‐ONIOM‐microiteration method over the entire reaction paths with comparable computational costs. © 2017 Wiley Periodicals, Inc.     

Enzyme mechanisms with hybrid quantum and molecular mechanical potentials. I. Theoretical considerations

The application of hybrid quantum mechanical and molecular mechanical (QM/MM) potentials to the study of chemical reactions in enzymes is outlined. The discussion is general and addresses the difficulties encountered in an enzyme QM/MM study. First, general criteria for determining whether a particular enzyme is an appropriate candidate for a QM/MM approach are outlined. Methods for obtaining starting structures are detailed. The importance of choosing appropriate levels of ab initio or semiempirical theory is emphasized. Approaches for interfacing the QM and MM regions are briefly discussed, with greater detail given to describing our CHARMM-GAMESS interface. Techniques for partitioning the system into QM and MM regions are explored. Link atom placement, as distant from reacting atoms as possible within the confines of computational efficiency, is examined in some detail. Methods for determining reaction paths are also discussed. © 1996 John Wiley & Sons, Inc.     

Toward a new approach for determination of solute's charge distribution to analyze interatomic electrostatic interactions in quantum mechanical/molecular mechanical simulations

We present an alternative approach to determine "density-dependent property"-derived charges for molecules in the condensed phase. In the case of a solution, it is essential to take into consideration the electron polarization of molecules in the active site of this system. The solute and solvent molecules in this site have to be described by a quantum mechanical technique and the others are allowed to be treated by a molecular mechanical method (QM/MM scheme). For calculations based on this scheme, using the forces and interaction energy as density-dependent property our charges from interaction energy and forces (CHIEF) approach can provide the atom-centered charges on the solute atoms. These charges reproduce well the electrostatic potentials around the solvent molecules and present properly the picture of the electron density of the QM subsystem in the solution system. Thus, the CHIEF charges can be considered as the atomic charges under the conditions of the QM/MM simulation, and then enable one to analyze electrostatic interactions between atoms in the QM and MM regions. This approach would give a view of the QM nuclei and electrons different from the conventional methods.     

Quantum mechanical/molecular mechanical study on the mechanism of the enzymatic Baeyer-Villiger reaction

We report a combined quantum mechanical/molecular mechanical (QM/MM) study on the mechanism of the enzymatic Baeyer-Villiger reaction catalyzed by cyclohexanone monooxygenase (CHMO). In QM/MM geometry optimizations and reaction path calculations, density functional theory (B3LYP/TZVP) is used to describe the QM region consisting of the substrate (cyclohexanone), the isoalloxazine ring of C4a-peroxyflavin, the side chain of Arg-329, and the nicotinamide ring and the adjacent ribose of NADP(+), while the remainder of the enzyme is represented by the CHARMM force field. QM/MM molecular dynamics simulations and free energy calculations at the semiempirical OM3/CHARMM level employ the same QM/MM partitioning. According to the QM/MM calculations, the enzyme-reactant complex contains an anionic deprotonated C4a-peroxyflavin that is stabilized by strong hydrogen bonds with the Arg-329 residue and the NADP(+) cofactor. The CHMO-catalyzed reaction proceeds via a Criegee intermediate having pronounced anionic character. The initial addition reaction has to overcome an energy barrier of about 9 kcal/mol. The formed Criegee intermediate occupies a shallow minimum on the QM/MM potential energy surface and can undergo fragmentation to the lactone product by surmounting a second energy barrier of about 7 kcal/mol. The transition state for the latter migration step is the highest point on the QM/MM energy profile. Gas-phase reoptimizations of the QM region lead to higher barriers and confirm the crucial role of the Arg-329 residue and the NADP(+) cofactor for the catalytic efficiency of CHMO. QM/MM calculations for the CHMO-catalyzed oxidation of 4-methylcyclohexanone reproduce and rationalize the experimentally observed (S)-enantioselectivity for this substrate, which is governed by the conformational preferences of the corresponding Criegee intermediate and the subsequent transition state for the migration step.     

Comparison of linear-scaling semiempirical methods and combined quantum mechanical/molecular mechanical methods for enzymic reactions. II. An energy decomposition analysis

QM/MM methods have been developed as a computationally feasible solution to QM simulation of chemical processes, such as enzyme-catalyzed reactions, within a more approximate MM representation of the condensed-phase environment. However, there has been no independent method for checking the quality of this representation, especially for highly nonisotropic protein environments such as those surrounding enzyme active sites. Hence, the validity of QM/MM methods is largely untested. Here we use the possibility of performing all-QM calculations at the semiempirical PM3 level with a linear-scaling method (MOZYME) to assess the performance of a QM/MM method (PM3/AMBER94 force field). Using two model pathways for the hydride-ion transfer reaction of the enzyme dihydrofolate reductase studied previously (Titmuss et al., Chem Phys Lett 2000, 320, 169-176), we have analyzed the reaction energy contributions (QM, QM/MM, and MM) from the QM/MM results and compared them with analogous-region components calculated via an energy partitioning scheme implemented into MOZYME. This analysis further divided the MOZYME components into Coulomb, resonance and exchange energy terms. For the model in which the MM coordinates are kept fixed during the reaction, we find that the MOZYME and QM/MM total energy profiles agree very well, but that there are significant differences in the energy components. Most significantly there is a large change (approximately 16 kcal/mol) in the MOZYME MM component due to polarization of the MM region surrounding the active site, and which arises mostly from MM atoms close to (<10 A) the active-site QM region, which is not modelled explicitly by our QM/MM method. However, for the model where the MM coordinates are allowed to vary during the reaction, we find large differences in the MOZYME and QM/MM total energy profiles, with a discrepancy of 52 kcal/mol between the relative reaction (product-reactant) energies. This is largely due to a difference in the MM energies of 58 kcal/mol, of which we can attribute approximately 40 kcal/mol to geometry effects in the MM region and the remainder, as before, to MM region polarization. Contrary to the fixed-geometry model, there is no correlation of the MM energy changes with distance from the QM region, nor are they contributed by only a few residues. Overall, the results suggest that merely extending the size of the QM region in the QM/MM calculation is not a universal solution to the MOZYME- and QM/MM-method differences. They also suggest that attaching physical significance to MOZYME Coulomb, resonance and exchange components is problematic. Although we conclude that it would be possible to reparameterize the QM/MM force field to reproduce MOZYME energies, a better way to account for both the effects of the protein environment and known deficiencies in semiempirical methods would be to parameterize the force field based on data from DFT or ab initio QM linear-scaling calculations. Such a force field could be used efficiently in MD simulations to calculate free energies.     

Conical intersections in solution: formulation, algorithm, and implementation with combined quantum mechanics/molecular mechanics method

The significance of conical intersections in photophysics, photochemistry, and photodissociation of polyatomic molecules in gas phase has been demonstrated by numerous experimental and theoretical studies. Optimization of conical intersections of small- and medium-size molecules in gas phase has currently become a routine optimization process, as it has been implemented in many electronic structure packages. However, optimization of conical intersections of small- and medium-size molecules in solution or macromolecules remains inefficient, even poorly defined, due to large number of degrees of freedom and costly evaluations of gradient difference and nonadiabatic coupling vectors. In this work, based on the sequential quantum mechanics and molecular mechanics (QM/MM) and QM/MM-minimum free energy path methods, we have designed two conical intersection optimization methods for small- and medium-size molecules in solution or macromolecules. The first one is sequential QM conical intersection optimization and MM minimization for potential energy surfaces; the second one is sequential QM conical intersection optimization and MM sampling for potential of mean force surfaces, i.e., free energy surfaces. In such methods, the region where electronic structures change remarkably is placed into the QM subsystem, while the rest of the system is placed into the MM subsystem; thus, dimensionalities of gradient difference and nonadiabatic coupling vectors are decreased due to the relatively small QM subsystem. Furthermore, in comparison with the concurrent optimization scheme, sequential QM conical intersection optimization and MM minimization or sampling reduce the number of evaluations of gradient difference and nonadiabatic coupling vectors because these vectors need to be calculated only when the QM subsystem moves, independent of the MM minimization or sampling. Taken together, costly evaluations of gradient difference and nonadiabatic coupling vectors in solution or macromolecules can be reduced significantly. Test optimizations of conical intersections of cyclopropanone and acetaldehyde in aqueous solution have been carried out successfully.     

Contributions of pauli repulsions to the energetics and physical properties computed in QM/MM methods

Conventional combined quantum mechanical/molecular mechanical (QM/MM) methods lack explicit treatment of Pauli repulsions between the quantum‐mechanical and molecular‐mechanical subsystems. Instead, classical Lennard‐Jones (LJ) potentials between QM and MM nuclei are used to model electronic Pauli repulsion and long‐range London dispersion, despite the fact that the latter two are inherently of quantum nature. Use of the simple LJ potential in QM/MM methods can reproduce minimal geometries and energies of many molecular clusters reasonably well, as compared to full QM calculations. However, we show here that the LJ potential cannot correctly describe subtle details of the electron density of the QM subsystem because of the neglect of Pauli repulsions between the QM and MM subsystems. The inaccurate electron density subsequently affects the calculation of electronic and magnetic properties of the QM subsystem. To explicitly consider Pauli interactions with QM/MM methods, we propose a method to use empirical effective potentials on the MM atoms. The test case of the binding energy and magnetic properties of a water dimer shows promising results for the general application of effective potentials to mimic Pauli repulsions in QM/MM calculations. © 2013 Wiley Periodicals, Inc.     

Variational calculation of quantum mechanical/molecular mechanical free energy with electronic polarization of solvent

Quantum mechanical/molecular mechanical (QM/MM) free energy calculation presents a significant challenge due to an excessive number of QM calculations. A useful approach for reducing the computational cost is that based on the mean field approximation to the QM subsystem. Here, we describe such a mean-field QM/MM theory for electronically polarizable systems by starting from the Hartree product ansatz for the total system and invoking a variational principle of free energy. The MM part is then recast to a classical polarizable model by introducing the charge response kernel. Numerical test shows that the potential of mean force (PMF) thus obtained agrees quantitatively with that obtained from a direct QM/MM calculation, indicating the utility of self-consistent mean-field approximation. Next, we apply the obtained method to prototypical reactions in several qualitatively different solvents and make a systematic comparison of polarization effects. The results show that in aqueous solution the PMF does not depend very much on the water models employed, while in nonaqueous solutions the PMF is significantly affected by explicit polarization. For example, the free energy barrier for a phosphoryl dissociation reaction in acetone and cyclohexane is found to increase by more than 10 kcal/mol when switching the solvent model from an empirical to explicitly polarizable one. The reason for this is discussed based on the parametrization of empirical nonpolarizable models.     

Modeling of serine protease prototype reactions with the flexible effective fragment potential quantum mechanical/molecular mechanical method

A complete cycle of chemical transformations for the serine protease prototype reaction is modeled following calculations with the flexible effective fragment quantum mechanical/molecular mechanical (QM/MM) method. The initial molecular model is based on the crystal structure of the trypsin–bovine pancreatic trypsin inhibitor complex including all atoms of the enzyme within approximately 15–18 Å of the oxygen center O of the catalytic serine residue. Several selections of the QM/MM partitioning are considered. Fractions of the side chains of the residues from the catalytic triad (serine, histidine and aspartic acid) and a central part of a model substrate around the C–N bond to be cleaved are included into the QM subsystem. The remaining part, or the MM subsystem, is represented by flexible chains of small effective fragments, whose potentials explicitly contribute to the Hamiltonian of the QM part, but the corresponding fragment–fragment interactions are described by the MM force fields. The QM/MM boundaries are extended over the C–C bonds of the peptides assigned to the QM subsystem in the enzyme, C–C and C–N bonds in model substrates. Multiple geometry optimizations have been performed by using the RHF/6-31G method in the QM part and OPLSAA or AMBER sets of MM parameters, resulting in a series of stationary points on the complex potential-energy surfaces. All structures generally accepted for the serine protease catalytic cycle have been located. Energies at the stationary points found have been recomputed at the MP2/6-31+G^* level for the QM part in the protein environment. Structural changes along the reaction path are analyzed with special attention to hydrogen-bonding networks. In the case of a model substrate selected as a short peptide CH₃(NHCO-CH₂)₂ – HN–CO–(CH₂–NHCO)CH₃ the computed energy profile for the acylation step shows too high activation energy barriers. The energetics of this rate-limiting step is considerably improved, if more realistic model for the substrate is considered, following the motifs of the ThrI11–GlyI12–ProI13-–CysI14–LysI15–AlaI16–ArgI17–IleI18–IleI19 sequence of the bovine pancreatic trypsin inhibitor.     

Exploring potential energy surfaces for chemical reactions: an overview of some practical methods

Potential energy surfaces form a central concept in the application of electronic structure methods to the study of molecular structures, properties, and reactivities. Recent advances in tools for exploring potential energy surfaces are surveyed. Methods for geometry optimization of equilibrium structures, searching for transition states, following reaction paths and ab initio molecular dynamics are discussed. For geometry optimization, topics include methods for large molecules, QM/MM calculations, and simultaneous optimization of the wave function and the geometry. Path optimization methods and dynamics based techniques for transition state searching and reaction path following are outlined. Developments in the calculation of ab initio classical trajectories in the Born-Oppenheimer and Car-Parrinello approaches are described.     

Catalytic mechanism of glycosyltransferases: hybrid quantum mechanical/molecular mechanical study of the inverting N-acetylglucosaminyltransferase I

The Golgi glycosyltransferase, N-acetylglucosaminyltransferase I (GnT-I), catalyzes the transfer of a GlcNAc residue from the donor UDP-GlcNAc to the C2-hydroxyl group of a mannose residue in the trimannosyl core of the Man5GlcNAc2-Asn-X oligosaccharide. The catalytic mechanism of GnT-I was investigated using a hybrid quantum mechanical/molecular mechanical (QM/MM) method with a QM part containing 88 atoms treated with density functional theory (DFT) at the BP/TZP level. The remaining parts of a GnT-I complex, altogether 5633 atoms, were modeled using the AMBER molecular force field. A theoretical model of a Michaelis complex was built using the X-ray structure of GnT-I in complex with the donor having geometrical features consistent with kinetic studies. The QM(DFT)/MM model identified a concerted SN2-type of transition state with D291 as the catalytic base for the reaction in the enzyme active site. The TS model features nearly simultaneous nucleophilic addition and dissociation steps accompanied by the transfer of the nucleophile proton Hb2 to the catalytic base D291. The structure of the TS model is characterized by the Ob2-C1 and C1-O1 bond distances of 1.912 and 2.542 A, respectively. The activation energy for the proposed reaction mechanism was estimated to be approximately 19 kcal mol-1. The calculated alpha-deuterium kinetic isotope effect of 1.060 is consistent with the proposed reaction mechanism. Theoretical results also identified interactions between the Hb6 and beta-phosphate oxygen of the UDP and a low-barrier hydrogen bond between the nucleophile and the catalytic base D291. It is proposed that these interactions contribute to a stabilization of TS. This modeling study provided detailed insight into the mechanism of the GlcNAc transfer catalyzed by GnT-I, which is the first step in the conversion of high mannose oligosaccharides to complex and hybrid N-glycan structures.     

Importance of accurate charges in molecular docking: quantum mechanical/molecular mechanical (QM/MM) approach

The extent to which accuracy of electric charges plays a role in protein-ligand docking is investigated through development of a docking algorithm, which incorporates quantum mechanical/molecular mechanical (QM/MM) calculations. In this algorithm, fixed charges of ligands obtained from force field parameterization are replaced by QM/MM calculations in the protein environment, treating only the ligands as the quantum region. The algorithm is tested on a set of 40 cocrystallized structures taken from the Protein Data Bank (PDB) and provides strong evidence that use of nonfixed charges is important. An algorithm, dubbed "Survival of the Fittest" (SOF) algorithm, is implemented to incorporate QM/MM charge calculations without any prior knowledge of native structures of the complexes. Using an iterative protocol, this algorithm is able in many cases to converge to a nativelike structure in systems where redocking of the ligand using a standard fixed charge force field exhibits nontrivial errors. The results demonstrate that polarization effects can play a significant role in determining the structures of protein-ligand complexes, and provide a promising start towards the development of more accurate docking methods for lead optimization applications.     

The dynamo library for molecular simulations using hybrid quantum mechanical and molecular mechanical potentials

The Dynamo module library has been developed for the simulation of molecular systems using hybrid quantum mechanical (QM) and molecular mechanical (MM) potentials. Dynamo is not a program package but is a library of Fortran 90 modules that can be employed by those interested in writing their own programs for performing molecular simulations. The library supports a range of different types of molecular calculation including geometry optimizations, reaction‐path determinations and molecular dynamics and Monte Carlo simulations. This article outlines the general structure and capabilities of the library and describes in detail Dynamo's semiempirical QM/MM hybrid potential. Results are presented to indicate three particular aspects of this implementation—the handling of long‐range nonbonding interactions, the nature of the boundary between the quantum mechanical and molecular mechanical atoms and how to perform path‐integral hybrid‐potential molecular dynamics simulations. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 1088–1100, 2000     

YinYang atom: a simple combined ab initio quantum mechanical molecular mechanical model

A simple interface is proposed for combined quantum mechanical (QM) molecular mechanical (MM) calculations for the systems where the QM and MM regions are connected through covalent bonds. Within this model, the atom that connects the two regions, called YinYang atom here, serves as an ordinary MM atom to other MM atoms and as a hydrogen-like atom to other QM atoms. Only one new empirical parameter is introduced to adjust the length of the connecting bond and is calibrated with the molecule propanol. This model is tested with the computation of equilibrium geometries and protonation energies for dozens of molecules. Special attention is paid on the influence of MM point charges on optimized geometry and protonation energy, and it is found that it is important to maintain local charge-neutrality in the MM region in order for the accurate calculation of the protonation and deprotonation energies. Overall the simple YinYang atom model yields comparable results to some other QM/MM models.