Chemical accuracy in QM/MM calculations on enzyme-catalysed reactions

Combined quantum mechanics/molecular mechanics (QM/MM) modelling has the potential to answer fundamental questions about enzyme mechanisms and catalysis. Calculations using QM/MM methods can now predict barriers for enzyme-catalysed reactions with unprecedented, near chemical accuracy, i.e. to within 1 kcal/mol in the best cases. Quantitative predictions from first-principles calculations were only previously possible for very small molecules. At this level, quantitative, reliable predictions can be made about the mechanisms of enzyme-catalysed reactions. This development signals a new era of computational biochemistry.     

The adaptive buffered force QM/MM method in the CP2K and AMBER software packages

The implementation and validation of the adaptive buffered force (AdBF) quantum‐mechanics/molecular‐mechanics (QM/MM) method in two popular packages, CP2K and AMBER are presented. The implementations build on the existing QM/MM functionality in each code, extending it to allow for redefinition of the QM and MM regions during the simulation and reducing QM‐MM interface errors by discarding forces near the boundary according to the buffered force‐mixing approach. New adaptive thermostats, needed by force‐mixing methods, are also implemented. Different variants of the method are benchmarked by simulating the structure of bulk water, water autoprotolysis in the presence of zinc and dimethyl‐phosphate hydrolysis using various semiempirical Hamiltonians and density functional theory as the QM model. It is shown that with suitable parameters, based on force convergence tests, the AdBF QM/MM scheme can provide an accurate approximation of the structure in the dynamical QM region matching the corresponding fully QM simulations, as well as reproducing the correct energetics in all cases. Adaptive unbuffered force‐mixing and adaptive conventional QM/MM methods also provide reasonable results for some systems, but are more likely to suffer from instabilities and inaccuracies. © 2015 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

QM/MM vibrational mode tracking

Vibrational spectroscopy is a powerful tool to investigate the structure and dynamics of biomolecules. When small subsystems of large molecules such as active centers of enzymes are studied, quantum chemical calculations based on quantum mechanics/molecular mechanics (QM/MM) coupling schemes are a valuable means to interpret the spectra. The goal of this work is a methodological pilot study on how to selectively and thus efficiently extract certain vibrational information for extended molecular systems described by QM/MM methods. This is achieved by an extension of the mode tracking algorithm and a comparison with the partial Hessian diagonalization approach. After validating the methodology for the CO stretching vibration of 2-butanone and a delocalized CO stretch in acetylacetone, the stretching and bending modes of the CO ligand in CO myoglobin are tracked. Such systems represent an ideal application for mode tracking, because only a few strongly localized vibrations are sought for, while the large remainder of the molecule is of interest only as far as it affects these local vibrations. This influence is treated exactly by mode tracking.     

Reliable treatment of electrostatics in combined QM/MM simulation of macromolecules

A robust approach for dealing with electrostatic interactions for spherical boundary conditions has been implemented in the QM/MM framework. The development was based on the generalized solvent boundary potential (GSBP) method proposed by Im et al. [J. Chem. Phys. 114, 2924 (2001)], and the specific implementation was applied to the self-consistent-charge density-functional tight-binding approach as the quantum mechanics (QM) level, although extension to other QM methods is straightforward. Compared to the popular stochastic boundary-condition scheme, the new protocol offers a balanced treatment between quantum mechanics/molecular mechanics (QM/MM) and MM/MM interactions; it also includes the effect of the bulk solvent and macromolecule atoms outside of the microscopic region at the Poisson-Boltzmann level. The new method was illustrated with application to the enzyme human carbonic anhydrase II and compared to stochastic boundary-condition simulations using different electrostatic treatments. The GSBP-based QM/MM simulations were most consistent with available experimental data, while conventional stochastic boundary simulations yielded various artifacts depending on different electrostatic models. The results highlight the importance of carefully treating electrostatics in QM/MM simulations of biomolecules and suggest that the commonly used truncation schemes should be avoided in QM/MM simulations, especially in simulations that involve extensive conformational samplings. The development of the GSBP-based QM/MM protocol has opened up the exciting possibility of studying chemical events in very complex biomolecular systems in a multiscale framework.     

Geometry optimization with QM/MM,ONIOM, and other combined methods. I. Microiterations and constraints

Hybrid energy methods such as QM/MM and ONIOM, that combine different levels of theory into one calculation, have been very successful in describing large systems. Geometry optimization methods can take advantage of the partitioning of these calculations into a region treated at a quantum mechanical (QM) level of theory and the larger, remaining region treated by an inexpensive method such as molecular mechanics (MM). A series of microiterations can be employed to fully optimize the MM region for each optimization step in the QM region. Cartesian coordinates are used for the MM region and are chosen so that the internal coordinates of the QM region remain constant during the microiterations. The coordinates of the MM region are augmented to permit rigid body translation and rotation of the QM region. This is essential if any atoms in the MM region are constrained, but it also improves the efficiency of unconstrained optimizations. Because of the microiterations, special care is needed for the optimization step in the QM region so that the system remains in the same local valley during the course of the optimization. The optimization methodology with microiterations, constraints, and step-size control are illustrated by calculations on bacteriorhodopsin and other systems.     

