Hydrolysis of phosphotriesters: a theoretical analysis of the enzymatic and solution mechanisms

A theoretical study on the alkaline hydrolysis of paraoxon, one of the most popular organophosphorus pesticides, in aqueous solution and in the active site of Pseudomonas diminuta phosphotriesterase (PTE) is presented. Simulations by means of hybrid quantum mechanics/molecular mechanics (QM/MM) potentials show that the hydrolysis of paraoxon takes place through an A(N)D(N) or associative mechanism both in solution and in the active site of PTE. The results correctly reproduce the magnitude of the activation free energies and can be used to rationalize the observed kinetic isotope effects (KIEs) for the hydrolysis of paraoxon in both media. Enzymatic hydrolysis of O,O-diethyl p-chlorophenyl phosphate, a phosphotriester having a leaving group with higher pK(a) than paraoxon, was also simulated. Hydrolysis of this phosphotriester by PTE follows a A(N)+D(N) mechanism with a pentacoordinate intermediate. Moreover, the leaving group of this new substrate coordinates to one of the zinc ions of the bimetallic active site in order to stabilize the large negative charge developed on the oxygen atom of the leaving group when the P-O bond is broken in the products state. To accommodate this new ligand in the coordination shell, carbamylated Lys169 must be displaced from one zinc ion to the other, which in turn affects the acidity of Asp301, a residue originally bound to the second zinc ion. This ability to displace some of the ligands of the coordination shell of the zinc centers would explain the promiscuity of this enzyme, which is capable of catalyzing hydrolysis of different substrate by means of different mechanisms.     

The implementation of a fast and accurate QM/MM potential method in Amber

Version 9 of the Amber simulation programs includes a new semi-empirical hybrid QM/MM functionality. This includes support for implicit solvent (generalized Born) and for periodic explicit solvent simulations using a newly developed QM/MM implementation of the particle mesh Ewald (PME) method. The code provides sufficiently accurate gradients to run constant energy QM/MM MD simulations for many nanoseconds. The link atom approach used for treating the QM/MM boundary shows improved performance, and the user interface has been rewritten to bring the format into line with classical MD simulations. Support is provided for the PM3, PDDG/PM3, PM3CARB1, AM1, MNDO, and PDDG/MNDO semi-empirical Hamiltonians as well as the self-consistent charge density functional tight binding (SCC-DFTB) method. Performance has been improved to the point where using QM/MM, for a QM system of 71 atoms within an explicitly solvated protein using periodic boundaries and PME requires less than twice the cpu time of the corresponding classical simulation.     

New Insights into the Jahn–Teller Effect through ab initio Quantum-Mechanical/Molecular-Mechanical Molecular Dynamics Simulations of CuII in Water

The Cu^II hydration shell structure has been studied by means of classical molecular dynamics (MD) simulations including three-body corrections and hybrid quantum-mechanical/molecular-mechanical (QM/MM) molecular dynamics (MD) simulations at the Hartree–Fock level. The copper(II ) ion is found to be six-fold coordinated and [Cu(H₂O)₆]²⁺ exhibits a distorted octahedral structure. The QM/MM MD approach reproduces correctly the experimentally observed Jahn–Teller effect but exhibits faster inversions (<200 fs) and a more complex behaviour than expected from experimental data. The dynamic Jahn–Teller effect causes the high lability of [Cu(H₂O)₆]²⁺ with a ligand-exchange rate constant some orders of magnitude higher than its neighbouring ions Ni^II and Zn^II. Nevertheless, no first-shell water exchange occurred during a 30-ps simulation. The structure of the hydrated ion is discussed in terms of radial distribution functions, coordination numbers, and various angular distributions and the dynamical properties as librational and vibrational motions and reorientational times were evaluated, which lead to detailed information about the first hydration shell. Second-shell water-exchange processes could be observed within the simulation time scale and yielded a mean ligand residence time of ≈20 ps.     

Classical and quantum mechanical/molecular mechanical molecular dynamics simulations of alanine dipeptide in water: comparisons with IR and vibrational circular dichroism spectra

