How does the cAMP-dependent protein kinase catalyze the phosphorylation reaction: an ab initio QM/MM study

We have carried out density functional theory QM/MM calculations on the catalytic subunit of cAMP-dependent protein kinase (PKA). The QM/MM calculations indicate that the phosphorylation reaction catalyzed by PKA is mainly dissociative, and Asp166 serves as the catalytic base to accept the proton delivered by the substrate peptide. Among the key interactions in the active site, the Mg(2+) ions, glycine rich loop, and Lys72 are found to stabilize the transition state through electrostatic interactions. On the other hand, Lys168, Asn171, Asp184, and the conserved waters bound to Mg(2+) ions do not directly contribute to lower the energy barrier of the phosphorylation reaction, and possible roles for these residues are proposed. The QM/MM calculations with different QM/MM partition schemes or different initial structures yield consistent results. In addition, we have carried out 12 ns molecular dynamics simulations on both wild type and K168A mutated PKA, respectively, to demonstrate that the catalytic role of Lys168 is to keep ATP and substrate peptide in the near-attack reactive conformation.     

Convergence in the QM‐only and QM/MM modeling of enzymatic reactions: A case study for acetylene hydratase

We report systematic quantum mechanics‐only (QM‐only) and QM/molecular mechanics (MM) calculations on an enzyme‐catalyzed reaction to assess the convergence behavior of QM‐only and QM/MM energies with respect to the size of the chosen QM region. The QM and MM parts are described by density functional theory (typically B3LYP/def2‐SVP) and the CHARMM force field, respectively. Extending our previous work on acetylene hydratase with QM regions up to 157 atoms (Liao and Thiel, J. Chem. Theory Comput. 2012, 8, 3793), we performed QM/MM geometry optimizations with a QM region M4 composed of 408 atoms, as well as further QM/MM single‐point calculations with even larger QM regions up to 657 atoms. A charge deletion analysis was conducted for the previously used QM/MM model ( M3a , with a QM region of 157 atoms) to identify all MM residues with strong electrostatic contributions to the reaction energetics (typically more than 2 kcal/mol), which were then included in M4 . QM/MM calculations with this large QM region M4 lead to the same overall mechanism as the previous QM/MM calculations with M3a , but there are some variations in the relative energies of the stationary points, with a mean absolute deviation (MAD) of 2.7 kcal/mol. The energies of the two relevant transition states are close to each other at all levels applied (typically within 2 kcal/mol), with the first (second) one being rate‐limiting in the QM/MM calculations with M3a ( M4 ). QM‐only gas‐phase calculations give a very similar energy profile for QM region M4 (MAD of 1.7 kcal/mol), contrary to the situation for M3a where we had previously found significant discrepancies between the QM‐only and QM/MM results (MAD of 7.9 kcal/mol). Extension of the QM region beyond M4 up to M7 (657 atoms) leads to only rather small variations in the relative energies from single‐point QM‐only and QM/MM calculations (MAD typically about 1–2 kcal/mol). In the case of acetylene hydratase, a model with 408 QM atoms thus seems sufficient to achieve convergence in the computed relative energies to within 1–2 kcal/mol.Copyright © 2013 Wiley Periodicals, Inc.     

QM/MM analysis suggests that Alkaline Phosphatase (AP) and nucleotide pyrophosphatase/phosphodiesterase slightly tighten the transition state for phosphate diester hydrolysis relative to solution: implication for catalytic promiscuity in the AP superfamily

