Combined QM/MM and path integral simulations of kinetic isotope effects in the proton transfer reaction between nitroethane and acetate ion in water

An integrated Feynman path integral-free energy perturbation and umbrella sampling (PI-FEP/UM) method has been used to investigate the kinetic isotope effects (KIEs) in the proton transfer reaction between nitroethane and acetate ion in water. In the present study, both nuclear and electronic quantum effects are explicitly treated for the reacting system. The nuclear quantum effects are represented by bisection sampling centroid path integral simulations, while the potential energy surface is described by a combined quantum mechanical and molecular mechanical (QM/MM) potential. The accuracy essential for computing KIEs is achieved by a FEP technique that transforms the mass of a light isotope into a heavy one, which is equivalent to the perturbation of the coordinates for the path integral quasiparticle in the bisection sampling scheme. The PI-FEP/UM method is applied to the proton abstraction of nitroethane by acetate ion in water through molecular dynamics simulations. The rule of the geometric mean and the Swain-Schaad exponents for various isotopic substitutions at the primary and secondary sites have been examined. The computed total deuterium KIEs are in accord with experiments. It is found that the mixed isotopic Swain-Schaad exponents are very close to the semiclassical limits, suggesting that tunneling effects do not significantly affect this property for the reaction between nitroethane and acetate ion in aqueous solution.     

A quantum generalization of intrinsic reaction coordinate using path integral centroid coordinates

We propose a generalization of the intrinsic reaction coordinate (IRC) for quantum many-body systems described in terms of the mass-weighted ring polymer centroids in the imaginary-time path integral theory. This novel kind of reaction coordinate, which may be called the "centroid IRC," corresponds to the minimum free energy path connecting reactant and product states with a least amount of reversible work applied to the center of masses of the quantum nuclei, i.e., the centroids. We provide a numerical procedure to obtain the centroid IRC based on first principles by combining ab initio path integral simulation with the string method. This approach is applied to NH(3) molecule and N(2)H(5) (-) ion as well as their deuterated isotopomers to study the importance of nuclear quantum effects in the intramolecular and intermolecular proton transfer reactions. We find that, in the intramolecular proton transfer (inversion) of NH(3), the free energy barrier for the centroid variables decreases with an amount of about 20% compared to the classical one at the room temperature. In the intermolecular proton transfer of N(2)H(5) (-), the centroid IRC is largely deviated from the "classical" IRC, and the free energy barrier is reduced by the quantum effects even more drastically.     

Transmission coefficient calculation for proton transfer in triosephosphate isomerase based on the reaction path potential method

A global potential energy surface has been constructed through interpolation of our recently developed reaction path potential for chemical reactions in enzymes which is derived from combined ab initio quantum mechanical and molecular mechanical calculations. It has been implemented for the activated molecular dynamics simulations of the initial proton transfer reaction catalyzed by triosephosphate isomerase. To examine the dynamical effects on the rate constants of the enzymatic reaction, the classical transmission coefficient kappa(t) is evaluated to be 0.47 with the reactive flux approach, demonstrating considerable deviations from transition state theory. In addition, the fluctuations of protein environments have small effects on the barrier recrossing, and the transmission coefficient kappa(t) strongly depends on the fluctuations of atoms near the active site of the enzyme.     

Simulation of tunneling in enzyme catalysis by combining a biased propagation approach and the quantum classical path method: application to lipoxygenase

The ability of using wave function propagation approaches to simulate isotope effects in enzymes is explored, focusing on the large H/D kinetic isotope effect of soybean lipoxygenase-1 (SLO-1). The H/D kinetic isotope effect (KIE) is calculated as the ratio of the rate constants for hydrogen and deuterium transfer. The rate constants are calculated from the time course of the H and D nuclear wave functions. The propagations are done using one-dimensional proton potentials generated as sections from the full multidimensional surface of the reacting system in the protein. The sections are obtained during a classical empirical valence bond (EVB) molecular dynamics simulation of SLO-1. Since the propagations require an extremely long time for treating realistic activation barriers, it is essential to use an effective biasing approach. Thus, we develop here an approach that uses the classical quantum path (QCP) method to evaluate the quantum free energy change associated with the biasing potential. This approach provides an interesting alternative to full QCP simulations and to other current approaches for simulating isotope effects in proteins. In particular, this approach can be used to evaluate the quantum mechanical transmission factor or other dynamical effects, while still obtaining reliable quantized activation free energies due to the QCP correction.     