We have implemented the combined quantum mechanical (QM)/molecular mechanical (MM) molecular dynamics (MD) simulations of alanine dipeptide in water along with the polarizable and nonpolarizable classical MD simulations with different models of water. For the QM/MM MD simulation, the alanine dipeptide is treated with the AM1 or PM3 approximations and the fluctuating solute dipole moment is calculated by the Mulliken population analysis. For the classical MD simulations, the solute is treated with the polarizable or nonpolarizable AMBER and polarizable CHARMM force fields and water is treated with the TIP3P, TIP4P, or TIP5P model. It is found that the relative populations of right-handed alpha-helix and extended beta and P(II) conformations in the simulation trajectory strongly depend on the simulation method. For the QM/MM MD simulations, the PM3/MM shows that the P(II) conformation is dominant, whereas the AM1/MM predicts that the dominant conformation is alpha(R). Polarizable CHARMM force field gives almost exclusively P(II) conformation and other force fields predict that both alpha-helical and extended (beta and P(II)) conformations are populated with varying extents. Solvation environment around the dipeptide is investigated by examining the radial distribution functions and numbers and lifetimes of hydrogen bonds. Comparing the simulated IR and vibrational circular dichroism spectra with experimental results, we concluded that the dipeptide adopts the P(II) conformation and PM3/MM, AMBER03 with TIP4P water, and AMBER polarizable force fields are acceptable for structure determination of the dipeptide considered in this paper.     

The CHARMM–TURBOMOLE interface for efficient and accurate QM/MM molecular dynamics,free energies,and excited state properties

The quantum mechanical (QM)/molecular mechanical (MM) interface between Chemistry at HARvard Molecular Mechanics (CHARMM) and TURBOMOLE is described. CHARMM provides an extensive set of simulation algorithms, like molecular dynamics (MD) and free energy perturbation, and support for mature nonpolarizable and Drude polarizable force fields. TURBOMOLE provides fast QM calculations using density functional theory or wave function methods and excited state properties. CHARMM–TURBOMOLE is well‐suited for extended QM/MM MD simulations using first principles methods with large (triple‐ζ) basis sets. We demonstrate these capabilities with a QM/MM simulation of Mg²⁺(aq), where the MM outer sphere water molecules are represented using the SWM4‐NDP Drude polarizable force field and the ion and inner coordination sphere are represented using QM PBE, PBE0, and MP2 methods. The relative solvation free energies of Mg²⁺ and Zn²⁺ were calculated using thermodynamic integration. We also demonstrate the features for excited state properties. We calculate the time‐averaged solution absorption spectrum of indole, the emission spectrum of the indole excited state, and the electronic circular dichroism spectrum of an oxacepham. © 2014 Wiley Periodicals, Inc.     

Coordination number of zinc ions in the phosphotriesterase active site by molecular dynamics and quantum mechanics

We have run several molecular dynamics (MD) simulations on zinc-containing phosphotriesterase (PTE) with two bound substrates, sarin and paraoxon, and with the substrate analog diethyl 4-methylbenzylphosphonate. A standard nonbonded model was employed to treat the zinc ions with the commonly used charge of +2. In all the trajectories, we observed a tightly bound water (TBW) molecule in the active site that was coordinated to the less buried zinc ion. The phosphoryl oxygen of the substrate/inhibitor was found to be coordinated to the same zinc ion so that, considering all ligands, the less buried zinc was hexa-coordinated. The hexa-coordination of this zinc ion was not seen in the deposited X-ray pdb files for PTE. Several additional MD simulations were then performed using different charges (+1, +1.5) on the zinc ions, along with ab initio and density functional theory (DFT) calculations, to evaluate the following possibilities: the crystal diffraction data were not correctly interpreted; the hexa-coordinated zinc ion in PTE is only present in solution and not in the crystal; and the hexa-coordinated zinc ion in PTE is an artifact of the force field used. A charge of +1.5 leads to a coordination number (CN) of 5 on both zinc ions, which is consistent with the results from ab initio and DFT calculations and with the latest high resolution X-ray crystal structure. The commonly used charge of +2 produces a CN of 6 on the less buried zinc. The CN on the more buried zinc ion is 5 when the substrate/inhibitor is present in the simulation, and increases to 6 when the substrate/inhibitor is removed prior to the simulation. The results of both of the MD and quantum mechanical calculations lead to the conclusion that the zinc ions in the PTE active site are both penta-coordinated, and that the MD simulations performed with the charge of +2 overestimate the CN of the zinc ions in the PTE active site. The overall protein structures in the simulations remain unaffected by the change in zinc charge from +2 to +1.5. The results also suggest that the charge +1.5 is the most appropriate for the molecular dynamics simulations on zinc-containing PTE when a nonbonded model is used and no global thermodynamic conclusion is sought. We also show that the standard nonbonded model is not able to properly treat the CN and energy at the same time. A preliminary, promising charge-transfer model is discussed with the use of the zinc charge of +1.5.     