Several members of the Alkaline Phosphatase (AP) superfamily exhibit a high level of catalytic proffciency and promiscuity in structurally similar active sites. A thorough characterization of the nature of transition state for different substrates in these enzymes is crucial for understanding the molecular mechanisms that govern those remarkable catalytic properties. In this work, we study the hydrolysis of a phosphate diester, MpNPP(-), in solution, two experimentally well-characterized variants of AP (R166S AP, R166S/E322Y AP) and wild type Nucleotide pyrophosphatase/phosphodiesterase (NPP) by QM/MM calculations in which the QM method is an approximate density functional theory previously parametrized for phosphate hydrolysis (SCC-DFTBPR). The general agreements found between these calculations and available experimental data for both solution and enzymes support the use of SCC-DFTBPR/MM for a semiquantitative analysis of the catalytic mechanism and nature of transition state in AP and NPP. Although phosphate diesters are cognate substrates for NPP but promiscuous substrates for AP, the calculations suggest that their hydrolysis reactions catalyzed by AP and NPP feature similar synchronous transition states that are slightly tighter in nature compared to that in solution, due in part to the geometry of the bimetallic zinc motif. Therefore, this study provides the first direct computational support to the hypothesis that enzymes in the AP superfamily catalyze cognate and promiscuous substrates via similar transition states to those in solution. Our calculations do not support the finding of recent QM/MM studies by López-Canut and co-workers, who suggested that the same diester substrate goes through a much looser transition state in NPP/AP than in solution, a result likely biased by the large structural distortion of the bimetallic zinc site in their simulations. Finally, our calculations for different phosphate diester orientations and phosphorothioate diesters highlight that the interpretation of thio-substitution experiments is not always straightforward.     

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.     

Development of massive multilevel molecular dynamics simulation program,platypus (PLATform for dYnamic protein unified simulation), for the elucidation of protein functions

A massively parallel program for quantum mechanical‐molecular mechanical (QM/MM) molecular dynamics simulation, called Platypus (PLATform for dYnamic Protein Unified Simulation), was developed to elucidate protein functions. The speedup and the parallelization ratio of Platypus in the QM and QM/MM calculations were assessed for a bacteriochlorophyll dimer in the photosynthetic reaction center (DIMER) on the K computer, a massively parallel computer achieving 10 PetaFLOPs with 705,024 cores. Platypus exhibited the increase in speedup up to 20,000 core processors at the HF/cc‐pVDZ and B3LYP/cc‐pVDZ, and up to 10,000 core processors by the CASCI(16,16)/6‐31G** calculations. We also performed excited QM/MM‐MD simulations on the chromophore of Sirius (SIRIUS) in water. Sirius is a pH‐insensitive and photo‐stable ultramarine fluorescent protein. Platypus accelerated on‐the‐fly excited‐state QM/MM‐MD simulations for SIRIUS in water, using over 4000 core processors. In addition, it also succeeded in 50‐ps (200,000‐step) on‐the‐fly excited‐state QM/MM‐MD simulations for the SIRIUS in water. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

A water-mediated and substrate-assisted catalytic mechanism for Sulfolobus solfataricus DNA polymerase IV

DNA polymerases are enzymes responsible for the synthesis of DNA from nucleotides. Understanding their molecular fundamentals is a prerequisite for elucidating their aberrant activities in diseases such as cancer. Here we have carried out ab initio quantum mechanical/molecular mechanical (QM/MM) studies on the nucleotidyl-transfer reaction catalyzed by the lesion-bypass DNA polymerase IV (Dpo4) from Sulfolobus solfataricus, with template guanine and Watson-Crick paired dCTP as the nascent base pair. The results suggested a novel water-mediated and substrate-assisted (WMSA) mechanism: the initial proton transfer to the alpha-phosphate of the substrate via a bridging crystal water molecule is the rate-limiting step, the nucleotidyl-transfer step is associative with a metastable pentacovalent phosphorane intermediate, and the pyrophosphate leaving is facilitated by a highly coordinated proton relay mechanism through mediation of water which neutralizes the evolving negative charge. The conserved carboxylates, which retain their liganding to the two Mg2+ ions during the reaction process, are found to be essential in stabilizing transition states. This WMSA mechanism takes specific advantage of the unique structural features of this low-fidelity lesion-bypass Y-family polymerase, which has a more spacious and solvent-exposed active site than replicative and repair polymerases.     

pKa calculations in solution and proteins with QM/MM free energy perturbation simulations: a quantitative test of QM/MM protocols

