Quantum mechanics and mixed quantum mechanics/molecular mechanics simulations of model nerve agents with acetylcholinesterase

 The accurate modeling of biological processes presents major computational difficulties owing to the inherent complexity of the macromolecular systems of interest. Simulations of biochemical reactivity tend to require highly computationally intensive quantum mechanical methods, but localized chemical effects tend to depend significantly on properties of the extended biological environment – a regime far more readily examined with lower-level classical empirical models. Mixed quantum/classical techniques are gaining in popularity as a means of bridging these competing requirements. Here we present results comparing two quantum mechanics/molecular mechanics implementations (the SIMOMM technique of Gordon et al. as implemented in GAMESS, and the ONIOM technique of Morokuma et al. found in Gaussian 98) as performed on the enzyme acetylcholinesterase and model nerve agents. This work represents part of the initial phase of a DoD HPCMP Challenge project in which we are attempting to reliably characterize the biochemical processes responsible for nerve agent activity and inhibition, thereby allowing predictions on compounds unrelated to those already studied. Received: 10 October 2001 / Accepted: 13 November 2002 / Published online: 1 April 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: M. M. Hurley e-mail: hurley@arl.army.mil     

Cytochrome P450CAM enzymatic catalysis cycle: a quantum mechanics/molecular mechanics study

The catalytic pathway of cytochrome P450cam is studied by means of a hybrid quantum mechanics/molecular mechanics method. Our results reveal an active role of the enzyme in the different catalytic steps. The protein initially controls the energy gap between the high- and low-spin states in the substrate binding process, allowing thermodynamic reduction by putidaredoxin reductase and molecular oxygen addition. A second electron reduction activates the delivery of protons to the active site through a selective interaction of Thr252 and the distal oxygen causing the O--O cleavage. Finally, the protein environment catalyzes the substrate hydrogen atom abstraction step with a remarkably low free energy barrier ( approximately 8 kcal/mol). Our results are consistent with the effect of mutations on the enzymatic efficacy and provide a satisfactory explanation for the experimental failure to trap the proposed catalytically competent species, a ferryl Fe(IV) heme.     

Structure and ultrafast dynamics of liquid water: a quantum mechanics/molecular mechanics molecular dynamics simulations study