Simulation of slow reaction with quantum character: Neutral hydrolysis of carboxylic ester

By computer simulation, using both quantum and classical dynamics, we determined the rate constant and the kinetic isotope effect of the rate-determining step in the neutral hydrolysis of p-methoxyphenyl dichloroacetate in aqueous solution. This step involves a proton transfer concerted with the formation of a C O bond. A method of biased sampling was used; the Gibbs free energy of the biased configuration from which proton transfer is likely to occur was determined by a combination of semiempirical quantum calculations and thermodynamic integration. The proton dynamics was modeled with the quantum-dynamical density matrix evolution method that includes nonadiabatic pathways. The proton dynamics is driven by a fluctuating proton potential that was derived from a classical molecular dynamics simulation of the system including solvent. The calculated rate constant of 3×10⁻² s⁻¹ agrees within the error of the calculation with the experimentally observed value of 2.78×10⁻³. The calculated pseudo-first-order kinetic isotope effect of 3.9 is in good agreement with the experimentally observed value of 3.2. The results show the feasibility of computational approaches to slow reactions in complex environments, where proton transfer with an essential quantum-dynamical nature is the rate-limiting step. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 886–895, 1999     

Hydride transfer catalyzed by xylose isomerase: mechanism and quantum effects

We have applied molecular dynamics umbrella-sampling simulation and ensemble-averaged variational transition state theory with multidimensional tunneling (EA-VTST/MT) to calculate the reaction rate of xylose-to- xylulose isomerization catalyzed by xylose isomerase in the presence of two Mg2+ ions. The calculations include determination of the free energy of activation profile and ensemble averaging in the transmission coefficient. The potential energy function is approximated by a combined QM/MM/SVB method involving PM3 for the quantum mechanical (QM) subsystem, CHARMM22 and TIP3P for the molecular mechanical (MM) environment, and a simple valence bond (SVB) local function of two bond distances for the hydride transfer reaction. The simulation confirms the essential features of a mechanism postulated on the basis of kinetics and X-ray data by Whitlow et al. (Whitlow, M.; Howard, A. J.; Finzel, B. C.; Poulos, T. L.; Winborne, E.; Gilliland, G. L. Proteins 1991, 9, 153) and Ringe, Petsko, and coworkers (Labie, A.; Allen, K.-N.; Petsko, G. A.; Ringe, D. Biochemistry 1994, 33, 5469). This mechanism involves a rate-determining 1,2-hydride shift with prior and post proton transfers. Inclusion of quantum mechanical vibrational energy is important for computing the free energy of activation, and quantum mechanical tunneling effects are essential for computing kinetic isotope effects (KIEs). It is found that 85% of the reaction proceeds by tunneling and 15% by overbarrier events. The computed KIE for the ratio of hydride to deuteride transfer is in good agreement with the experimental results. The molecular dynamics simulations reveal that proton and hydride transfer reactions are assisted by breathing motions of the mobile Mg2+ ion in the active site, providing evidence for concerted motion of Mg2+ during the hydride transfer step.     

A theoretical study on excited state double proton transfer reaction of a 7-azaindole dimer: an ab initio potential energy surface and its empirical valence bond model

Double proton transfer (DPT) reaction of a 7-azaindole dimer in the first ππ* electronically excited state was studied theoretically. We investigated the reaction mechanism through constructing a full dimensional empirical valence bond potential energy function (PEF) based on potential energies evaluated by ab initio molecular orbital methods, and carrying out quantum dynamics calculations with the PEF. Potential energy surfaces of the DPT obtained at the multi-reference perturbation level of theory favors a concerted DPT mechanism, although a stepwise channel is suggested to open for an excited initial vibrational state. Reduced two dimensional quantum dynamics calculations for a reaction surface Hamiltonian of DPT coordinates were performed. Time constants of the reaction were evaluated to be on the order of picoseconds, which is consistent with experiments. On the other hand, the computed kinetic isotope effect deviates from experimental evidence, suggesting the importance of intermolecular stretching motion, which is not explicit in the present calculations for the quantum effect.     

Hybrid quantum/classical molecular dynamics for a proton transfer reaction coupled to a dissipative bath