Fragmentation-based QM/MM simulations: length dependence of chain dynamics and hydrogen bonding of polyethylene oxide and polyethylene in aqueous solutions

Molecular fragmentation quantum mechanics (QM) calculations have been combined with molecular mechanics (MM) to construct the fragmentation QM/MM method for simulations of dilute solutions of macromolecules. We adopt the electrostatics embedding QM/MM model, where the low-cost generalized energy-based fragmentation calculations are employed for the QM part. Conformation energy calculations, geometry optimizations, and Born-Oppenheimer molecular dynamics simulations of poly(ethylene oxide), PEO(n) (n = 6-20), and polyethylene, PE(n) ( n = 9-30), in aqueous solution have been performed within the framework of both fragmentation and conventional QM/MM methods. The intermolecular hydrogen bonding and chain configurations obtained from the fragmentation QM/MM simulations are consistent with the conventional QM/MM method. The length dependence of chain conformations and dynamics of PEO and PE oligomers in aqueous solutions is also investigated through the fragmentation QM/MM molecular dynamics simulations.     

Calculation of wave‐functions with frozen orbitals in mixed quantum mechanics/molecular mechanics methods. II. Application of the local basis equation

The application of the local basis equation (Ferenczy and Adams, J. Chem. Phys. 2009 , 130, 134108) in mixed quantum mechanics/molecular mechanics (QM/MM) and quantum mechanics/quantum mechanics (QM/QM) methods is investigated. This equation is suitable to derive local basis nonorthogonal orbitals that minimize the energy of the system and it exhibits good convergence properties in a self‐consistent field solution. These features make the equation appropriate to be used in mixed QM/MM and QM/QM methods to optimize orbitals in the field of frozen localized orbitals connecting the subsystems. Calculations performed for several properties in divers systems show that the method is robust with various choices of the frozen orbitals and frontier atom properties. With appropriate basis set assignment, it gives results equivalent with those of a related approach [G. G. Ferenczy previous paper in this issue] using the Huzinaga equation. Thus, the local basis equation can be used in mixed QM/MM methods with small size quantum subsystems to calculate properties in good agreement with reference Hartree–Fock–Roothaan results. It is shown that bond charges are not necessary when the local basis equation is applied, although they are required for the self‐consistent field solution of the Huzinaga equation based method. Conversely, the deformation of the wave‐function near to the boundary is observed without bond charges and this has a significant effect on deprotonation energies but a less pronounced effect when the total charge of the system is conserved. The local basis equation can also be used to define a two layer quantum system with nonorthogonal localized orbitals surrounding the central delocalized quantum subsystem. © 2013 Wiley Periodicals, Inc.     

Theoretical study of vitamin properties from combined QM-MM methods: Comparison of chemical shifts and energy

The combination of quantum mechanics (QM) and molecular mechanics (MM) methods has become an alternative tool for many applications for which pure QM and MM are not suitable. The QM-MM method has been used for different types of problems, for example, structural biology, surface phenomena, and the liquid phase. In this paper, we have implemented these methods for vitamins, an important kind of biological molecule, and then compared results. The calculations were done by the full ab initio method (HF/3–21 g and HF/6–31 g) and QM-MM (ONIOM) method with HF(3–21 g)/AM1/UFF; then, we found that the geometry obtained by the QM-MM method is very accurate and this rapid method can be used in place of time consuming ab initio methods for large molecules. A comparison of energy values in the QM-MM and QM methods is given. We compare chemical shifts and conclude that the QM-MM method is a perturbed full QM method. The text was submitted by the authors in English.     

An efficient method for the calculation of quantum mechanics/molecular mechanics free energies

The combination of quantum mechanics (QM) with molecular mechanics (MM) offers a route to improved accuracy in the study of biological systems, and there is now significant research effort being spent to develop QM/MM methods that can be applied to the calculation of relative free energies. Currently, the computational expense of the QM part of the calculation means that there is no single method that achieves both efficiency and rigor; either the QM/MM free energy method is rigorous and computationally expensive, or the method introduces efficiency-led assumptions that can lead to errors in the result, or a lack of generality of application. In this paper we demonstrate a combined approach to form a single, efficient, and, in principle, exact QM/MM free energy method. We demonstrate the application of this method by using it to explore the difference in hydration of water and methane. We demonstrate that it is possible to calculate highly converged QM/MM relative free energies at the MP2/aug-cc-pVDZ/OPLS level within just two days of computation, using commodity processors, and show how the method allows consistent, high-quality sampling of complex solvent configurational change, both when perturbing hydrophilic water into hydrophobic methane, and also when moving from a MM Hamiltonian to a QM/MM Hamiltonian. The results demonstrate the validity and power of this methodology, and raise important questions regarding the compatibility of MM and QM/MM forcefields, and offer a potential route to improved compatibility.     

Mixed ab initio QM/MM modeling using frozen orbitals and tests with alanine dipeptide and tetrapeptide