A quantum mechanics/molecular mechanics molecular dynamics simulation was performed for liquid water to investigate structural and dynamical properties of this peculiar liquid. The most important region containing a central reference molecule and all nearest surrounding molecules (first coordination shell) was treated by Hartree-Fock (HF), post-Hartree-Fock [second-order Moller-Plesset perturbation theory (MP2)], and hybrid density functional B3LYP [Becke's three parameter functional (B3) with the correlation functional of Lee, Yang, and Parr (LYP)] methods. In addition, another HF-level simulation (2HF) included the full second coordination shell. Site to site interactions between oxygen-oxygen, oxygen-hydrogen, and hydrogen-hydrogen atoms of all ab initio methods were compared to experimental data. The absence of a second peak and the appearance of a shoulder instead in the gO-O graph obtained from the 2HF simulation is notable, as this feature has been observed so far only for pressurized or heated water. Dynamical data show that the 2HF procedure compensates some of the deficiency of the HF one-shell simulation, reducing the difference between correlated (MP2) and HF results. B3LYP apparently leads to too rigid structures and thus to an artificial slow down of the dynamics.     

MoBioTools: A toolkit to setup quantum mechanics/molecular mechanics calculations

We present a toolkit that allows for the preparation of QM/MM input files from a conformational ensemble of molecular geometries. The package is currently compatible with trajectory and topology files in Amber, CHARMM, GROMACS and NAMD formats, and has the possibility to generate QM/MM input files for Gaussian (09 and 16), Orca (≥4.0), NWChem and (Open)Molcas. The toolkit can be used in command line, so that no programming experience is required, although it presents some features that can also be employed as a python application programming interface. We apply the toolkit in four situations in which different electronic-structure properties of organic molecules in the presence of a solvent or a complex biological environment are computed: the reduction potential of the nucleobases in acetonitrile, an energy decomposition analysis of tyrosine interacting with water, the absorption spectrum of an azobenzene derivative integrated into a voltage-gated ion channel, and the absorption and emission spectra of the luciferine/luciferase complex. These examples show that the toolkit can be employed in a manifold of situations for both the electronic ground state and electronically excited states. It also allows for the automatic correction of the active space in the case of CASSCF calculations on an ensemble of geometries, as it is shown for the azobenzene derivative photoswitch case.     

Catalytic mechanism of RNA backbone cleavage by ribonuclease H from quantum mechanics/molecular mechanics simulations

We use quantum mechanics/molecular mechanics simulations to study the cleavage of the ribonucleic acid (RNA) backbone catalyzed by ribonuclease H. This protein is a prototypical member of a large family of enzymes that use two-metal catalysis to process nucleic acids. By combining Hamiltonian replica exchange with a finite-temperature string method, we calculate the free energy surface underlying the RNA-cleavage reaction and characterize its mechanism. We find that the reaction proceeds in two steps. In a first step, catalyzed primarily by magnesium ion A and its ligands, a water molecule attacks the scissile phosphate. Consistent with thiol-substitution experiments, a water proton is transferred to the downstream phosphate group. The transient phosphorane formed as a result of this nucleophilic attack decays by breaking the bond between the phosphate and the ribose oxygen. In the resulting intermediate, the dissociated but unprotonated leaving group forms an alkoxide coordinated to magnesium ion B. In a second step, the reaction is completed by protonation of the leaving group, with a neutral Asp132 as a likely proton donor. The overall reaction barrier of ～15 kcal mol(-1), encountered in the first step, together with the cost of protonating Asp132, is consistent with the slow measured rate of ～1-100/min. The two-step mechanism is also consistent with the bell-shaped pH dependence of the reaction rate. The nonmonotonic relative motion of the magnesium ions along the reaction pathway agrees with X-ray crystal structures. Proton-transfer reactions and changes in the metal ion coordination emerge as central factors in the RNA-cleavage reaction.     

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.     

Protein/ligand binding free energies calculated with quantum mechanics/molecular mechanics

The calculation of binding affinities for flexible ligands has hitherto required the availability of reliable molecular mechanics parameters for the ligands, a restriction that can in principle be lifted by using a mixed quantum mechanics/molecular mechanics (QM/MM) representation in which the ligand is treated quantum mechanically. The feasibility of this approach is evaluated here, combining QM/MM with the Poisson-Boltzmann/surface area model of continuum solvation and testing the method on a set of 47 benzamidine derivatives binding to trypsin. The experimental range of the absolute binding energy (DeltaG = -3.9 to -7.6 kcal/mol) is reproduced well, with a root-mean-square (RMS) error of 1.2 kcal/mol. When QM/MM is applied without reoptimization to the very different ligands of FK506 binding protein the RMS error is only 0.7 kcal/mol. The results show that QM/MM is a promising new avenue for automated docking and scoring of flexible ligands. Suggestions are made for further improvements in accuracy.     

Hybrid quantum mechanics/molecular mechanics-based molecular dynamics simulation of acid-catalyzed dehydration of polyols in liquid water

We use the conversion of protonated glycerol to acrolein for a case study of the mechanism of acid-catalyzed dehydration of polyols in aqueous environments. We employ hybrid Quamtum Mechanics/Molecular Mechanics Molecular Dynamics (QM/MM MD) simulations with biased sampling and perform free energy calculations for the elementary steps of the reaction. We investigate the effects of solvent dynamics and in particular the role of quantum mechanical water in the dehydration mechanism. We present results supporting a mechanism that proceeds via water-mediated proton transfers and thus through an enol intermediate. We find that the first dehydration may take place by two, low-energy pathways requiring, respectively, 20.9 and 18.8 kcal/mol of activation free energy. The second dehydration requires 19.9 kcal/mol of activation free energy while for the overall reaction we compute a free energy change of -8 kcal/mol.     

Mechanisms of enzymatic hydrolysis of nucleoside triphosphates by quantum and molecular mechanics

Results of simulation of the mechanism of hydrolysis of adenosine triphosphate and guanosine triphosphate in protein matrices, as well as of deprotonated methyl triphosphate in water clusters by quantum and molecular mechanics with separation of the reaction system into conformationally flexible effective fragments are discussed.     

Proton storage site in bacteriorhodopsin: new insights from quantum mechanics/molecular mechanics simulations of microscopic pK(a) and infrared spectra

Identifying the group that acts as the proton storage/loading site is a challenging but important problem for understanding the mechanism of proton pumping in biomolecular proton pumps, such as bacteriorhodopsin (bR) and cytochrome c oxidase. Recent experimental studies of bR propelled the idea that the proton storage/release group (PRG) in bR is not an amino acid but a water cluster embedded in the protein. We argue that this idea is at odds with our knowledge of protein electrostatics, since invoking the water cluster as the PRG would require the protein to raise the pK(a) of a hydronium by almost 11 pK(a) units, which is difficult considering known cases of pK(a) shifts in proteins. Our recent quantum mechanics/molecular mechanics (QM/MM) simulations suggested an alternative "intermolecular proton bond" model in which the stored proton is shared between two conserved Glu residues (194 and 204). Here we show that this model leads to microscopic pK(a) values consistent with available experimental data and the functional requirement of a PRG. Extensive QM/MM simulations also show that, independent of a number of technical issues, such as the influence of QM region size, starting X-ray structure, and nuclear quantum effects, the "intermolecular proton bond" model is qualitatively consistent with available spectroscopic data. Potential of mean force calculations show explicitly that the stored proton strongly prefers the pair of Glu residues over the water cluster. The results and analyses help highlight the importance of considering protein electrostatics and provide arguments for why the "intermolecular proton bond" model is likely applicable to the PRG in biomolecular proton pumps in general.     

Density functional self-consistent quantum mechanics/molecular mechanics theory for linear and nonlinear molecular properties: Applications to solvated water and formaldehyde

A combined quantum mechanics/molecular mechanics (QM/MM) method is described, where the polarization between the solvent and solute is accounted for using a self-consistent scheme linear in the solvent polarization. The QM/MM method is implemented for calculation of energies and molecular response properties including the calculation of linear and quadratic response functions using the density-functional theory (DFT) and the Hartree-Fock (HF) theory. Sample calculations presented for ground-state energies, first-order ground-state properties, excitation energies, first-order excited state properties, polarizabilities, first-hyperpolarizabilities, and two-photon absorptions strengths of formaldehyde suggests that DFT may in some cases be a sufficiently reliable alternative to high-level theory, such as coupled-cluster (CC) theory, in modeling solvent shifts, whereas results obtained with the HF wave function deviate significantly from the CC results. Calculations carried out on water gives results that also are comparable with CC calculations in accuracy for ground-state and first-order properties. However, to obtain such accuracy an exchange-correlation functional capable of describing the diffuse Rydberg states must be chosen.     

Thermal Gaussian molecular dynamics for quantum dynamics simulations of many-body systems: application to liquid para-hydrogen

A new method, here called thermal Gaussian molecular dynamics (TGMD), for simulating the dynamics of quantum many-body systems has recently been introduced [I. Georgescu and V. A. Mandelshtam, Phys. Rev. B 82, 094305 (2010)]. As in the centroid molecular dynamics (CMD), in TGMD the N-body quantum system is mapped to an N-body classical system. The associated both effective Hamiltonian and effective force are computed within the variational Gaussian wave-packet approximation. The TGMD is exact for the high-temperature limit, accurate for short times, and preserves the quantum canonical distribution. For a harmonic potential and any form of operator A?, it provides exact time correlation functions C(AB)(t) at least for the case of B, a linear combination of the position, x, and momentum, p, operators. While conceptually similar to CMD and other quantum molecular dynamics approaches, the great advantage of TGMD is its computational efficiency. We introduce the many-body implementation and demonstrate it on the benchmark problem of calculating the velocity time auto-correlation function for liquid para-hydrogen, using a system of up to N = 2592 particles.     

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

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

Structural forms of green fluorescent protein by quantum mechanics/molecular mechanics calculations

The molecular modeling of structural forms of the green fluorescent protein (GFP) with the Ser65Thr single-site mutation was performed by the quantum mechanics/molecular mechanics (QM/MM) method. Two model systems were constructed based on the crystallographic structure from the Protein Data Bank (PDB entry code 1EMA.) The model systems differ in the initial protonation state of the side chain of the amino acid residue Glu222 near the chromophore. The atomic coordinates of the protein macromolecule corresponding to the equilibrium geometric configurations were determined by total energy minimization using the QM/MM method within the density functional theory approximation PBE0/cc-pVDZ for the quantum subsystem that consists of the chromophore, a water molecule, and the side chains of Arg96, Glu222, and Ser205, and with the parameters of the AMBER force field for the molecular mechanics subsystem. In the analysis of the results, particular attention was given to the hydrogen bond redistribution in the chromophore-containing region of the protein caused by a change in the protonation state of the chromophore. The results obtained from the model containing the initially protonated side chain of Glu222 suggest a new interpretation of the photophysical processes in the green fluorescent protein.     

Substrate hydroxylation in methane monooxygenase: quantitative modeling via mixed quantum mechanics/molecular mechanics techniques

Using broken-symmetry unrestricted density functional theory quantum mechanical (QM) methods in concert with mixed quantum mechanics/molecular mechanics (QM/MM) methods, the hydroxylation of methane and substituted methanes by intermediate Q in methane monooxygenase hydroxylase (MMOH) has been quantitatively modeled. This protocol allows the protein environment to be included throughout the calculations and its effects (electrostatic, van der Waals, strain) upon the reaction to be accurately evaluated. With the current results, recent kinetic data for CH3X (X = H, CH3, OH, CN, NO2) substrate hydroxylation in MMOH (Ambundo, E. A.; Friesner, R. A.; Lippard, S. J. J. Am. Chem. Soc. 2002, 124, 8770-8771) can be rationalized. Results for methane, which provide a quantitative test of the protocol, including a substantial kinetic isotope effect (KIE), are in reasonable agreement with experiment. Specific features of the interaction of each of the substrates with MMO are illuminated by the QM/MM modeling, and the resulting effects upon substrate binding are quantitatively incorporated into the calculations. The results as a whole point to the success of the QM/MM methodology and enhance our understanding of MMOH catalytic chemistry. We also identify systematic errors in the evaluation of the free energy of binding of the Michaelis complexes of the substrates, which most likely arise from inadequate sampling and/or the use of harmonic approximations to evaluate the entropy of the complex. More sophisticated sampling methods will be required to achieve greater accuracy in this aspect of the calculation.     

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.     

GridMol2.0: Implementation and application of linear-scale quantum mechanics methods and molecular visualization

GridMol is a “one-stop” platform for molecular structure building, scientific computing, and molecular visualization aided by a high-performance computing environment. GridMol version 2.0 introduces two unique features: the first is fragment-based linear-scaling quantum chemistry methods, such as molecular fractionation with conjugate caps and fragment molecular orbital methods; the second is that GridMol enables users to visualize molecular geometries along a geometry optimization and an intrinsic reaction coordinate calculation. Compared with version 1.0, fragment-based linear-scaling quantum chemistry methods implemented in GridMol version 2.0 can be used as a useful tool for performing quantum calculations for large molecular systems to explore the mechanisms involved in protein-ligand or targeted drug interactions.