A hybrid quantum/classical molecular dynamics approach is applied to a proton transfer reaction represented by a symmetric double well system coupled to a dissipative bath. In this approach, the proton is treated quantum mechanically and all bath modes are treated classically. The transition state theory rate constant is obtained from the potential of mean force, which is generated along a collective reaction coordinate with umbrella sampling techniques. The transmission coefficient, which accounts for dynamical recrossings of the dividing surface, is calculated with a reactive flux approach combined with the molecular dynamics with quantum transitions surface hopping method. The hybrid quantum/classical results agree well with numerically exact results in the spatial-diffusion-controlled regime, which is most relevant for proton transfer in proteins. This hybrid quantum/classical approach has already been shown to be computationally practical for studying proton transfer in large biological systems. These results have important implications for future applications to hydrogen transfer reactions in solution and proteins.     

Solvent polarization and kinetic isotope effects in nitroethane deprotonation and implications to the nitroalkane oxidase reaction

We have carried out a mixed molecular dynamics and centroid path integral simulation using a combined quantum mechanical and molecular mechanical (QM/MM) potential to study the anomalous Br?nsted relationship between rates and equilibria for deprotonation of nitroalkanes in water, which is known as the nitroalkane anomaly. The deprotonation process is catalyzed by nitroalkane oxidase. Our results show that the difference in solvent polarization effects for the TS and products is a major factor for the differential solvent effects on rate and equilibrium of nitroalkane deprotonation. This is due to poor charge delocalization as a result of slow rehybridization compared to bond breaking. Although solvent effects do not affect significantly the computed kinetic isotope effects in comparison with the gas-phase value, there is slight solvent-induced increase in tunneling. The present results suggest that an effective means by which the transition state can be stabilized in the enzyme nitroalkane oxidase is to facilitate the Calpha rehybridization.     

Car-Parrinello and path integral molecular dynamics study of the hydrogen bond in the chloroacetic acid dimer system

We have studied the double proton transfer (DPT) reaction in the cyclic dimer of chloroacetic acid using both classical and path integral Car-Parrinello molecular dynamics. We also attempt to quantify the errors in the potential energy surface that arise from the use of a pure density functional. In the classical dynamics a clear reaction mechanism can be identified, where asynchronized DPT arises due to coupling between the O-H stretching oscillator and several low energy intermolecular vibrational modes. This mechanism is considerably altered when quantum tunneling is permitted in the simulation. The introduction of path integrals leads to considerable changes in the thermally averaged molecular geometry, leading to shorter and more centered hydrogen bond linkages.     

Simulations of the large kinetic isotope effect and the temperature dependence of the hydrogen atom transfer in lipoxygenase

Elucidating the role of nuclear quantum mechanical (NQM) effects in enzyme catalysis is a topic of significant current interest. Despite the great experimental progress in this field it is important to have theoretical approaches capable of evaluating and analyzing nuclear quantum mechanical contributions to catalysis. In this study, we use the catalytic reaction of lipoxygenase, which is characterized by an extremely large kinetic isotope effect, as a challenging test case for our simulation approach. This is done by applying the quantum classical path (QCP) method with an empirical valence bond potential energy surface. Our computational strategy evaluates the relevant NQM corrections and reproduces the large observed kinetic isotope effect and the temperature dependence of the H atom transfer reaction while being less successful with the D atom transfer reaction. However, the main point of our study is not so much to explore the temperature dependence of the isotope effect but rather to develop and validate an approach for calculations of nuclear quantum mechanical contributions to activation free energies. Here, we find that the deviation between the calculated and observed activation free energies is small for both H and D at all investigated temperatures. The present study also explores the nature of the reorganization energy in the enzyme and solution reactions. It is found that the outer-sphere reorganization energy is extremely small. This reflects the fact that the considered reaction involves a very small charge transfer. The implication of this finding is discussed in the framework of the qualitative vibronic model. The main point of the present study is, however, that the rigorous QCP approach provides a reliable computational tool for evaluating NQM contributions to catalysis even when the given reaction includes large tunneling contributions. Interestingly, our results indicate that the NQM effects in the lipoxygenase reaction are similar in the enzyme and in the reference solution reactions, and thus do not contribute to catalysis. We also reached similar conclusions in studies of other enzymes.     

Quantum dynamics of hydride transfer catalyzed by bimetallic electrophilic catalysis: synchronous motion of Mg(2+) and H(-) in xylose isomerase