MM and QM/MM Modeling of Threonyl-tRNA Synthetase: Model Testing and Simulations

Aminoacyl-tRNA synthetases are centrally important enzymes in protein synthesis. We have investigated threonyl-tRNA synthetase from E. coli, complexed with reactants, using molecular mechanics and combined quantum mechanical/molecular mechanical (QM/MM) techniques. These modeling methods have the potential to provide molecular level understanding of enzyme catalytic processes. Modeling of this enzyme presents a number of challenges. The procedure of system preparation and testing is described in detail. For example, the number of metal ions at the active site, and their positions, were investigated. Molecular dynamics simulations suggest that the system is most stable when it contains only one magnesium ion, and the zinc ion is removed. Two different QM/MM methods were tested in models based on the findings of MM molecular dynamics simulations. AM1/CHARMM calculations resulted in unrealistic structures for the phosphates in this system. This is apparently due to an error of AM1. PM3/CHARMM calculations proved to be more suitable for this enzyme system. These results will provide a useful basis for future modeling investigations of the enzyme mechanism and dynamics.     

Investigation of the metal binding site in methionine aminopeptidase by density functional theory

All methionine aminopeptidases exhibit the same conserved metal binding site. The structure of this site with either Co²⁺ions or Zn²⁺ions was investigated using density functional theory. The calculations showed that the structure of the site was not influenced by the identity of the metal ions. This was the case for both of the systems studied; one based on the X-ray structure of the human methionine aminopeptidase type 2 (hMetAP-2) and the other based on the X-ray structure of the E. colimethionine aminopeptidase type 1 (eMetAP-1). Another important structural issue is the identity of the bridging oxygen, which is part of either a water molecule or a hydroxide ion. Within the site of hMetAP-2 the results strongly indicate that a hydroxide ion bridges the metal ions. By contrast, the nature of the oxygen bridging the metal ions within the metal binding site of eMetAP-1 cannot be determined based on the results here, due to the similar structural results obtained with a bridging water molecule and a bridging hydroxide ion.     

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.     

A valence bond model for aqueous Cu(II) and Zn(II) ions in the AMOEBA polarizable force field

A general molecular mechanics (MM) model for treating aqueous Cu²⁺ and Zn²⁺ ions was developed based on valence bond (VB) theory and incorporated into the atomic multipole optimized energetics for biomolecular applications (AMOEBA) polarizable force field. Parameters were obtained by fitting MM energies to that computed by ab initio methods for gas‐phase tetra‐ and hexa‐aqua metal complexes. Molecular dynamics (MD) simulations using the proposed AMOEBA‐VB model were performed for each transition metal ion in aqueous solution, and solvent coordination was evaluated. Results show that the AMOEBA‐VB model generates the correct square‐planar geometry for gas‐phase tetra‐aqua Cu²⁺ complex and improves the accuracy of MM model energetics for a number of ligation geometries when compared to quantum mechanical (QM) computations. On the other hand, both AMOEBA and AMOEBA‐VB generate results for Zn²⁺–water complexes in good agreement with QM calculations. Analyses of the MD trajectories revealed a six‐coordination first solvation shell for both Cu²⁺ and Zn²⁺ ions in aqueous solution, with ligation geometries falling in the range reported by previous studies. © 2012 Wiley Periodicals, Inc.     

Antibiotic deactivation by a dizinc beta-lactamase: mechanistic insights from QM/MM and DFT studies