The accuracy of biological simulations depends, in large part, on the treatment of electrostatics. Due to the availability of accurate experimental values, calculation of pKa provides stringent evaluation of computational methods. The generalized solvent boundary potential (GSBP) and Ewald summation electrostatic treatments were recently implemented for combined quantum mechanical and molecular mechanics (QM/MM) simulations by our group. These approaches were tested by calculating pKa shifts due to differences in electronic structure and electrostatic environment; the shifts were determined for a series of small molecules in solution, using various electrostatic treatments, and two residues (His 31, Lys 102) in the M102K T4-lysozyme mutant with large pKa shifts, using the GSBP approach. The calculations utilized a free energy perturbation scheme with the QM/MM potential function involving the self-consistent charge density functional tight binding (SCC-DFTB) and CHARMM as the QM and MM methods, respectively. The study of small molecules demonstrated that inconsistent electrostatic models produced results that were difficult to correct in a robust manner; by contrast, extended electrostatics, GSBP, and Ewald simulations produced consistent results once a bulk solvation contribution was carefully chosen. In addition to the electrostatic treatment, the pKa shifts were also sensitive to the level of the QM method and the scheme of treating QM/MM Coulombic interactions; however, simple perturbative corrections based on SCC-DFTB/CHARMM trajectories and higher level single point energy calculations were found to give satisfactory results. Combining all factors gave a root-mean-square difference of 0.7 pKa units for the relative pKa values of the small molecules compared to experiment. For the residues in the lysozyme, an accurate pKa shift was obtained for His 31 with multiple nanosecond simulations. For Lys 102, however, the pKa shift was estimated to be too large, even after more than 10 nanosecond simulations for each lambda window; the difficulty was due to the significant, but slow, reorganization of the protein and water structure when Lys 102 was protonated. The simulations support that Lys 102 is deprotonated in the X-ray structure and the protein is highly destabilized when this residue is protonated.     

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.     

Systematic QM/MM investigation of factors that affect the cytochrome P450-catalyzed hydrogen abstraction of camphor

The hydrogen abstraction reaction of camphor in cytochrome P450(cam) has been investigated in the native enzyme environment by combined quantum mechanical/molecular mechanical (QM/MM) calculations and in the gas phase by density functional calculations. This work has been motivated by contradictory published QM/MM results. In an attempt to pinpoint the origin of these discrepancies, we have systematically studied the factors that may affect the computed barriers, including the QM/MM setup, the optimization procedures, and the choice of QM region, basis set, and protonation states. It is found that the ChemShell and QSite programs used in the published QM/MM calculations yield similar results at given geometries, and that the discrepancies mainly arise from two technical issues (optimization protocols and initial system preparation) that need to be well controlled in QM/MM work. In the course of these systematic investigations, new mechanistic insights have been gained. The crystallographic water 903 placed near the oxo atom of Compound I lowers the hydrogen abstraction barrier by ca. 4 kcal/mol, and thus acts as a catalyst for this reaction. Spin density may appear at the A-propionate side chain of the heme if the carboxylate group is not properly screened, which might be expected to happen during protein dynamics, but not in static equilibrium situations. There is no clear correlation between the computed A-propionate spin density and the hydrogen abstraction barrier, and hence, no support for a previously proposed side-chain mediated transition state stabilization mechanism. Standard QM/MM optimizations yield an A-propionate environment close to the X-ray structure only for protonated Asp297, and not for deprotonated Asp297, but the computed barriers are similar in both cases. An X-ray like A-propionate environment can also be obtained when deprotonated Asp297 is included in the QM region and His355 is singly protonated, but this Compound II-type species with a closed-shell porphyrin ring has a higher hydrogen abstraction barrier and should thus not be mechanistically relevant.     

A quantum mechanical/molecular mechanical study on the catalysis of the pyridoxal 5'-phosphate-dependent enzyme L-serine dehydratase