Xylose isomerase exhibits a bridged-bimetallic active-site motif in which the substrate is bound to two metals connected by a glutamate bridge, and X-ray crystallographic studies suggest that metal movement is involved in the hydride transfer rate-controlling catalytic step. Here we report classical/quantal dynamical simulations of this step that provide new insight into the metal motion. The potential energy surface is calculated by treating xylose with semiempirical molecular orbital theory augmented by a simple valence bond potential and the rest of the system by molecular mechanics. The rate constant for the hydride-transfer step was calculated by ensemble-averaged dynamical simulations including both variational transition-state theory for determination of the statistically averaged dynamical bottleneck and optimized multidimensional tunneling calculations. The dynamics calculations include 25 317 atoms, with quantized vibrational free energy in 89 active-site degrees of freedom, and with 32 atoms moving through static secondary zone transition-state configurations in the quantum tunneling simulation. Our simulations show that the average Mg-Mg distance R increases monotonically as a function of the hydride-transfer progress variable z. The range of the average R along the reaction path is consistent with the X-ray structure, thus providing a dynamical demonstration of the postulated role of Mg in catalysis. We also predicted the primary deuterium kinetic isotope effect (KIE) for the chemical step. We calculated a KIE of 3.8 for xylose at 298 K, which is consistent with somewhat smaller experimentally observed KIEs for glucose substrate at higher temperatures. More than half of our KIE is due to tunneling; neglecting quantum effects on the reaction coordinate reduces the calculated KIE to 1.8.     

Efficient estimators for quantum instanton evaluation of the kinetic isotope effects: application to the intramolecular hydrogen transfer in pentadiene

The quantum instanton approximation is used to compute kinetic isotope effects for intramolecular hydrogen transfer in cis-1,3-pentadiene. Due to the importance of skeleton motions, this system with 13 atoms is a simple prototype for hydrogen transfer in enzymatic reactions. The calculation is carried out using thermodynamic integration with respect to the mass of the isotopes and a path integral Monte Carlo evaluation of relevant thermodynamic quantities. Efficient "virial" estimators are derived for the logarithmic derivatives of the partition function and the delta-delta correlation functions. These estimators require significantly fewer Monte Carlo samples since their statistical error does not increase with the number of discrete time slices in the path integral. The calculation treats all 39 degrees of freedom quantum mechanically and uses an empirical valence bond potential based on a molecular mechanics force field.     

Proton-coupled electron transfer in soybean lipoxygenase

The proton-coupled electron transfer reaction catalyzed by soybean lipoxygenase-1 is studied with a multistate continuum theory that represents the transferring hydrogen nucleus as a quantum mechanical wave function. The inner-sphere reorganization energy of the iron cofactor is calculated with density functional theory, and the outer-sphere reorganization energy of the protein is calculated with the frequency-resolved cavity model for conformations obtained with docking simulations. Both classical and quantum mechanical treatments of the proton donor-acceptor vibrational motion are presented. The temperature dependence of the calculated rates and kinetic isotope effects is in agreement with the experimental data. The weak temperature dependence of the rates is due to the relatively small free energy barrier arising from a balance between the reorganization energy and the reaction free energy. The unusually high deuterium kinetic isotope effect of 81 is due to the small overlap of the reactant and product proton vibrational wave functions and the dominance of the lowest energy reactant and product vibronic states in the tunneling process. The temperature dependence of the kinetic isotope effect is strongly influenced by the proton donor-acceptor distance with the dominant contribution to the overall rate. This dominant proton donor-acceptor distance is significantly smaller than the equilibrium donor-acceptor distance and is determined by a balance between the larger coupling and the smaller Boltzmann probability as the distance decreases. Thus, the proton donor-acceptor vibrational motion plays a vital role in decreasing the dominant donor-acceptor distance relative to its equilibrium value to facilitate the proton-coupled electron transfer reaction.     

Nonadiabatic proton/deuteron transfer within the benzophenone-triethylamine triplet contact radical ion pair: exploration of the influence of structure upon reaction

The dynamics of proton transfer within a variety of substituted benzophenone-triethylamine triplet contact radical ion pairs are examined in the solvents acetonitrile and dimethylformamide. The correlation of the proton-transfer rate constants with DeltaG reveals an inverted region. The kinetic deuterium isotope effects are also examined. The solvent and isotope dependence of the transfer processes are analyzed within the context of the Lee-Hynes model for nonadiabatic proton transfer. Theoretical analysis of the experimental data suggests that the reaction path for proton/deuteron transfer involves tunneling, and the origin of the inverted region is attributed to a curved tunneling path.     

Activation free energy for proton transfer in solution