Hybrid quantum mechanical/molecular mechanical (QM/MM) methods and density functional theory (DFT) were used to investigate the initial ring-opening step in the hydrolysis of moxalactam catalyzed by the dizinc L1 beta-lactamase from Stenotrophomonas maltophilia. Anchored at the enzyme active site via direct metal binding as suggested by a recent X-ray structure of an enzyme-product complex (Spencer, J.; et al. J. Am. Chem. Soc. 2005, 127, 14439), the substrate is well aligned with the nucleophilic hydroxide that bridges the two zinc ions. Both QM/MM and DFT results indicate that the addition of the hydroxide nucleophile to the carbonyl carbon in the substrate lactam ring leads to a metastable intermediate via a dominant nucleophilic addition barrier. The potential of mean force obtained by SCC-DFTB/MM simulations and corrected by DFT/MM calculations yields a reaction free energy barrier of 23.5 kcal/mol, in reasonable agreement with the experimental value of 18.5 kcal/mol derived from kcat of 0.15 s(-1). It is further shown that zinc-bound Asp120 plays an important role in aligning the nucleophile, but accepts the hydroxide proton only after the nucleophilic addition. The two zinc ions are found to participate intimately in the catalysis, consistent with the proposed mechanism. In particular, the Zn(1) ion is likely to serve as an "oxyanion hole" in stabilizing the carbonyl oxygen, while the Zn(2) ion acts as an electrophilic catalyst to stabilize the anionic nitrogen leaving group.     

Phosphorylation reaction in cAPK protein kinase-free energy quantum mechanical/molecular mechanics simulations

We present results of a theoretical analysis of the phosphorylation reaction in cAMP-dependent protein kinase using a combined quantum mechanical and molecular mechanics (QM/MM) approach. Detailed analysis of the reaction pathway is provided using a novel QM/MM implementation of the nudged elastic band method, finite temperature fluctuations of the protein environment are taken into account using free energy calculations, and an analysis of hydrogen bond interactions is performed on the basis of calculated frequency shifts. The late transfer of the substrate proton to the conserved aspartate (D166), the activation free energy of 15 kcal/mol, and the slight exothermic (-3 kcal/mol) character of the reaction are all consistent with the experimental data. The near attack conformation of D166 in the reactant state is maintained by interactions with threonine-201, asparagine-177, and most notably by a conserved water molecule serving as a strong structural link between the primary metal ion and the D166. The secondary Mg ion acts as a Lewis acid, attacking the beta-gamma bridging oxygen of ATP. This interaction, along with a strong hydrogen bond between the D166 and the substrate, contributes to the stabilization of the transition state. Lys-168 maintains a hydrogen bond to a transferring phosphoryl group throughout a reaction process. This interaction increases in the product state and contributes to its stabilization.     

Smooth solvation method for d-orbital semiempirical calculations of biological reactions. 2. Application to transphosphorylation thio effects in solution

Density-functional and semiempirical quantum methods and continuum dielectric and explicit solvation models are applied to study the role of solvation on the stabilization of native and thio-substituted transphosphorylation reactions. Extensive comparison is made between results obtained from the different methods. For the semiempirical methods, explicit solvation was treated using a hybrid quantum mechanical/molecular mechanical (QM/MM) approach and the implicit solvation was treated using a recently developed smooth solvation model implemented into a d-orbital semiempirical framework (MNDO/d-SCOSMO) within CHARMM. The different quantum and solvation methods were applied to the transesterification of a 3'-ribose,5'-methyl phosphodiester that serves as a nonenzymatic model for the self-cleavage reaction catalyzed by the hammerhead and hairpin ribozymes. Thio effects were studied for a double sulfur substitution at the nonbridging phosphoryl oxygen positions. The reaction profiles of both the native and double sulfur-substituted reactions from the MNDO/d-SCOSMO calculations were similar to the QM/MM results and consistent with the experimentally observed trends. These results underscore the need for a d-orbital semiempirical representation for phosphorus and sulfur for the study of experimentally observed thio effects in enzymatic and nonenzymatic phosphoryl transfer reactions. One of the major advantages of the present approach is that it can be applied to model chemical reactions at a significantly lower computational cost than either the density-functional calculations with implicit solvation or the semiempirical QM/MM simulations with explicit solvent.     

CPMD/GULP QM/MM interface for modeling periodic solids: Implementation and its application in the study of Y‐zeolite supported Rhn clusters

We report here the development of hybrid quantum mechanics/molecular mechanics (QM/MM) interface between the plane‐wave density functional theory based CPMD code and the empirical force‐field based GULP code for modeling periodic solids and surfaces. The hybrid QM/MM interface is based on the electrostatic coupling between QM and MM regions. The interface is designed for carrying out full relaxation of all the QM and MM atoms during geometry optimizations and molecular dynamics simulations, including the boundary atoms. Both Born–Oppenheimer and Car–Parrinello molecular dynamics schemes are enabled for the QM part during the QM/MM calculations. This interface has the advantage of parallelization of both the programs such that the QM and MM force evaluations can be carried out in parallel to model large systems. The interface program is first validated for total energy conservation and parallel scaling performance is benchmarked. Oxygen vacancy in α‐cristobalite is then studied in detail and the results are compared with a fully QM calculation and experimental data. Subsequently, we use our implementation to investigate the structure of rhodium cluster (Rh_n; n = 2 to 6) formed from Rh(C₂H₄)₂ complex adsorbed within a cavity of Y‐zeolite in a reducible atmosphere of H₂ gas. © 2016 Wiley Periodicals, Inc.     