The catalytic mechanism of a pyridoxal 5'-phosphate-dependent enzyme, l-serine dehydratase, has been investigated using ab initio quantum mechanical/molecular mechanical (QM/MM) methods. New insights into the chemical steps have been obtained, including the chemical role of the substrate carboxyl group in the Schiff base formation step and a proton-relaying mechanism involving the phosphate of the cofactor in the beta-hydroxyl-leaving step. The latter step is of no barrier and follows sequentially after the elimination of the alpha-proton, leading to a single but sequential alpha, beta-elimination step. The rate-limiting transition state is specifically stabilized by the enzyme environment. At this transition state, charges are localized on the substrate carboxyl group, as well as on the amino group of Lys41. Specific interactions of the enzyme environment with these groups are able to lower the activation barrier significantly. One major difficulty associated with studies of complicated enzymatic reactions using ab initio QM/MM models is the appropriate choices of reaction coordinates. In this study, we have made use of efficient semiempirical models and pathway optimization techniques to overcome this difficulty.     

Theoretical perspectives on the reaction mechanism of serine proteases: the reaction free energy profiles of the acylation process

The reaction mechanism of serine proteases (trypsin), which catalyze peptide hydrolysis, is studied theoretically by ab initio QM/MM electronic structure calculations combined with Molecular Dynamics-Free Energy Perturbation calculations. We have calculated the entire reaction free energy profiles of the first reaction step of this enzyme (acylation process). The present calculations show that the rate-determining step of the acylation is the formation of the tetrahedral intermediate, and the breakdown of this intermediate has a small energy barrier. The calculated activation free energy for the acylation is approximately 17.8 kcal/mol at QM/MM MP2/(aug)-cc-pVDZ//HF/6-31(+)G/AMBER level, and this reaction is an exothermic process. MD simulations of the enzyme-substrate (ES) complex and the free enzyme in aqueous phase show that the substrate binding induces slight conformational changes around the active site, which favor the alignment of the reactive fragments (His57, Asp102, and Ser195) together in a reactive orientation. It is also shown that the proton transfer from Ser195 to His57 and the nucleophilic attack of Ser195 to the carbonyl carbon of the scissile bond of the substrate occur in a concerted manner. In this reaction, protein environment plays a crucial role to lowering the activation free energy by stabilizing the tetrahedral intermediate compared to the ES complex. The polarization energy calculations show that the enzyme active site is in a very polar environment because of the polar main chain contributions of protein. Also, the ground-state destabilization effect (steric strain) is not a major catalytic factor. The most important catalytic factor of stabilizing the tetrahedral intermediate is the electrostatic interaction between the active site and particular regions of protein: the main chain NH groups in Gly193 and Ser195 (so-called oxyanion hole region) stabilize negative charge generated on the carbonyl oxygen of the scissile bond, and the main chain carbonyl groups in Ile212 approximately Ser214 stabilize a positive charge generated on the imidazole ring of His57.     

Recent advances in quantum mechanical/molecular mechanical calculations of enzyme catalysis: hydrogen tunnelling in liver alcohol dehydrogenase and inhibition of elastase by α-ketoheterocycles

 Hybrid quantum mechanical (QM)/molecular mechanical (MM) calculations are used to study two aspects of enzyme catalysis, Kinetic isotope effects associated with the hydride ion transfer step in the reduction of benzyl alcohol by liver alcohol dehydrogenase are studied by employing variational transition-state theory and optimised multidimensional tunnelling. With the smaller QM region, described at the Hartree–Fock ab initio level, together with a parameterised zinc atom charge, good agreement with experiment is obtained. A comparison is made with the proton transfer in methylamine dehydrogenase. The origin of the large range in pharmacological activity shown by a series of α-ketoheterocycle inhibitors of the serine protease, elastase, is investigated by both force field and QM/MM calculations. Both models point to two different inhibition mechanisms being operative. Initial QM/MM calculations suggest that these are binding, and reaction to form a tetrahedral intermediate, the latter process occurring for only the more potent set of inhibitors. Recieved 3 October 2001 / Accepted: 6 September 2002 / Published online: 31 January 2003 Contribution to the Proceedings of the Symposium on Combined QM/MM Methods at the 222nd National Meeting of the American Chemical Society, 2001 Correspondence to: I. H. Hillier Acknowledgements. We thank EPSRC and BBSRC for support of the research and D.G. Truhlar for the use of the POLYRATE code.     