Methodology is discussed for mixed ab initio quantum mechanics/molecular mechanics modeling of systems where the quantum mechanics (QM) and molecular mechanics (MM) regions are within the same molecule. The ab initio QM calculations are at the restricted Hartree–Fock level using the pseudospectral method of the Jaguar program while the MM part is treated with the OPLS force fields implemented in the IMPACT program. The interface between the QM and MM regions, in particular, is elaborated upon, as it is dealt with by “breaking” bonds at the boundaries and using Boys-localized orbitals found from model molecules in place of the bonds. These orbitals are kept frozen during QM calculations. Results from tests of the method to find relative conformational energies and geometries of alanine dipeptides and alanine tetrapeptides are presented along with comparisons to pure QM and pure MM calculations. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1468–1494, 1999     

GARLEEK: Adding an extra flavor to ONIOM

The ONIOM method, developed in the group of Keiji Morokuma, is one of the most successful examples of quantum mechanics/molecular mechanics (QM/MM) treatments, and of multilayer methods in general. Its implementation in the Gaussian program package is in particular widely used. This implementation has access to the wide variety of QM methods available in Gaussian, but is limited to only three specific force fields. The current article presents the GARLEEK interface, which expands the availability of molecular mechanics methods to the wide variety of force fields available in MM packages. The focus is in the simple installation and use. Two examples of the performance of the interface with selected systems are provided. GARLEEK is MIT-licensed and freely available at https://github.com/insilichem/garleek . © 2018 Wiley Periodicals, Inc.     

Multiscale electrostatic embedding simulations for modeling structure and dynamics of molecules in solution: A tutorial review

The main concepts and important technical details of electrostatic embedding quantum mechanics/molecular mechanics (QM/MM) simulations are explained and illustrated with the intent of assisting newcomers in performing and gauging the accuracy of such simulations, focused on smaller molecules in solution. Beginners are advised on how to increase the reliability and accuracy of the simulations through benchmarking. Central considerations on methodologies for QM/MM Molecular Dynamics (MD) simulations are presented, alongside technical fundamentals regarding the construction and manipulation of simulation systems using the python-based Atomic Simulation Environment (ASE). A worked example of QM/MM Born–Oppenheimer MD is included, and a flowchart summarizing the most salient decisions and tasks within the methodology is presented.     

Mechanism of proteolysis in matrix metalloproteinase‐2 revealed by QM/MM modeling

The mechanism of enzymatic peptide hydrolysis in matrix metalloproteinase‐2 (MMP‐2) was studied at atomic resolution through quantum mechanics/molecular mechanics (QM/MM) simulations. An all‐atom three‐dimensional molecular model was constructed on the basis of a crystal structure from the Protein Data Bank (ID: 1QIB), and the oligopeptide Ace‐Gln‐Gly～Ile‐Ala‐Gly‐Nme was considered as the substrate. Two QM/MM software packages and several computational protocols were employed to calculate QM/MM energy profiles for a four‐step mechanism involving an initial nucleophilic attack followed by hydrogen bond rearrangement, proton transfer, and C? N bond cleavage. These QM/MM calculations consistently yield rather low overall barriers for the chemical steps, in the range of 5–10 kcal/mol, for diverse QM treatments (PBE0, B3LYP, and BB1K density functionals as well as local coupled cluster treatments) and two MM force fields (CHARMM and AMBER). It, thus, seems likely that product release is the rate‐limiting step in MMP‐2 catalysis. This is supported by an exploration of various release channels through QM/MM reaction path calculations and steered molecular dynamics simulations. © 2015 Wiley Periodicals, Inc.     

Optimizing efficiency of perturbative Monte Carlo method

We introduce error weighting functions into the perturbative Monte Carlo method for use with a hybrid ab initio quantum mechanics/molecular mechanics (QM/MM) potential. The perturbative Monte Carlo approach introduced earlier provides a means to reduce the number of full SCF calculations in simulations using a QM/MM potential by evoking perturbation theory to calculate energy changes due to displacements of an MM molecule. The use of weighting functions, introduced here, allows an optimal number of MM molecule displacements to occur between the performance of the full self-consistent field calculations. This will allow the ab initio QM/MM approach to be applied to systems that require more accurate treatment of the QM and/or MM regions. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1632–1638, 1998     

Prediction of the mechanisms of enzyme-catalysed reactions using hybrid quantum mechanical/molecular mechanical methods

The use of hybrid methods, involving both quantum mechanics and molecular mechanics, to model the mechanism of enzyme-catalysed reactions, is discussed. Two alternative approaches to treating the electrostatic interactions between the quantum mechanical and molecular mechanical regions are studied, involving either the inclusion of this term in the electronic Hamiltonian (QM/MM), or evaluating it purely classically (MO + MM). In the latter scheme, possible problems of using force fields that are standard for macromolecular modelling are identified. The use of QM/MM schemes to investigate the mechanism of the enzymes thymidine phosphorylase (ThdPase) and protein tyrosine phosphatase (PTP) is described. For both systems, transition states have been identified using a PM3 Hamiltonian. For ThdPase, concerted motion of the enzyme during the course of the reaction is suggested and, for PTP, a two-step dephosphorylation reaction is indicated, both with quite low barriers.