Increasing the time step with mass scaling in Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations

Born‐Oppenheimer ab initio QM/MM molecular dynamics simulation with umbrella sampling is a state‐of‐the‐art approach to calculate free energy profiles of chemical reactions in complex systems. To further improve its computational efficiency, a mass‐scaling method with the increased time step in MD simulations has been explored and tested. It is found that by increasing the hydrogen mass to 10 amu, a time step of 3 fs can be employed in ab initio QM/MM MD simulations. In all our three test cases, including two solution reactions and one enzyme reaction, the resulted reaction free energy profiles with 3 fs time step and mass scaling are found to be in excellent agreement with the corresponding simulation results using 1 fs time step and the normal mass. These results indicate that for Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, the mass‐scaling method can significantly reduce its computational cost while has little effect on the calculated free energy profiles. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Role of Mg2+ in hammerhead ribozyme catalysis from molecular simulation

Molecular dynamics simulations have been performed to investigate the role of Mg2+ in the full-length hammerhead ribozyme cleavage reaction. In particular, the aim of this work is to characterize the binding mode and conformational events that give rise to catalytically active conformations and stabilization of the transition state. Toward this end, a series of eight 12 ns molecular dynamics simulations have been performed with different divalent metal binding occupations for the reactant, early and late transition state using recently developed force field parameters for metal ions and reactive intermediates in RNA catalysis. In addition, hybrid QM/MM calculations of the early and late transition state were performed to study the proton-transfer step in general acid catalysis that is facilitated by the catalytic Mg2+ ion. The simulations suggest that Mg2+ is profoundly involved in the hammerhead ribozyme mechanism both at structural and catalytic levels. Binding of Mg2+ in the active site plays a key structural role in the stabilization of stem I and II and to facilitate formation of near attack conformations and interactions between the nucleophile and G12, the implicated general base catalyst. In the transition state, Mg2+ binds in a bridging position where it stabilizes the accumulated charge of the leaving group while interacting with the 2'OH of G8, the implicated general acid catalyst. The QM/MM simulations provide support that, in the late transition state, the 2'OH of G8 can transfer a proton to the leaving group while directly coordinating the bridging Mg2+ ion. The present study provides evidence for the role of Mg2+ in hammerhead ribozyme catalysis. The proposed simulation model reconciles the interpretation of available experimental structural and biochemical data, and provides a starting point for more detailed investigation of the chemical reaction path with combined QM/MM methods.     

Mobility of the active site bound paraoxon and sarin in zinc-phosphotriesterase by molecular dynamics simulation and quantum chemical calculation.

The kinetic data published on phosphotriesterase (PTE), with various complexed metals, clearly indicates that the P=O and P=S bonds of phosphotriester and thiophosphotriester substrates, respectively, are strongly polarized by one or both of the active site complexed metal ions. However, this observation is not consistent with the three-dimensional X-ray crystal structure of zinc-substituted PTE with active site bound substrate analogue diethyl 4-methylbenzylphosphonate. In this structure, the distance between the phosphoryl oxygen and the nearest zinc is 3.4 A, a distance too large to afford strong polarization. In the present paper, the geometry and mobility of various PTE active site-substrate complexes are examined by performing both molecular dynamics (MD) simulations and quantum mechanical calculations. Two known substrates are considered, paraoxon and sarin, although their turnover rates vary about 100-fold. The results indicate that PTE forms a complex with either substrate in which the phosphoryl oxygen becomes strongly coordinated with the less buried zinc atom. It is shown that the geometry of the active site is changed when the protein is immersed in a water bath and relaxed by MD. The most substantial conformational change is the opening of the gateway in a pocket where the location of the leaving group is expected. The opening is observed for the pure enzyme as well as for the enzyme/substrate complexes and it ranges from 11 to 18 A. It is also shown that the pockets, in which the substrate substituents are localized, exhibit different flexibility and interact with the substrate with coordinated conformational adjustments.