Path integral molecular dynamics methods are employed to compute the free energy for proton transfer reactions for strongly hydrogen bonded systems in a polar solvent. The free energy profile is calculated using several different techniques, including: integration of the mean force acting on the proton path with its centroid constrained at different values, the integral form of the free energy calculation in the constrained-reaction-coordinate-dynamics ensemble and direct simulation of the unconstrained dynamics. The results show that estimates of the free energy barrier obtained by harmonic extrapolation are likely to be in error. Both quantum and classical results for the free energy are obtained and compared with simulations using adiabatic quantum dynamics. Comparison of the quantum and classical results show that there are quantum corrections to the solvent contributions to the free energy.     

A computational study of the intramolecular deprotonation of a carbon acid in aqueous solution

Proton transfer reactions are the rate-limiting steps in many biological and synthetic chemical processes, often requiring complex cofactors or catalysts to overcome the generally unfavourable thermodynamic process of carbanion intermediate formation. It has been suggested that quantum tunnelling processes enhance the kinetics of some of these reactions, which when coupled to protein motions may be an important consideration for enzyme catalysis. To obtain a better fundamental and quantitative understanding of these proton transfer mechanisms, a computational analysis of the intramolecular proton transfer from a carbon acid in the small molecule, 4-nitropentanoic acid, in aqueous solution is presented. Potential-energy surfaces from gas-phase, implicit and QM/MM (quantum mechanical/molecular mechanical) explicit solvation quantum chemistry models are compared, and the potential of mean force, for the full reaction coordinate, using umbrella-sampling molecular dynamics is analysed. Semi-classical multidimensional tunnelling corrections are also used to estimate the quantum tunnelling contributions and to understand the origin of the primary deuterium kinetic isotope effects (KIEs). The computational results are found to be in excellent agreement with the KIEs and the energetics obtained experimentally.     

Hybrid quantum/classical path integral approach for simulation of hydrogen transfer reactions in enzymes

A hybrid quantum/classical path integral Monte Carlo (QC-PIMC) method for calculating the quantum free energy barrier for hydrogen transfer reactions in condensed phases is presented. In this approach, the classical potential of mean force along a collective reaction coordinate is calculated using umbrella sampling techniques in conjunction with molecular dynamics trajectories propagated according to a mapping potential. The quantum contribution is determined for each configuration along the classical trajectory with path integral Monte Carlo calculations in which the beads move according to an effective mapping potential. This type of path integral calculation does not utilize the centroid constraint and can lead to more efficient sampling of the relevant region of conformational space than free-particle path integral sampling. The QC-PIMC method is computationally practical for large systems because the path integral sampling for the quantum nuclei is performed separately from the classical molecular dynamics sampling of the entire system. The utility of the QC-PIMC method is illustrated by an application to hydride transfer in the enzyme dihydrofolate reductase. A comparison of this method to the quantized classical path and grid-based methods for this system is presented.     

A combined quantum mechanical and molecular mechanical study of the reaction mechanism and alpha-amino acidity in alanine racemase

Combined quantum mechanical/molecular mechanical simulations have been carried out to investigate the origin of the carbon acidity enhancement in the alanine racemization reaction catalyzed by alanine racemase (AlaR). The present study shows that the enhancement of carbon acidity of alpha-amino acids by the cofactor pyridoxal 5'-phosphate (PLP) with an unusual, unprotonated pyridine is mainly due to solvation effects, in contrast to the intrinsic electron-withdrawing stabilization by the pyridinium ion to form a quinonoid intermediate. Alanine racemase further lowers the alpha-proton acidity and provides an overall 14-17 kcal/mol transition-state stabilization. The second key finding of this study is that the mechanism of racemization of an alanine zwitterion in water is altered from an essentially concerted process to a stepwise reaction by formation of an external aldimine adduct with the PLP cofactor. Finally, we have used a centroid path integral method to determine the intrinsic kinetic isotope effects for the two proton abstraction reactions, which are somewhat greater than the experimental estimates.     

Proton Transfer Studied Using a Combined Ab Initio Reactive Potential Energy Surface with Quantum Path Integral Methodology

The rates of intramolecular proton transfer are calculated on a full-dimensional reactive electronic potential energy surface that incorporates high level ab initio calculations along the reaction path and by using classical Transition State theory, Path-Integral Quantum Transition State Theory, and the Quantum Instanton approach. The specific example problem studied is malonaldehyde. Estimates of the kinetic isotope effect using the latter two methods are found to be in reasonable agreement with each other. Improvements and extensions of this practical, yet chemically accurate framework for the calculations of quantized, reactive dynamics are also discussed.