Whether proton transition to the triphosphate tail of ATP occurs at protein kinase environment: A Car‐Parrinello ab initio molecular dynamics study

Protonation state of the triphosphate tail of ATP (adenosine triphosphate) in protein environment is a fundamental issue, which has significant impact on the mechanism investigation of biochemical processes with ATP involved. Proton transition from surroundings (water molecule coordinating to magnesium, H_W; amino group of Lys, H_L) to the ATP tail in the catalytic core of protein kinase found recently disproved the commonly accepted deprotonation state of ATP tail. In this account, Car‐Parrinello ab initio molecular dynamics (CP‐AIMD) method has been employed to examine whether the proton transition occurs. To provide a comparison basis for the dynamics simulations, static quantum mechanics (QM), and combined quantum mechanics and molecular mechanics (QM/MM) calculations have also been carried out. Consistent results have been obtained that complete transition of hydrogen from the surroundings to the triphosphate tail of ATP is not allowed. The most dominant conformations correspond to the ones with H_W bonding to O(W) and H‐bonding to O(ATP), [O(W)‐H_W···O(ATP)], H_L bonding to N(Lys) and H‐bonding to O(ATP), [N(Lys)‐H_L···O(ATP)]. Metastable structures with one hydrogen atom bonding with two heavy atoms (hydrogen acceptors) were also located by our dynamic simulations. This bonding mode can satisfy the hungering for hydrogen of the two heavy atoms simultaneously. © 2008 Wiley Periodicals, Inc. Int J Quantum Chem, 2008     

Role of the catalytic triad and oxyanion hole in acetylcholinesterase catalysis: an ab initio QM/MM study

The initial step of the acylation reaction catalyzed by acetylcholinesterase (AChE) has been studied by a combined ab initio quantum mechanical/molecular mechanical (QM/MM) approach. The reaction proceeds through the nucleophilic addition of the Ser203 O to the carbonyl C of acetylcholine, and the reaction is facilitated by simultaneous proton transfer from Ser203 to His447. The calculated potential energy barrier at the MP2(6-31+G) QM/MM level is 10.5 kcal/mol, consistent with the experimental reaction rate. The third residue of the catalytic triad, Glu334, is found to be essential in stabilizing the transition state through electrostatic interactions. The oxyanion hole, formed by peptidic NH groups from Gly121, Gly122, and Ala204, is also found to play an important role in catalysis. Our calculations indicate that, in the AChE-ACh Michaelis complex, only two hydrogen bonds are formed between the carbonyl oxygen of ACh and the peptidic NH groups of Gly121 and Gly122. As the reaction proceeds, the distance between the carbonyl oxygen of ACh and NH group of Ala204 becomes smaller, and the third hydrogen bond is formed both in the transition state and in the tetrahedral intermediate.     

All electron quantum chemical calculation of the entire enzyme system confirms a collective catalytic device in the chorismate mutase reaction

To elucidate the catalytic power of enzymes, we analyzed the reaction profile of Claisen rearrangement of Bacillus subtilis chorismate mutase (BsCM) by all electron quantum chemical calculations using the fragment molecular orbital (FMO) method. To the best of our knowledge, this is the first report of ab initio-based quantum chemical calculations of the entire enzyme system, where we provide a detailed analysis of the catalytic factors that accomplish transition-state stabilization (TSS). FMO calculations deliver an ab initio-level estimate of the intermolecular interaction between the substrate and the amino acid residues of the enzyme. To clarify the catalytic role of Arg90, we calculated the reaction profile of the wild-type BsCM as well as Lys90 and Cit90 mutant BsCMs. Structural refinement and the reaction path determination were performed at the ab initio QM/MM level, and FMO calculations were applied to the QM/MM refined structures. Comparison between three types of reactions established two collective catalytic factors in the BsCM reaction: (1) the hydrogen bonds connecting the Glu78-Arg90-substrate cooperatively control the stability of TS relative to the ES complex and (2) the positive charge on Arg90 polarizes the substrate in the TS region to gain more electrostatic stabilization.     

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.